function [L,U,P] = LUfacP(A);
   %
   %  Input:   A  an n by n matrix. (A is unaltered on return)
   %       
   %           A  is first overwritten with the LU decompositon
   %              L in strict lower triangle, U in upper triangle
   %              and then L,U are explicitly constructed such that
   %              
   %                       PA = LU
   %
   %  Output:  L  a unit lower triangular n by n matrix
   %
   %           U  an n by n upper triangular matrix
   %
   %           P  an n by n permutation matrix
   %
   %
   %
   %
   %  D.C. Sorensen 
   %  3 Oct 00
   %-----------------------------------------------------------------
   %     
      [n,n] = size(A);
      p = zeros(n,1);
      for k = 1:n-1,
   %
   %    Find index of max elt. below diagonal, adjust for offset
   %    and  record in p(k)    
   %
          [tk,ik] = max(abs(A(k:n,k))); p(k) = ik+k-1;
   %
   %    Interchange row k with row p(k)
   %
          t = A(k,k:n); A(k,k:n) = A(p(k),k:n); A(p(k),k:n) = t;
   %
   %    Scale by 1/A(k,k) to determine the kth column of L
   %
          A(k+1:n,k) = A(k+1:n,k)/A(k,k);
   %
   %    Update lower (n-k) by (n-k) block
   %
   
          for j = k+1:n,
              for i = k+1:n,
                     A(i,j) = A(i,j) - A(i,k)*A(k,j);
              end
          end
      end
   %
   %  Apply permutations successively to the identity to get P
   %  and to the multipliers (lower triangle of overwritten A)     
   %  Note: the k-th permutation is applied to columns 1 : k-1.
   %
      P = eye(n);
      t = P(1,:); P(1,:) = P(p(1),:); P(p(1),:) = t;

      for k = 2:n-1;
          t = A(k,1:k-1); A(k,1:k-1) = A(p(k),1:k-1); A(p(k),1:k-1) = t;
          t = P(k,:); P(k,:) = P(p(k),:); P(p(k),:) = t;
      end 
      L = eye(n) + tril(A,-1);
      U = triu(A);