function [A,pvt,errflag] = ge_pp(A)
%
% Performs Gaussian elimination with partial pivoting on a 
% square matrix A. 
%
% Input: 
%  A     = square matrix
%
% Outputs: 
%     A  = matrix with the upper triangular factor of A in the 
%          upper triangular part and the multipliers in strict lower 
%          triangular part. 
%   pvt  = vector with integer entries containing pivot information. 
%
% errflag = error flag: 
%           0 => elimination completed with no error
%           k>0 => elimination stopped with singularity indicated 
%                  by zero pivot at step k
%

n = size(A,1);
if n ~= size(A,2), error('\n The matrix must be square.\n'); end
pvt = [1:1:n-1]';

maxA = max(max(abs(A)));

errflag = 0; 

% Elimination loop. 
for k = 1:n-1,
% Determine the kth pivot element.
  [max_elt,pvt_row] = max(abs(A(k:n,k)));
  if maxA\max_elt < 100*eps,
     fprintf('\n Elimination stopped with near-singularity indicated at step %d.\n', k);
     errflag = k;
     break; 
  end
  pvt_row = pvt_row + k - 1;
% Swap rows if necessary. 
  if pvt_row > k,
    pvt(k,1) = pvt_row;
    temp = A(k,k:n); 
    A(k,k:n) = A(pvt_row,k:n);
    A(pvt_row,k:n) = temp;
  end
% Form the kth column of multipliers. 
  A(k+1:n,k) = -A(k,k)\A(k+1:n,k);
% Perform the kth elimination step. 
  A(k+1:n,k+1:n) = A(k+1:n,k+1:n) + A(k+1:n,k)*A(k,k+1:n);
end

if maxA\abs(A(n,n)) < 100*eps, 
     fprintf('\n Elimination completed but near-singularity indicated at step %d.\n', n);
     errflag = n;
end

    
