Mercurial > octave

comparison libinterp/corefcn/lu.cc @ 28154:e7fe6703a81f stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
doc: Update documentation for lu(). * lu.cc (Flu): Rename input "thres" to "tresh" for clarity. Use uppercase letters for L,U,P.
author Rik <rik@octave.org>
date Thu, 12 Mar 2020 16:56:18 -0700
parents bd51beb6205e
children 09c071328135 0a5b15007766
comparison
equal deleted inserted replaced
28152:4609d001daee 28154:e7fe6703a81f
63 doc: /* -*- texinfo -*- 63 doc: /* -*- texinfo -*-
64 @deftypefn {} {[@var{L}, @var{U}] =} lu (@var{A}) 64 @deftypefn {} {[@var{L}, @var{U}] =} lu (@var{A})
65 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} lu (@var{A}) 65 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} lu (@var{A})
66 @deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} lu (@var{S}) 66 @deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} lu (@var{S})
67 @deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}, @var{R}] =} lu (@var{S}) 67 @deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}, @var{R}] =} lu (@var{S})
68 @deftypefnx {} {[@dots{}] =} lu (@var{S}, @var{thres}) 68 @deftypefnx {} {[@dots{}] =} lu (@var{S}, @var{thresh})
69 @deftypefnx {} {@var{y} =} lu (@dots{}) 69 @deftypefnx {} {@var{y} =} lu (@dots{})
70 @deftypefnx {} {[@dots{}] =} lu (@dots{}, "vector") 70 @deftypefnx {} {[@dots{}] =} lu (@dots{}, "vector")
71 @cindex LU decomposition 71 @cindex LU decomposition
72 Compute the LU@tie{}decomposition of @var{A}. 72 Compute the LU@tie{}decomposition of @var{A}.
73 73
74 If @var{A} is full then subroutines from @sc{lapack} are used, and if 74 If @var{A} is full then subroutines from @sc{lapack} are used, and if
75 @var{A} is sparse then @sc{umfpack} is used. 75 @var{A} is sparse then @sc{umfpack} is used.
76 76
77 The result is returned in a permuted form, according to the optional return 77 The result is returned in a permuted form, according to the optional return
78 value @var{P}. For example, given the matrix @code{a = [1, 2; 3, 4]}, 78 value @var{P}. For example, given the matrix @code{@var{A} = [1, 2; 3, 4]},
79 79
80 @example 80 @example
81 [l, u, p] = lu (@var{a}) 81 [@var{L}, @var{U}, @var{P}] = lu (@var{A})
82 @end example 82 @end example
83 83
84 @noindent 84 @noindent
85 returns 85 returns
86 86
87 @example 87 @example
88 @group 88 @group
89 l = 89 L =
90 90
91 1.00000 0.00000 91 1.00000 0.00000
92 0.33333 1.00000 92 0.33333 1.00000
93 93
94 u = 94 U =
95 95
96 3.00000 4.00000 96 3.00000 4.00000
97 0.00000 0.66667 97 0.00000 0.66667
98 98
99 p = 99 P =
100 100
101 0 1 101 0 1
102 1 0 102 1 0
103 @end group 103 @end group
104 @end example 104 @end example
115 Called with a fifth output argument and a sparse input matrix, @code{lu} 115 Called with a fifth output argument and a sparse input matrix, @code{lu}
116 attempts to use a scaling factor @var{R} on the input matrix such that 116 attempts to use a scaling factor @var{R} on the input matrix such that
117 @code{@var{P} * (@var{R} \ @var{A}) * @var{Q} = @var{L} * @var{U}}. 117 @code{@var{P} * (@var{R} \ @var{A}) * @var{Q} = @var{L} * @var{U}}.
118 This typically leads to a sparser and more stable factorization. 118 This typically leads to a sparser and more stable factorization.
119 119
120 An additional input argument @var{thres}, that defines the pivoting 120 An additional input argument @var{thresh} that defines the pivoting
121 threshold can be given. @var{thres} can be a scalar, in which case 121 threshold can be given. @var{thresh} can be a scalar, in which case
122 it defines the @sc{umfpack} pivoting tolerance for both symmetric and 122 it defines the @sc{umfpack} pivoting tolerance for both symmetric and
123 unsymmetric cases. If @var{thres} is a 2-element vector, then the first 123 unsymmetric cases. If @var{thresh} is a 2-element vector, then the first
124 element defines the pivoting tolerance for the unsymmetric @sc{umfpack} 124 element defines the pivoting tolerance for the unsymmetric @sc{umfpack}
125 pivoting strategy and the second for the symmetric strategy. By default, 125 pivoting strategy and the second for the symmetric strategy. By default,
126 the values defined by @code{spparms} are used ([0.1, 0.001]). 126 the values defined by @code{spparms} are used ([0.1, 0.001]).
127 127
128 Given the string argument @qcode{"vector"}, @code{lu} returns the values 128 Given the string argument @qcode{"vector"}, @code{lu} returns the values
145 145
146 if (nargin < 1 || (issparse && nargin > 3) || (! issparse && nargin > 2)) 146 if (nargin < 1 || (issparse && nargin > 3) || (! issparse && nargin > 2))
147 print_usage (); 147 print_usage ();
148 148
149 bool vecout = false; 149 bool vecout = false;
150 Matrix thres; 150 Matrix thresh;
151 int n = 1; 151 int n = 1;
152 152
153 while (n < nargin) 153 while (n < nargin)
154 { 154 {
155 if (args(n).is_string ()) 155 if (args(n).is_string ())
162 error ("lu: unrecognized string argument"); 162 error ("lu: unrecognized string argument");
163 } 163 }
164 else 164 else
165 { 165 {
166 if (! issparse) 166 if (! issparse)
167 error ("lu: can not define pivoting threshold THRES for full matrices"); 167 error ("lu: can not define pivoting threshold THRESH for full matrices");
168 168
169 Matrix tmp = args(n++).matrix_value (); 169 Matrix tmp = args(n++).matrix_value ();
170 if (tmp.numel () == 1) 170 if (tmp.numel () == 1)
171 { 171 {
172 thres.resize (1,2); 172 thresh.resize (1,2);
173 thres(0) = tmp(0); 173 thresh(0) = tmp(0);
174 thres(1) = tmp(0); 174 thresh(1) = tmp(0);
175 } 175 }
176 else if (tmp.numel () == 2) 176 else if (tmp.numel () == 2)
177 thres = tmp; 177 thresh = tmp;
178 else 178 else
179 error ("lu: THRES must be a 1 or 2-element vector"); 179 error ("lu: THRESH must be a 1- or 2-element vector");
180 } 180 }
181 } 181 }
182 182
183 octave_value_list retval; 183 octave_value_list retval;
184 bool scale = (nargout == 5); 184 bool scale = (nargout == 5);
203 "lu: function may fail when called with less than 4 output arguments and a sparse input"); 203 "lu: function may fail when called with less than 4 output arguments and a sparse input");
204 204
205 ColumnVector Qinit (nc); 205 ColumnVector Qinit (nc);
206 for (octave_idx_type i = 0; i < nc; i++) 206 for (octave_idx_type i = 0; i < nc; i++)
207 Qinit(i) = i; 207 Qinit(i) = i;
208 octave::math::sparse_lu<SparseMatrix> fact (m, Qinit, thres, 208 octave::math::sparse_lu<SparseMatrix> fact (m, Qinit, thresh,
209 false, true); 209 false, true);
210 210
211 if (nargout < 2) 211 if (nargout < 2)
212 retval(0) = fact.Y (); 212 retval(0) = fact.Y ();
213 else 213 else
237 } 237 }
238 } 238 }
239 else 239 else
240 { 240 {
241 retval.resize (scale ? 5 : 4); 241 retval.resize (scale ? 5 : 4);
242 octave::math::sparse_lu<SparseMatrix> fact (m, thres, scale); 242 octave::math::sparse_lu<SparseMatrix> fact (m, thresh, scale);
243 243
244 retval(0) = octave_value (fact.L (), 244 retval(0) = octave_value (fact.L (),
245 MatrixType (MatrixType::Lower)); 245 MatrixType (MatrixType::Lower));
246 retval(1) = octave_value (fact.U (), 246 retval(1) = octave_value (fact.U (),
247 MatrixType (MatrixType::Upper)); 247 MatrixType (MatrixType::Upper));
272 272
273 ColumnVector Qinit (nc); 273 ColumnVector Qinit (nc);
274 for (octave_idx_type i = 0; i < nc; i++) 274 for (octave_idx_type i = 0; i < nc; i++)
275 Qinit(i) = i; 275 Qinit(i) = i;
276 octave::math::sparse_lu<SparseComplexMatrix> fact (m, Qinit, 276 octave::math::sparse_lu<SparseComplexMatrix> fact (m, Qinit,
277 thres, false, 277 thresh, false,
278 true); 278 true);
279 279
280 if (nargout < 2) 280 if (nargout < 2)
281 retval(0) = fact.Y (); 281 retval(0) = fact.Y ();
282 else 282 else
305 } 305 }
306 } 306 }
307 else 307 else
308 { 308 {
309 retval.resize (scale ? 5 : 4); 309 retval.resize (scale ? 5 : 4);
310 octave::math::sparse_lu<SparseComplexMatrix> fact (m, thres, 310 octave::math::sparse_lu<SparseComplexMatrix> fact (m, thresh,
311 scale); 311 scale);
312 312
313 retval(0) = octave_value (fact.L (), 313 retval(0) = octave_value (fact.L (),
314 MatrixType (MatrixType::Lower)); 314 MatrixType (MatrixType::Lower));
315 retval(1) = octave_value (fact.U (), 315 retval(1) = octave_value (fact.U (),