Operation Count:

n-1 divisions to compute the numbers a_i1/a₁₁; n(n-1) multiplications and n(n-1) subtractions to do the elimination. Total: approximately 2 n² operations. [Note: We are interested in the operation count for large n. Hence it will not matter much whether it is 2 n² or 2n(n-1). The coefficient of the largest term n² is still 2.]

MATLAB code for this step (after rows have been exchanged so that the (1,1) element is nonzero):

% Eliminate first column using first row.
for i=2:n,
  mult = A(i,1)/A(1,1);
  for k=1:n,
    A(i,k) = A(i,k) - mult*A(1,k);
  end;
  b(i) = b(i) - mult*b(1);
end;

Note that the inner for loop can be replaced by the statement:

 
A(i,:) = A(i,:) - mult*A(1,:);   % This does the same thing as the inner 
                                 % loop above.

Step 2. If [a\tilde]₂₂ = 0, exchange rows so that the (2,2) element is nonzero. For each row i = 3, ¼, n, replace row i by itself minus [a\tilde]_i2/ [a\tilde]₂₂ times row 2.

Operation Count: n-2 divisions to compute the numbers [a\tilde]_i2/ [a\tilde]₂₂; (n-1)(n-2) multiplications and (n-1)(n-2) subtractions to do the elimination. Total: approximately 2 (n-1)² operations.

MATLAB code for this step (after rows have been exchanged so that the (2,2) element is nonzero):

for i=3:n,
  mult = A(i,2)/A(2,2);
  A(i,2:n) = A(i,2:n) - mult*A(2,2:n);  % Note that we work only with
                                        % columns 2 through n.
end;

Step 3. If [[a\tilde]\tilde]₃₃ = 0, exchange rows so that the (3,3) element is nonzero. For each row i = 4, ¼, n, replace row i by itself minus [[a\tilde]\tilde]_i3/ [[a\tilde]\tilde]₃₃ times row 3.

Operation Count: 2( n-2 )².

Etc.

Total Operation Count for Gaussian Elimination:

MATLAB Code:

%  Given an n by n nonsingular matrix A and an n-vector b, does Gaussian
%  elimination (with pivoting) to bring A to upper triangular form.

for j=1:n-1,
  [cmax,jmax] = max(abs(A(j:n,j))); % Find index of maximum element in column j.
  temp = A(j,:);                    % Exchange rows j and j+jmax-1.
  A(j,:) = A(j+jmax-1,:); 
  A(j+jmax-1,:) = temp;
  tempb = b(j);                   % Also exchange elements j and j+jmax-1 of b.
  b(j) = b(j+jmax-1);
  b(j+jmax-1) = tempb;
  for i=j+1:n,                      % Eliminate in column j using row j.
    mult = A(i,j)/A(j,j);
    A(i,j:n) = A(i,j:n) - mult*A(j,j:n);
    b(i) = b(i) - mult*b(j);
  end;
end;

The matrix is now upper triangular. We still must backsolve to find the solution of the linear system:

For i = n,n-1, ¼, 1,

The MATLAB code for backsolving is:

% Solve the linear system Ax=b once A has been reduced to upper triangular
% form.
x = zeros(n,1);      % Force x to be a column vector.
x(n) = b(n)/A(n,n);
for i=n-1:-1:1,
  x(i) = (b(i) - A(i,i+1:n)*x(i+1:n))/A(i,i);
end;