Mercurial > octave
view scripts/linear-algebra/linsolve.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) 2013-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{x} =} linsolve (@var{A}, @var{b}) ## @deftypefnx {} {@var{x} =} linsolve (@var{A}, @var{b}, @var{opts}) ## @deftypefnx {} {[@var{x}, @var{R}] =} linsolve (@dots{}) ## Solve the linear system @code{A*x = b}. ## ## With no options, this function is equivalent to the left division operator ## @w{(@code{x = A \ b})} or the matrix-left-divide function ## @w{(@code{x = mldivide (A, b)})}. ## ## Octave ordinarily examines the properties of the matrix @var{A} and chooses ## a solver that best matches the matrix. By passing a structure @var{opts} ## to @code{linsolve} you can inform Octave directly about the matrix @var{A}. ## In this case Octave will skip the matrix examination and proceed directly ## to solving the linear system. ## ## @strong{Warning:} If the matrix @var{A} does not have the properties listed ## in the @var{opts} structure then the result will not be accurate AND no ## warning will be given. When in doubt, let Octave examine the matrix and ## choose the appropriate solver as this step takes little time and the result ## is cached so that it is only done once per linear system. ## ## Possible @var{opts} fields (set value to true/false): ## ## @table @asis ## @item LT ## @var{A} is lower triangular ## ## @item UT ## @var{A} is upper triangular ## ## @item UHESS ## @var{A} is upper Hessenberg (currently makes no difference) ## ## @item SYM ## @var{A} is symmetric or complex Hermitian (currently makes no difference) ## ## @item POSDEF ## @var{A} is positive definite ## ## @item RECT ## @var{A} is general rectangular (currently makes no difference) ## ## @item TRANSA ## Solve @code{A'*x = b} if true rather than @code{A*x = b} ## @end table ## ## The optional second output @var{R} is the inverse condition number of ## @var{A} (zero if matrix is singular). ## @seealso{mldivide, matrix_type, rcond} ## @end deftypefn function [x, R] = linsolve (A, b, opts) if (nargin < 2) print_usage (); endif if (! (isnumeric (A) && isnumeric (b))) error ("linsolve: A and B must be numeric"); endif trans_A = false; ## Process any opts if (nargin > 2) if (! isstruct (opts)) error ("linsolve: OPTS must be a structure"); endif if (isfield (opts, "TRANSA") && opts.TRANSA) trans_A = true; endif if (isfield (opts, "POSDEF") && opts.POSDEF) A = matrix_type (A, "positive definite"); endif if (isfield (opts, "LT") && opts.LT) A = matrix_type (A, "lower"); elseif (isfield (opts, "UT") && opts.UT) A = matrix_type (A, "upper"); endif endif ## This way is faster as the transpose is not calculated in Octave, ## but forwarded as a flag option to BLAS. if (trans_A) x = A' \ b; else x = A \ b; endif if (nargout > 1) if (issquare (A)) R = rcond (A); else R = 0; endif endif endfunction %!test %! n = 10; %! A = rand (n); %! x = rand (n, 1); %! b = A * x; %! assert (linsolve (A, b), A \ b); %! assert (linsolve (A, b, struct ()), A \ b); %!test %! n = 10; %! A = triu (gallery ("condex", n)); %! x = rand (n, 1); %! b = A' * x; %! opts.UT = true; %! opts.TRANSA = true; %! assert (linsolve (A, b, opts), A' \ b); %!error <Invalid call> linsolve () %!error <Invalid call> linsolve (1) %!error linsolve (1,2,3) %!error <A and B must be numeric> linsolve ({1},2) %!error <A and B must be numeric> linsolve (1,{2}) %!error <OPTS must be a structure> linsolve (1,2,3)