octave_packages/m/linear-algebra/rref.m

   1 ## Copyright (C) 2000-2012 Paul Kienzle
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {Function File} {} rref (@var{A})
  21 ## @deftypefnx {Function File} {} rref (@var{A}, @var{tol})
  22 ## @deftypefnx {Function File} {[@var{r}, @var{k}] =} rref (@dots{})
  23 ## Return the reduced row echelon form of @var{A}.  @var{tol} defaults
  24 ## to @code{eps * max (size (@var{A})) * norm (@var{A}, inf)}.
  25 ##
  26 ## Called with two return arguments, @var{k} returns the vector of
  27 ## "bound variables", which are those columns on which elimination
  28 ## has been performed.
  29 ##
  30 ## @end deftypefn
  31
  32 ## Author: Paul Kienzle <pkienzle@users.sf.net>
  33 ##         (based on an anonymous source from the public domain)
  34
  35 function [A, k] = rref (A, tol)
  36
  37   if (nargin < 1 || nargin > 2)
  38     print_usage ();
  39   endif
  40
  41   if (ndims (A) > 2)
  42     error ("rref: expecting matrix argument");
  43   endif
  44
  45   [rows, cols] = size (A);
  46
  47   if (nargin < 2)
  48     if (isa (A, "single"))
  49       tol = eps ("single") * max (rows, cols) * norm (A, inf ("single"));
  50     else
  51       tol = eps * max (rows, cols) * norm (A, inf);
  52     endif
  53   endif
  54
  55   used = zeros (1, cols);
  56   r = 1;
  57   for c = 1:cols
  58     ## Find the pivot row
  59     [m, pivot] = max (abs (A(r:rows,c)));
  60     pivot = r + pivot - 1;
  61
  62     if (m <= tol)
  63       ## Skip column c, making sure the approximately zero terms are
  64       ## actually zero.
  65       A (r:rows, c) = zeros (rows-r+1, 1);
  66     else
  67       ## keep track of bound variables
  68       used (1, c) = 1;
  69
  70       ## Swap current row and pivot row
  71       A ([pivot, r], c:cols) = A ([r, pivot], c:cols);
  72
  73       ## Normalize pivot row
  74       A (r, c:cols) = A (r, c:cols) / A (r, c);
  75
  76       ## Eliminate the current column
  77       ridx = [1:r-1, r+1:rows];
  78       A (ridx, c:cols) = A (ridx, c:cols) - A (ridx, c) * A(r, c:cols);
  79
  80       ## Check if done
  81       if (r++ == rows)
  82         break;
  83       endif
  84     endif
  85   endfor
  86   k = find (used);
  87
  88 endfunction
  89
  90 %!test
  91 %! a = [1];
  92 %! [r k] = rref(a);
  93 %! assert(r, [1], 2e-8);
  94 %! assert(k, [1], 2e-8);
  95
  96 %!test
  97 %! a = [1 3; 4 5];
  98 %! [r k] = rref(a);
  99 %! assert(rank(a), rank(r), 2e-8);
 100 %! assert(r, eye(2), 2e-8);
 101 %! assert(k == [1, 2] || k == [2, 1]);
 102
 103
 104 %!test
 105 %! a = [1 3; 4 5; 7 9];
 106 %! [r k] = rref(a);
 107 %! assert(rank(a), rank(r), 2e-8);
 108 %! assert(r, eye(3)(:,1:2), 2e-8);
 109 %! assert(k, [1 2], 2e-8);
 110
 111 %!test
 112 %! a = [1 2 3; 2 4 6; 7 2 0];
 113 %! [r k] = rref(a);
 114 %! assert(rank(a), rank(r), 2e-8);
 115 %! assert(r, [1 0 (3-7/2); 0 1 (7/4); 0 0 0], 2e-8);
 116 %! assert(k, [1 2], 2e-8);
 117
 118 %!test
 119 %! a = [1 2 1; 2 4 2.01; 2 4 2.1];
 120 %! tol = 0.02;
 121 %! [r k] = rref(a, tol);
 122 %! assert(rank(a, tol), rank(r, tol), 2e-8);
 123 %! tol = 0.2;
 124 %! [r k] = rref(a, tol);
 125 %! assert(rank(a, tol), rank(r, tol), 2e-8);
 126
 127 %!error rref();
 128