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