Mercurial > octave-antonio
annotate scripts/linear-algebra/rref.m @ 10549:95c3e38098bf
Untabify .m scripts
| author | Rik <code@nomad.inbox5.com> |
|---|---|
| date | Fri, 23 Apr 2010 11:28:50 -0700 |
| parents | 16f53d29049f |
| children | c776f063fefe |
| rev | line source |
|---|---|
| 9245 | 1 ## Copyright (C) 2000, 2006, 2007, 2008, 2009 Paul Kienzle |
| 5827 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 5827 | 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 | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 5827 | 18 |
| 19 ## -*- texinfo -*- | |
| 20 ## @deftypefn {Function File} {[@var{r}, @var{k}] =} rref (@var{a}, @var{tol}) | |
| 21 ## | |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
22 ## Returns the reduced row echelon form of @var{a}. @var{tol} defaults |
| 5827 | 23 ## to @code{eps * max (size (@var{a})) * norm (@var{a}, inf)}. |
| 24 ## | |
| 25 ## Called with two return arguments, @var{k} returns the vector of | |
| 26 ## "bound variables", which are those columns on which elimination | |
| 27 ## has been performed. | |
| 28 ## | |
| 29 ## @end deftypefn | |
| 30 | |
| 31 ## Author: Paul Kienzle <pkienzle@users.sf.net> | |
| 32 ## (based on a anonymous source from the public domain) | |
| 33 | |
| 34 function [A, k] = rref (A, tolerance) | |
| 35 | |
| 36 if (nargin < 1 || nargin > 2) | |
| 37 print_usage (); | |
| 38 endif | |
| 39 | |
| 40 if (ndims (A) > 2) | |
| 41 error ("rref: expecting matrix argument"); | |
| 42 endif | |
| 43 | |
| 44 [rows, cols] = size (A); | |
| 45 | |
| 46 if (nargin < 2) | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
47 if (isa (A, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
48 tolerance = eps ("single") * max (rows, cols) * norm (A, inf ("single")); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
49 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
50 tolerance = eps * max (rows, cols) * norm (A, inf); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
51 endif |
| 5827 | 52 endif |
| 53 | |
| 54 used = zeros (1, cols); | |
| 55 r = 1; | |
| 56 for c = 1:cols | |
| 57 ## Find the pivot row | |
| 58 [m, pivot] = max (abs (A(r:rows,c))); | |
| 59 pivot = r + pivot - 1; | |
| 60 | |
| 61 if (m <= tolerance) | |
| 62 ## Skip column c, making sure the approximately zero terms are | |
| 63 ## actually zero. | |
| 64 A (r:rows, c) = zeros (rows-r+1, 1); | |
| 65 else | |
| 66 ## keep track of bound variables | |
| 67 used (1, c) = 1; | |
| 68 | |
| 69 ## Swap current row and pivot row | |
| 70 A ([pivot, r], c:cols) = A ([r, pivot], c:cols); | |
| 71 | |
| 72 ## Normalize pivot row | |
| 73 A (r, c:cols) = A (r, c:cols) / A (r, c); | |
| 74 | |
| 75 ## Eliminate the current column | |
| 76 ridx = [1:r-1, r+1:rows]; | |
| 77 A (ridx, c:cols) = A (ridx, c:cols) - A (ridx, c) * A(r, c:cols); | |
| 78 | |
| 79 ## Check if done | |
| 80 if (r++ == rows) | |
| 10549 | 81 break; |
| 5827 | 82 endif |
| 83 endif | |
| 84 endfor | |
| 85 k = find (used); | |
| 86 | |
| 87 endfunction |