  function [L,U,piv] = GEpiv(A);
%
% Pre:
%   A      n-by-n
%
% Post:
%   L      n-by-n unit lower triangular with |L(i,j)|<=1.
%   U      n-by-n upper triangular
%   piv    integer n-vector that is a permutation of 1:n.
%
%          A(piv,:) = LU  
 
    [n,n] = size(A); 
	piv = 1:n;
    for k=1:n-1
       [maxv,r] = max(abs(A(k:n,k)));
       q = r+k-1;
	   piv([k q]) = piv([q k]);
	   A([k q],:) = A([q k],:);
       if A(k,k) ~= 0
          A(k+1:n,k) = A(k+1:n,k)/A(k,k);
          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
    end
	L = eye(n,n) + tril(A,-1);
    U = triu(A);

         
