function [P, L, U, pivcol, sign] = plu(A) % plu Rectangular PA=LU factorization *with row exchanges*. % % [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 PA = LU. % U is the same size as A. % P and L are square, with as many rows as A. % sign = det(P); it is 1 or -1. % % See also elim, slu, lu, rref, partic, nulbasis, colbasis. [m, n] = size(A); P = eye(m, m); L = eye(m, m); U = zeros(m, n); pivcol = []; tol = sqrt(eps); 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 end