1######################################################################## 2## 3## Copyright (C) 2013-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{x} =} linsolve (@var{A}, @var{b}) 28## @deftypefnx {} {@var{x} =} linsolve (@var{A}, @var{b}, @var{opts}) 29## @deftypefnx {} {[@var{x}, @var{R}] =} linsolve (@dots{}) 30## Solve the linear system @code{A*x = b}. 31## 32## With no options, this function is equivalent to the left division operator 33## @w{(@code{x = A \ b})} or the matrix-left-divide function 34## @w{(@code{x = mldivide (A, b)})}. 35## 36## Octave ordinarily examines the properties of the matrix @var{A} and chooses 37## a solver that best matches the matrix. By passing a structure @var{opts} 38## to @code{linsolve} you can inform Octave directly about the matrix @var{A}. 39## In this case Octave will skip the matrix examination and proceed directly 40## to solving the linear system. 41## 42## @strong{Warning:} If the matrix @var{A} does not have the properties listed 43## in the @var{opts} structure then the result will not be accurate AND no 44## warning will be given. When in doubt, let Octave examine the matrix and 45## choose the appropriate solver as this step takes little time and the result 46## is cached so that it is only done once per linear system. 47## 48## Possible @var{opts} fields (set value to true/false): 49## 50## @table @asis 51## @item LT 52## @var{A} is lower triangular 53## 54## @item UT 55## @var{A} is upper triangular 56## 57## @item UHESS 58## @var{A} is upper Hessenberg (currently makes no difference) 59## 60## @item SYM 61## @var{A} is symmetric or complex Hermitian (currently makes no difference) 62## 63## @item POSDEF 64## @var{A} is positive definite 65## 66## @item RECT 67## @var{A} is general rectangular (currently makes no difference) 68## 69## @item TRANSA 70## Solve @code{A'*x = b} if true rather than @code{A*x = b} 71## @end table 72## 73## The optional second output @var{R} is the inverse condition number of 74## @var{A} (zero if matrix is singular). 75## @seealso{mldivide, matrix_type, rcond} 76## @end deftypefn 77 78function [x, R] = linsolve (A, b, opts) 79 80 if (nargin < 2 || nargin > 3) 81 print_usage (); 82 endif 83 84 if (! (isnumeric (A) && isnumeric (b))) 85 error ("linsolve: A and B must be numeric"); 86 endif 87 88 trans_A = false; 89 90 ## Process any opts 91 if (nargin > 2) 92 if (! isstruct (opts)) 93 error ("linsolve: OPTS must be a structure"); 94 endif 95 if (isfield (opts, "TRANSA") && opts.TRANSA) 96 trans_A = true; 97 endif 98 if (isfield (opts, "POSDEF") && opts.POSDEF) 99 A = matrix_type (A, "positive definite"); 100 endif 101 if (isfield (opts, "LT") && opts.LT) 102 A = matrix_type (A, "lower"); 103 elseif (isfield (opts, "UT") && opts.UT) 104 A = matrix_type (A, "upper"); 105 endif 106 endif 107 108 ## This way is faster as the transpose is not calculated in Octave, 109 ## but forwarded as a flag option to BLAS. 110 if (trans_A) 111 x = A' \ b; 112 else 113 x = A \ b; 114 endif 115 116 if (nargout > 1) 117 if (issquare (A)) 118 R = rcond (A); 119 else 120 R = 0; 121 endif 122 endif 123 124endfunction 125 126 127%!test 128%! n = 10; 129%! A = rand (n); 130%! x = rand (n, 1); 131%! b = A * x; 132%! assert (linsolve (A, b), A \ b); 133%! assert (linsolve (A, b, struct ()), A \ b); 134 135%!test 136%! n = 10; 137%! A = triu (gallery ("condex", n)); 138%! x = rand (n, 1); 139%! b = A' * x; 140%! opts.UT = true; 141%! opts.TRANSA = true; 142%! assert (linsolve (A, b, opts), A' \ b); 143 144%!error linsolve () 145%!error linsolve (1) 146%!error linsolve (1,2,3) 147%!error <A and B must be numeric> linsolve ({1},2) 148%!error <A and B must be numeric> linsolve (1,{2}) 149%!error <OPTS must be a structure> linsolve (1,2,3) 150