%  gauss_el.m   Gaussian elimination with partial pivoting

function [U,f] = gauss_el(A,b,display)

%  On entry: A is a square matrix,
%            b is a column vector of the same dimension,
%            display is 1 if step-by-step display desired, 0 otherwise.
%  On exit:  U and f are the result of applying Gaussian elimination
%            with partial pivoting to the system Ax=b.

n = length(b);
if size(A) ~= [n,n] | size(b) ~= [n,1],
  error('mismatched dimension')
end;

for j=1:n-1,
  [cmax,jmax] = max(abs(A(j:n,j)));  % Find index of max 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;
  if display,                        % Display results after each step
    A, b, pause,                     % if desired.
  end;
  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);
    if display,                      % Display results after each step 
      A, b, pause,                   % if desired.
    end;
  end;
end;
U = A;
f = b;