Mercurial > octave-nkf

comparison scripts/linear-algebra/linsolve.m @ 17500:b66f068e4468

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 0a8c35ae5ce1
comparison
equal deleted inserted replaced
17499:c3a3532e3d98 17500:b66f068e4468
12 ## 12 ##
13 ## You should have received a copy of the GNU General Public License 13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; If not, see <http://www.gnu.org/licenses/>. 14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
15 15
16 ## -*- texinfo -*- 16 ## -*- texinfo -*-
17 ## @deftypefn{Function File}{[@var{X}, @var{R}] =} csaps(@var{A}, @var{B}, @var{options}=[]) 17 ## @deftypefn {Function File} {@var{x} =} linsolve (@var{A}, @var{b})
18 ## @deftypefnx {Function File} {@var{x} =} linsolve (@var{A}, @var{b}, @var{opts})
19 ## @deftypefnx {Function File} {[@var{x}, @var{R}] =} linsolve (@dots{})
20 ## Solve the linear system @code{A*x = b}.
18 ## 21 ##
19 ## Solve a linear system@* 22 ## With no options, this function is equivalent to the left division operator
23 ## @w{(@code{x = A \ b})} or the matrix-left-divide function
24 ## @w{(@code{x = mldivide (A, b)})}.
20 ## 25 ##
21 ## With no options, this is the same as @code{A \ B} 26 ## Octave ordinarily examines the properties of the matrix @var{A} and chooses
27 ## a solver that best matches the matrix. By passing a structure @var{opts}
28 ## to @code{linsolve} you can inform Octave directly about the matrix @var{A}.
29 ## In this case Octave will skip the matrix examination and proceed directly
30 ## to solving the linear system.
22 ## 31 ##
23 ## Possible option fields (set to true/false): 32 ## @strong{Warning:} If the matrix @var{A} does not have the properties
33 ## listed in the @var{opts} structure then the result will not be accurate
34 ## AND no warning will be given. When in doubt, let Octave examine the matrix
35 ## and choose the appropriate solver as this step takes little time and the
36 ## result is cached so that it is only done once per linear system.
37 ##
38 ## Possible @var{opts} fields (set value to true/false):
39 ##
24 ## @table @asis 40 ## @table @asis
25 ## @item @var{LT} 41 ## @item LT
26 ## A is lower triangular 42 ## @var{A} is lower triangular
27 ## @item @var{UT} 43 ##
28 ## A is upper triangular 44 ## @item UT
29 ## @item @var{UHESS} 45 ## @var{A} is upper triangular
30 ## A is upper Hessenberg (currently makes no difference) 46 ##
31 ## @item @var{SYM} 47 ## @item UHESS
32 ## A is symmetric (currently makes no difference) 48 ## @var{A} is upper Hessenberg (currently makes no difference)
33 ## @item @var{POSDEF} 49 ##
34 ## A is positive definite 50 ## @item SYM
35 ## @item @var{RECT} 51 ## @var{A} is symmetric or complex Hermitian (currently makes no difference)
36 ## A is general rectangular (currently makes no difference) 52 ##
37 ## @item @var{TRANSA} 53 ## @item POSDEF
38 ## Compute @code{transpose(A) \ B} 54 ## @var{A} is positive definite
55 ##
56 ## @item RECT
57 ## @var{A} is general rectangular (currently makes no difference)
58 ##
59 ## @item TRANSA
60 ## Solve @code{A'*x = b} by @code{transpose (A) \ b}
39 ## @end table 61 ## @end table
40 ## 62 ##
41 ## The optional second output @var{R} is the inverse condition number of @var{A} (zero if matrix is singular) 63 ## The optional second output @var{R} is the inverse condition number of
64 ## @var{A} (zero if matrix is singular).
65 ## @seealso{mldivide, matrix_type, rcond}
42 ## @end deftypefn 66 ## @end deftypefn
43 67
44 ## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> 68 ## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu>
45 69
46 function [X, R] = linsolve(A, B, options) 70 function [x, R] = linsolve (A, b, opts)
47 71
48 trans_A = false; 72 if (nargin < 2 || nargin > 3)
73 print_usage ();
74 endif
49 75
50 #process any options 76 if (! (isnumeric (A) && isnumeric (b)))
51 if nargin > 2 77 error ("linsolve: A and B must be numeric");
52 if ~isstruct(options)
53 error('Third input must be a structure')
54 endif 78 endif
55 if isfield(options, 'TRANSA') && options.TRANSA 79
56 trans_A = true; 80 ## Process any opts
57 A = A'; 81 if (nargin > 2)
82 if (! isstruct (opts))
83 error ("linsolve: OPTS must be a structure");
84 endif
85 trans_A = false;
86 if (isfield (opts, "TRANSA") && opts.TRANSA)
87 trans_A = true;
88 A = A';
89 endif
90 if (isfield (opts, "POSDEF") && opts.POSDEF)
91 A = matrix_type (A, "positive definite");
92 endif
93 if (isfield (opts, "LT") && opts.LT)
94 if (trans_A)
95 A = matrix_type (A, "upper");
96 else
97 A = matrix_type (A, "lower");
98 endif
99 endif
100 if (isfield (opts, "UT") && opts.UT)
101 if (trans_A)
102 A = matrix_type (A, "lower");
103 else
104 A = matrix_type (A, "upper");
105 endif
106 endif
58 endif 107 endif
59 if isfield(options, 'POSDEF') && options.POSDEF 108
60 A = matrix_type (A, 'positive definite'); 109 x = A \ b;
61 endif 110
62 if isfield(options, 'LT') && options.LT 111 if (nargout > 1)
63 if trans_A 112 if (issquare (A))
64 A = matrix_type (A, 'upper'); 113 R = rcond (A);
65 else 114 else
66 A = matrix_type (A, 'lower'); 115 R = 0;
67 endif 116 endif
68 endif 117 endif
69 if isfield(options, 'UT') && options.UT 118 endfunction
70 if trans_A
71 A = matrix_type (A, 'lower');
72 else
73 A = matrix_type (A, 'upper');
74 endif
75 endif
76 endif
77 119
78 X = A \ B;
79 120
80 if nargout > 1 121 %!test
81 if issquare(A) 122 %! n = 4;
82 R = 1 ./ cond(A); 123 %! A = triu (rand (n));
83 else 124 %! x = rand (n, 1);
84 R = 0; 125 %! b = A' * x;
85 endif 126 %! opts.UT = true;
86 endif 127 %! opts.TRANSA = true;
128 %! assert (linsolve (A, b, opts), A' \ b);
87 129
88 %!shared n, A, B, x, opts 130 %!error linsolve ()
89 %! n = 4; A = triu(rand(n)); x = rand(n, 1); B = A' * x; 131 %!error linsolve (1)
90 %! opts.UT = true; opts.TRANSA = true; 132 %!error linsolve (1,2,3)
91 %!assert (linsolve(A, B, opts), A' \ B); 133 %!error <A and B must be numeric> linsolve ({1},2)
92 134 %!error <A and B must be numeric> linsolve (1,{2})
135 %!error <OPTS must be a structure> linsolve (1,2,3)