Mercurial > octave
view scripts/linear-algebra/rref.m @ 33580:80346999b171 bytecode-interpreter tip
build: Fix typo in test/compile/module.mk (bug #65658).
* test/compile/module.mk: Fix typo in path to file.
author | A.R. Burgers <arburgers@gmail.com> |
---|---|
date | Mon, 13 May 2024 11:33:36 +0200 |
parents | 2e484f9f1f18 |
children |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2000-2024 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## This file is part of Octave. ## ## Octave is free software: you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{r} =} rref (@var{A}) ## @deftypefnx {} {@var{r} =} rref (@var{A}, @var{tol}) ## @deftypefnx {} {[@var{r}, @var{k}] =} rref (@dots{}) ## Return the reduced row echelon form of @var{A}. ## ## @var{tol} defaults to ## @code{eps * max (size (@var{A})) * norm (@var{A}, inf)}. ## ## The optional return argument @var{k} contains the vector of ## "bound variables", which are those columns on which elimination has been ## performed. ## ## @end deftypefn function [A, k] = rref (A, tol) if (nargin < 1) print_usage (); endif if (ndims (A) > 2) error ("rref: A must be a 2-dimensional matrix"); endif [rows, cols] = size (A); if (nargin < 2) if (isa (A, "single")) tol = eps ("single") * max (rows, cols) * norm (A, inf ("single")); else tol = eps * max (rows, cols) * norm (A, inf); endif endif used = zeros (1, cols); r = 1; for c = 1:cols ## Find the pivot row [m, pivot] = max (abs (A(r:rows,c))); pivot = r + pivot - 1; if (m <= tol) ## Skip column c, making sure the approximately zero terms are ## actually zero. A(r:rows, c) = zeros (rows-r+1, 1); else ## keep track of bound variables used(1, c) = 1; ## Swap current row and pivot row A([pivot, r], c:cols) = A([r, pivot], c:cols); ## Normalize pivot row A(r, c:cols) = A(r, c:cols) / A(r, c); ## Eliminate the current column ridx = [1:r-1, r+1:rows]; A(ridx, c:cols) = A(ridx, c:cols) - A(ridx, c) * A(r, c:cols); ## Check if done if (r++ == rows) break; endif endif endfor if (nargout > 1) k = find (used); endif endfunction %!test %! a = [1]; %! [r k] = rref (a); %! assert (r, [1], 2e-8); %! assert (k, [1], 2e-8); %!test %! a = [1 3; 4 5]; %! [r k] = rref (a); %! assert (rank (a), rank (r), 2e-8); %! assert (r, eye (2), 2e-8); %! assert (k == [1, 2] || k == [2, 1]); %!test %! a = [1 3; 4 5; 7 9]; %! [r k] = rref (a); %! assert (rank (a), rank (r), 2e-8); %! assert (r, eye (3)(:,1:2), 2e-8); %! assert (k, [1 2], 2e-8); %!test %! a = [1 2 3; 2 4 6; 7 2 0]; %! [r k] = rref (a); %! assert (rank (a), rank (r), 2e-8); %! assert (r, [1 0 (3-7/2); 0 1 (7/4); 0 0 0], 2e-8); %! assert (k, [1 2], 2e-8); %!test %! a = [1 2 1; 2 4 2.01; 2 4 2.1]; %! tol = 0.02; %! [r k] = rref (a, tol); %! assert (rank (a, tol), rank (r, tol), 2e-8); %! tol = 0.2; %! [r k] = rref (a, tol); %! assert (rank (a, tol), rank (r, tol), 2e-8); ## Test input validation %!error <Invalid call> rref ()