function [p,A] = plu (A)
%
% function [p,A] = plu (A)
%
% Perform LU decomposition with scaled partial pivoting.
% Upon return the coefficients of L and U replace those 
% of the input n-by-n nonsingular matrix A. The row interchanges
% performed are recorded in the vector p.
%

n = size(A,1);

% find scales, initialize permutation vector p

s = max(abs(A'))';
p = 1:n;

% LU decomposition with partial pivoting

for k = 1:n-1

   % find row index of relative maximum in column k
   [val,q] = max ( abs(A(k:n,k)) ./ s(k:n) );
   q = q + k-1;

   % interchange rows k and q and record this in p
   t = A(k,1:n); ts = s(k); tp = p(k);
   A(k,1:n) = A(q,1:n); s(k) = s(q); p(k) = p(q);
   A(q,1:n) = t; s(q) = ts; p(q) = tp;

   % Compute the corresponding column of L
   A(k+1:n,k) = A(k+1:n,k) / A(k,k);

   % Update submatrix by outer product
   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
