PLU.M
上传用户:sfyaiting
上传日期:2009-10-25
资源大小:320k
文件大小:1k
源码类别:

GPS编程

开发平台:

Matlab

  1. function [P,L,U,pivcol,sign] = plu(A)
  2. %PLU Pivoting, rectangular, LU factorization.
  3. % [P,L,U] = PLU(A), for a rectangular matrix A, uses Gaussian elimination
  4. % to compute a permutation matrix P, a lower triangular matrix L and 
  5. % an upper trapezoidal matrix U so that L*U = P*A.
  6. % U is the same size as A.  P and L are square, with as many rows as A.
  7. %
  8. % See also SLU, LU, REF, SOLVE, NULL, BASIS.
  9. [m,n] = size(A);
  10. P = eye(m,m);
  11. L = eye(m,m);
  12. U = zeros(m,n);
  13. pivcol = [];
  14. tol = 1.e-6;
  15. sign = 1;
  16. p = 1;
  17. for k = 1:min(m,n)
  18.    [r, p] = findpiv(A(k:m,p:n),k,p,tol);
  19.    if r ~= k
  20.       A([r k],1:n) = A([k r],1:n);
  21.       if k > 1
  22.          L([r k],1:k-1) = L([k r],1:k-1); 
  23.       end
  24.       P([r k],1:m) = P([k r],1:m);
  25.       sign = -sign;
  26.    end
  27.    if abs(A(k,p)) >= tol
  28.       pivcol = [pivcol p];
  29.       for i = k+1:m
  30.          L(i,k) = A(i,p)/A(k,p);
  31.          for j = k+1:n
  32.             A(i,j) = A(i,j) - L(i,k)*A(k,j);
  33.          end
  34.       end
  35.    end
  36.    for j = k:n
  37.       U(k,j) = A(k,j) * (abs(A(k,j)) >= tol);
  38.    end
  39.    if p < n
  40.       p = p+1;
  41.    end
  42. end
  43. if nargout < 4
  44.    nopiv = 1:n;
  45.    nopiv(pivcol) = [];
  46.    if ~isempty(pivcol)
  47.       disp('Pivots in columns:'), disp(pivcol);
  48.    end
  49.    if ~isempty(nopiv)
  50.       disp('No pivots in columns:'), disp(nopiv);
  51.    end
  52.    rank = length(pivcol);
  53.    if rank > 0
  54.       roworder = P*(1:m)';
  55.       disp('Pivots in rows:'), disp(roworder(1:rank)'); 
  56.    end
  57. end