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	Thu Sep 26 08:30:26 2013 -0700
+++ b/doc/interpreter/linalg.txi	Thu Sep 26 09:38:51 2013 -0700
@@ -96,6 +96,8 @@
 
 @DOCSTRING(inv)
 
+@DOCSTRING(linsolve)
+
 @DOCSTRING(matrix_type)
 
 @DOCSTRING(norm)
--- a/scripts/linear-algebra/linsolve.m	Thu Sep 26 08:30:26 2013 -0700
+++ b/scripts/linear-algebra/linsolve.m	Thu Sep 26 09:38:51 2013 -0700
@@ -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)