Mercurial > hg > octave-nkf
changeset 17500:b66f068e4468
linsolve.m: Use Octave coding conventions.
* doc/interpreter/linalg.txi: Add linsolve to manual.
* scripts/linear-algebra/linsolve.m: Redo docstring. Use Octave coding conventions.
Add %!error input validation tests.
author | Rik <rik@octave.org> |
---|---|
date | Thu, 26 Sep 2013 09:38:51 -0700 |
parents | c3a3532e3d98 |
children | 55680de6a897 |
files | doc/interpreter/linalg.txi scripts/linear-algebra/linsolve.m |
diffstat | 2 files changed, 102 insertions(+), 57 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/interpreter/linalg.txi +++ b/doc/interpreter/linalg.txi @@ -96,6 +96,8 @@ @DOCSTRING(inv) +@DOCSTRING(linsolve) + @DOCSTRING(matrix_type) @DOCSTRING(norm)
--- a/scripts/linear-algebra/linsolve.m +++ b/scripts/linear-algebra/linsolve.m @@ -14,79 +14,122 @@ ## along with this program; If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File}{[@var{X}, @var{R}] =} csaps(@var{A}, @var{B}, @var{options}=[]) +## @deftypefn {Function File} {@var{x} =} linsolve (@var{A}, @var{b}) +## @deftypefnx {Function File} {@var{x} =} linsolve (@var{A}, @var{b}, @var{opts}) +## @deftypefnx {Function File} {[@var{x}, @var{R}] =} linsolve (@dots{}) +## Solve the linear system @code{A*x = b}. ## -## Solve a linear system@* +## 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)})}. ## -## With no options, this is the same as @code{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. ## -## Possible option fields (set to true/false): +## @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 @var{LT} -## A is lower triangular -## @item @var{UT} -## A is upper triangular -## @item @var{UHESS} -## A is upper Hessenberg (currently makes no difference) -## @item @var{SYM} -## A is symmetric (currently makes no difference) -## @item @var{POSDEF} -## A is positive definite -## @item @var{RECT} -## A is general rectangular (currently makes no difference) -## @item @var{TRANSA} -## Compute @code{transpose(A) \ B} +## @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} by @code{transpose (A) \ b} ## @end table ## -## The optional second output @var{R} is the inverse condition number of @var{A} (zero if matrix is singular) +## 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 ## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> -function [X, R] = linsolve(A, B, options) - -trans_A = false; +function [x, R] = linsolve (A, b, opts) -#process any options -if nargin > 2 - if ~isstruct(options) - error('Third input must be a structure') + if (nargin < 2 || nargin > 3) + print_usage (); + endif + + if (! (isnumeric (A) && isnumeric (b))) + error ("linsolve: A and B must be numeric"); endif - if isfield(options, 'TRANSA') && options.TRANSA - trans_A = true; - A = A'; + + ## Process any opts + if (nargin > 2) + if (! isstruct (opts)) + error ("linsolve: OPTS must be a structure"); + endif + trans_A = false; + if (isfield (opts, "TRANSA") && opts.TRANSA) + trans_A = true; + A = A'; + endif + if (isfield (opts, "POSDEF") && opts.POSDEF) + A = matrix_type (A, "positive definite"); + endif + if (isfield (opts, "LT") && opts.LT) + if (trans_A) + A = matrix_type (A, "upper"); + else + A = matrix_type (A, "lower"); + endif + endif + if (isfield (opts, "UT") && opts.UT) + if (trans_A) + A = matrix_type (A, "lower"); + else + A = matrix_type (A, "upper"); + endif + endif endif - if isfield(options, 'POSDEF') && options.POSDEF - A = matrix_type (A, 'positive definite'); - endif - if isfield(options, 'LT') && options.LT - if trans_A - A = matrix_type (A, 'upper'); + + x = A \ b; + + if (nargout > 1) + if (issquare (A)) + R = rcond (A); else - A = matrix_type (A, 'lower'); + R = 0; endif endif - if isfield(options, 'UT') && options.UT - if trans_A - A = matrix_type (A, 'lower'); - else - A = matrix_type (A, 'upper'); - endif - endif -endif +endfunction -X = A \ B; -if nargout > 1 - if issquare(A) - R = 1 ./ cond(A); - else - R = 0; - endif -endif +%!test +%! n = 4; +%! A = triu (rand (n)); +%! x = rand (n, 1); +%! b = A' * x; +%! opts.UT = true; +%! opts.TRANSA = true; +%! assert (linsolve (A, b, opts), A' \ b); -%!shared n, A, B, x, opts -%! n = 4; A = triu(rand(n)); x = rand(n, 1); B = A' * x; -%! opts.UT = true; opts.TRANSA = true; -%!assert (linsolve(A, B, opts), A' \ B); - +%!error linsolve () +%!error 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)