splu.m

function [P, L, U, sign] = splu(A)

% splu  Square PA=LU factorization *with row exchanges*.
%
% [P, L, U] = splu(A), for a square, invertible matrix A,
% uses Gaussian elimination to compute a permutation
% matrix P, a lower triangular matrix L and 
% an upper triangular matrix U so that P*A = L*U.
% P, L and U are the same size as A.
% sign = det(P); it is 1 or -1.
%
% See also slu, lu, rref, partic, nulbasis, colbasis.

[m, n] = size(A);
if m ~= n
   error('Matrix must be square.')
end
P = eye(n, n);
L = eye(n, n);
U = zeros(n, n);
tol = sqrt(eps);
sign = 1;

for k = 1:n
   if abs(A(k, k)) < tol
      for r = k:n
         if abs(A(r, k)) >= tol
            break
         end
         if r == n
            if nargout == 4
               sign = 0;
               return
            else
               disp('A is singular within tolerance.')
               error(['No pivot in column ' int2str(k) '.'])
            end
         end
      end
      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:n) = P([k r], 1:n);
      sign = -sign;
   end
   for i = k+1:n
      L(i, k) = A(i, k) / A(k, k);
      for j = k+1:n
         A(i, j) = A(i, j) - L(i, k)*A(k, j);
      end
   end
   for j = k:n
      U(k, j) = A(k, j) * (abs(A(k, j)) >= tol);
   end
end

if nargout < 4
   roworder = P*(1:n)';
   disp('Pivots in rows:'), disp(roworder'); end
end