home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / MATFUN.DI$ / RREF.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.7 KB  |  59 lines
  1. function [A,jb] = rref(A,tol)
  2. %RREF    Reduced row echelon form.
  3. %    R = RREF(A) produces the reduced row echelon form of A.
  4. %
  5. %    [R,jb] = RREF(A) also returns a vector, jb, so that:
  6. %        r = length(jb) is this algorithm's idea of the rank of A,
  7. %        x(jb) are the bound variables in a linear system, Ax = b,
  8. %        A(:,jb) is a basis for the range of A,
  9. %        R(1:r,jb) is the r-by-r identity matrix.
  10. %
  11. %    [R,jb] = RREF(A,TOL) uses the given tolerance in the rank tests.
  12. %
  13. %    See also RREFMOVIE, RANK, ORTH, NULL, QR, SVD.
  14.  
  15. %    CBM, 11/24/85, 1/29/90, 7/12/92.
  16. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  17.  
  18. [m,n] = size(A);
  19.  
  20. % Does it appear that elements of A are ratios of small integers?
  21. [num, den] = rat(A);
  22. rats = all(all(A==num./den));
  23.  
  24. % Compute the default tolerance if none was provided.
  25. if (nargin < 2), tol = max(m,n)*eps*norm(A,'inf'); end
  26.  
  27. % Loop over the entire matrix.
  28. i = 1;
  29. j = 1;
  30. jb = [];
  31. while (i <= m) & (j <= n)
  32.    % Find value and index of largest element in the remainder of column j.
  33.    [p,k] = max(abs(A(i:m,j))); k = k+i-1;
  34.    if (p <= tol)
  35.       % The column is negligible, zero it out.
  36.       A(i:m,j) = zeros(m-i+1,1);
  37.       j = j + 1;
  38.    else
  39.       % Remember column index
  40.       jb = [jb j];
  41.       % Swap i-th and k-th rows.
  42.       A([i k],j:n) = A([k i],j:n);
  43.       % Divide the pivot row by the pivot element.
  44.       A(i,j:n) = A(i,j:n)/A(i,j);
  45.       % Subtract multiples of the pivot row from all the other rows.
  46.       for k = [1:i-1 i+1:m]
  47.          A(k,j:n) = A(k,j:n) - A(k,j)*A(i,j:n);
  48.       end
  49.       i = i + 1;
  50.       j = j + 1;
  51.    end
  52. end
  53.  
  54. % Return "rational" numbers if appropriate.
  55. if rats
  56.     [num,den] = rat(A);
  57.     A=num./den;
  58. end
  59.