plu.m

function [P,L,U,pivcol,sign] = plu(A)
%PLU	Pivoting, rectangular, LU factorization.
%	[P,L,U] = PLU(A), for a rectangular matrix A, uses Gaussian elimination
%	to compute a permutation matrix P, a lower triangular matrix L and 
%	an upper trapezoidal matrix U so that L*U = P*A.
%	U is the same size as A.  P and L are square, with as many rows as A.
%
%	See also SLU, LU, REF, SOLVE, NULL, BASIS.

[m,n] = size(A);
P = eye(m,m);
L = eye(m,m);
U = zeros(m,n);
pivcol = [];
tol = 1.e-6;
sign = 1;

p = 1;
for k = 1:min(m,n)
   [r, p] = findpiv(A(k:m,p:n),k,p,tol);
   if r ~= k
      A([r k],1:n) = A([k r],1:n);
      if k > 1
         L([r k],1:k-1) = L([k r],1:k-1); 
      end
      P([r k],1:m) = P([k r],1:m);
      sign = -sign;
   end
   if abs(A(k,p)) >= tol
      pivcol = [pivcol p];
      for i = k+1:m
         L(i,k) = A(i,p)/A(k,p);
         for j = k+1:n
            A(i,j) = A(i,j) - L(i,k)*A(k,j);
         end
      end
   end
   for j = k:n
      U(k,j) = A(k,j) * (abs(A(k,j)) >= tol);
   end
   if p < n
      p = p+1;
   end
end

if nargout < 4
   nopiv = 1:n;
   nopiv(pivcol) = [];
   if ~isempty(pivcol)
      disp('Pivots in columns:'), disp(pivcol);
   end
   if ~isempty(nopiv)
      disp('No pivots in columns:'), disp(nopiv);
   end
   rank = length(pivcol);
   if rank > 0
      roworder = P*(1:m)';
      disp('Pivots in rows:'), disp(roworder(1:rank)'); 
   end
end