function [A, ipivt, iflag] = lu_pp( A );
% 
% LU_PP  computes the LU--decomposition with partial
% pivoting of a matrix A 
% 
% Usage
%       [A, ipivt, iflag] = lu_pp( A )
% 
% input:
%        A:     the n by n matrix A
% 
% output:
%        A:     the LU-decomposition of  A
% 
%        ipivt: pivot information
% 
%        iflag: error flag
%               iflag = 0  Row reductions could be performed,
%                          A is upper triangular
%               iflag = 1  dimension of A or of b is not correct
% 
%
% Matthias Heinkenschloss
% Feb 24, 2001


iflag = 0;

% get size of A and check dimensions
[m,n]    = size(A);
if ( m ~= n )
   iflag = 1;
   return
end

% Start LU--decomposition

for k = 1:n-1
%  Find pivot index
   [amax, i ] = max(abs(A(k:n,k)));
   i          = i + k - 1;
   ipivt(k)   = i;
     

%  Interchange rows if necessary
   if( k ~= i ) 
       tmp      = A(k,k:n);
       A(k,k:n) = A(i,k:n);
       A(i,k:n) = tmp;
   end 

   if( amax > 0 )  % if amax == 0 subcolumn A(k:n,k) is zero
       for i = k+1:n
%          compute multipliers
           A(i,k) = -A(i,k)/A(k,k);

%          Perform row elimination
           A(i,k+1:n) = A(i,k+1:n) + A(i,k) * A(k,k+1:n);
       end
   end
end

%end of lu_pp.m