Mercurial > octave

comparison src/corefcn/lu.cc @ 15039:e753177cde93

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory * __contourc__.cc, __dispatch__.cc, __lin_interpn__.cc, __pchip_deriv__.cc, __qp__.cc, balance.cc, besselj.cc, betainc.cc, bsxfun.cc, cellfun.cc, colloc.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, det.cc, dlmread.cc, dot.cc, eig.cc, fft.cc, fft2.cc, fftn.cc, filter.cc, find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc, hess.cc, hex2num.cc, inv.cc, kron.cc, lookup.cc, lsode.cc, lu.cc, luinc.cc, matrix_type.cc, max.cc, md5sum.cc, mgorth.cc, nproc.cc, pinv.cc, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, spparms.cc, sqrtm.cc, str2double.cc, strfind.cc, sub2ind.cc, svd.cc, syl.cc, time.cc, tril.cc, typecast.cc: Move functions from DLD-FUNCTIONS/ to corefcn/ directory. Include "defun.h", not "defun-dld.h". Change docstring to refer to these as "Built-in Functions". * build-aux/mk-opts.pl: Generate options code with '#include "defun.h"'. Change option docstrings to refer to these as "Built-in Functions". * corefcn/module.mk: List of functions to build in corefcn/ dir. * DLD-FUNCTIONS/config-module.awk: Update to new build system. * DLD-FUNCTIONS/module-files: Remove functions which are now in corefcn/ directory. * src/Makefile.am: Update to build "convenience library" in corefcn/. Octave program now links against all other libraries + corefcn libary. * src/find-defun-files.sh: Strip $srcdir from filename. * src/link-deps.mk: Add REGEX and FFTW link dependencies for liboctinterp. * type.m, which.m: Change failing tests to use 'amd', still a dynamic function, rather than 'dot', which isn't.
author Rik <rik@octave.org>
date Fri, 27 Jul 2012 15:35:00 -0700
parents src/DLD-FUNCTIONS/lu.cc@5ae9f0f77635
children
comparison
equal deleted inserted replaced
15038:ab18578c2ade 15039:e753177cde93
1 /*
2
3 Copyright (C) 1996-2012 John W. Eaton
4
5 This file is part of Octave.
6
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20
21 */
22
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26
27 #include "CmplxLU.h"
28 #include "dbleLU.h"
29 #include "fCmplxLU.h"
30 #include "floatLU.h"
31 #include "SparseCmplxLU.h"
32 #include "SparsedbleLU.h"
33
34 #include "defun.h"
35 #include "error.h"
36 #include "gripes.h"
37 #include "oct-obj.h"
38 #include "utils.h"
39 #include "ov-re-sparse.h"
40 #include "ov-cx-sparse.h"
41
42 template <class MT>
43 static octave_value
44 get_lu_l (const base_lu<MT>& fact)
45 {
46 MT L = fact.L ();
47 if (L.is_square ())
48 return octave_value (L, MatrixType (MatrixType::Lower));
49 else
50 return L;
51 }
52
53 template <class MT>
54 static octave_value
55 get_lu_u (const base_lu<MT>& fact)
56 {
57 MT U = fact.U ();
58 if (U.is_square () && fact.regular ())
59 return octave_value (U, MatrixType (MatrixType::Upper));
60 else
61 return U;
62 }
63
64 DEFUN (lu, args, nargout,
65 "-*- texinfo -*-\n\
66 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} lu (@var{A})\n\
67 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} lu (@var{A})\n\
68 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} lu (@var{S})\n\
69 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}, @var{R}] =} lu (@var{S})\n\
70 @deftypefnx {Built-in Function} {[@dots{}] =} lu (@var{S}, @var{thres})\n\
71 @deftypefnx {Built-in Function} {@var{y} =} lu (@dots{})\n\
72 @deftypefnx {Built-in Function} {[@dots{}] =} lu (@dots{}, \"vector\")\n\
73 @cindex LU decomposition\n\
74 Compute the LU@tie{}decomposition of @var{A}. If @var{A} is full\n\
75 subroutines from\n\
76 @sc{lapack} are used and if @var{A} is sparse then @sc{umfpack} is used. The\n\
77 result is returned in a permuted form, according to the optional return\n\
78 value @var{P}. For example, given the matrix @code{a = [1, 2; 3, 4]},\n\
79 \n\
80 @example\n\
81 [l, u, p] = lu (@var{a})\n\
82 @end example\n\
83 \n\
84 @noindent\n\
85 returns\n\
86 \n\
87 @example\n\
88 @group\n\
89 l =\n\
90 \n\
91 1.00000 0.00000\n\
92 0.33333 1.00000\n\
93 \n\
94 u =\n\
95 \n\
96 3.00000 4.00000\n\
97 0.00000 0.66667\n\
98 \n\
99 p =\n\
100 \n\
101 0 1\n\
102 1 0\n\
103 @end group\n\
104 @end example\n\
105 \n\
106 The matrix is not required to be square.\n\
107 \n\
108 When called with two or three output arguments and a spare input matrix,\n\
109 @code{lu} does not attempt to perform sparsity preserving column\n\
110 permutations. Called with a fourth output argument, the sparsity\n\
111 preserving column transformation @var{Q} is returned, such that\n\
112 @code{@var{P} * @var{A} * @var{Q} = @var{L} * @var{U}}.\n\
113 \n\
114 Called with a fifth output argument and a sparse input matrix,\n\
115 @code{lu} attempts to use a scaling factor @var{R} on the input matrix\n\
116 such that\n\
117 @code{@var{P} * (@var{R} \\ @var{A}) * @var{Q} = @var{L} * @var{U}}.\n\
118 This typically leads to a sparser and more stable factorization.\n\
119 \n\
120 An additional input argument @var{thres}, that defines the pivoting\n\
121 threshold can be given. @var{thres} can be a scalar, in which case\n\
122 it defines the @sc{umfpack} pivoting tolerance for both symmetric and\n\
123 unsymmetric cases. If @var{thres} is a 2-element vector, then the first\n\
124 element defines the pivoting tolerance for the unsymmetric @sc{umfpack}\n\
125 pivoting strategy and the second for the symmetric strategy. By default,\n\
126 the values defined by @code{spparms} are used ([0.1, 0.001]).\n\
127 \n\
128 Given the string argument \"vector\", @code{lu} returns the values of @var{P}\n\
129 and @var{Q} as vector values, such that for full matrix, @code{@var{A}\n\
130 (@var{P},:) = @var{L} * @var{U}}, and @code{@var{R}(@var{P},:) * @var{A}\n\
131 (:, @var{Q}) = @var{L} * @var{U}}.\n\
132 \n\
133 With two output arguments, returns the permuted forms of the upper and\n\
134 lower triangular matrices, such that @code{@var{A} = @var{L} * @var{U}}.\n\
135 With one output argument @var{y}, then the matrix returned by the @sc{lapack}\n\
136 routines is returned. If the input matrix is sparse then the matrix @var{L}\n\
137 is embedded into @var{U} to give a return value similar to the full case.\n\
138 For both full and sparse matrices, @code{lu} loses the permutation\n\
139 information.\n\
140 @end deftypefn")
141 {
142 octave_value_list retval;
143 int nargin = args.length ();
144 bool issparse = (nargin > 0 && args(0).is_sparse_type ());
145 bool scale = (nargout == 5);
146
147 if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5))
148 || (!issparse && (nargin > 2 || nargout > 3)))
149 {
150 print_usage ();
151 return retval;
152 }
153
154 bool vecout = false;
155 Matrix thres;
156
157 int n = 1;
158 while (n < nargin && ! error_state)
159 {
160 if (args (n).is_string ())
161 {
162 std::string tmp = args(n++).string_value ();
163
164 if (! error_state )
165 {
166 if (tmp.compare ("vector") == 0)
167 vecout = true;
168 else
169 error ("lu: unrecognized string argument");
170 }
171 }
172 else
173 {
174 Matrix tmp = args(n++).matrix_value ();
175
176 if (! error_state )
177 {
178 if (!issparse)
179 error ("lu: can not define pivoting threshold THRES for full matrices");
180 else if (tmp.nelem () == 1)
181 {
182 thres.resize (1,2);
183 thres(0) = tmp(0);
184 thres(1) = tmp(0);
185 }
186 else if (tmp.nelem () == 2)
187 thres = tmp;
188 else
189 error ("lu: expecting 2-element vector for THRES");
190 }
191 }
192 }
193
194 octave_value arg = args(0);
195
196 octave_idx_type nr = arg.rows ();
197 octave_idx_type nc = arg.columns ();
198
199 int arg_is_empty = empty_arg ("lu", nr, nc);
200
201 if (issparse)
202 {
203 if (arg_is_empty < 0)
204 return retval;
205 else if (arg_is_empty > 0)
206 return octave_value_list (5, SparseMatrix ());
207
208 ColumnVector Qinit;
209 if (nargout < 4)
210 {
211 Qinit.resize (nc);
212 for (octave_idx_type i = 0; i < nc; i++)
213 Qinit (i) = i;
214 }
215
216 if (arg.is_real_type ())
217 {
218 SparseMatrix m = arg.sparse_matrix_value ();
219
220 switch (nargout)
221 {
222 case 0:
223 case 1:
224 case 2:
225 {
226 SparseLU fact (m, Qinit, thres, false, true);
227
228 if (nargout < 2)
229 retval(0) = fact.Y ();
230 else
231 {
232 PermMatrix P = fact.Pr_mat ();
233 SparseMatrix L = P.transpose () * fact.L ();
234 retval(1) = octave_value (fact.U (),
235 MatrixType (MatrixType::Upper));
236
237 retval(0) = octave_value (L,
238 MatrixType (MatrixType::Permuted_Lower,
239 nr, fact.row_perm ()));
240 }
241 }
242 break;
243
244 case 3:
245 {
246 SparseLU fact (m, Qinit, thres, false, true);
247
248 if (vecout)
249 retval(2) = fact.Pr_vec ();
250 else
251 retval(2) = fact.Pr_mat ();
252
253 retval(1) = octave_value (fact.U (),
254 MatrixType (MatrixType::Upper));
255 retval(0) = octave_value (fact.L (),
256 MatrixType (MatrixType::Lower));
257 }
258 break;
259
260 case 4:
261 default:
262 {
263 SparseLU fact (m, thres, scale);
264
265 if (scale)
266 retval(4) = fact.R ();
267
268 if (vecout)
269 {
270 retval(3) = fact.Pc_vec ();
271 retval(2) = fact.Pr_vec ();
272 }
273 else
274 {
275 retval(3) = fact.Pc_mat ();
276 retval(2) = fact.Pr_mat ();
277 }
278 retval(1) = octave_value (fact.U (),
279 MatrixType (MatrixType::Upper));
280 retval(0) = octave_value (fact.L (),
281 MatrixType (MatrixType::Lower));
282 }
283 break;
284 }
285 }
286 else if (arg.is_complex_type ())
287 {
288 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
289
290 switch (nargout)
291 {
292 case 0:
293 case 1:
294 case 2:
295 {
296 SparseComplexLU fact (m, Qinit, thres, false, true);
297
298 if (nargout < 2)
299 retval(0) = fact.Y ();
300 else
301 {
302 PermMatrix P = fact.Pr_mat ();
303 SparseComplexMatrix L = P.transpose () * fact.L ();
304 retval(1) = octave_value (fact.U (),
305 MatrixType (MatrixType::Upper));
306
307 retval(0) = octave_value (L,
308 MatrixType (MatrixType::Permuted_Lower,
309 nr, fact.row_perm ()));
310 }
311 }
312 break;
313
314 case 3:
315 {
316 SparseComplexLU fact (m, Qinit, thres, false, true);
317
318 if (vecout)
319 retval(2) = fact.Pr_vec ();
320 else
321 retval(2) = fact.Pr_mat ();
322
323 retval(1) = octave_value (fact.U (),
324 MatrixType (MatrixType::Upper));
325 retval(0) = octave_value (fact.L (),
326 MatrixType (MatrixType::Lower));
327 }
328 break;
329
330 case 4:
331 default:
332 {
333 SparseComplexLU fact (m, thres, scale);
334
335 if (scale)
336 retval(4) = fact.R ();
337
338 if (vecout)
339 {
340 retval(3) = fact.Pc_vec ();
341 retval(2) = fact.Pr_vec ();
342 }
343 else
344 {
345 retval(3) = fact.Pc_mat ();
346 retval(2) = fact.Pr_mat ();
347 }
348 retval(1) = octave_value (fact.U (),
349 MatrixType (MatrixType::Upper));
350 retval(0) = octave_value (fact.L (),
351 MatrixType (MatrixType::Lower));
352 }
353 break;
354 }
355 }
356 else
357 gripe_wrong_type_arg ("lu", arg);
358 }
359 else
360 {
361 if (arg_is_empty < 0)
362 return retval;
363 else if (arg_is_empty > 0)
364 return octave_value_list (3, Matrix ());
365
366 if (arg.is_real_type ())
367 {
368 if (arg.is_single_type ())
369 {
370 FloatMatrix m = arg.float_matrix_value ();
371
372 if (! error_state)
373 {
374 FloatLU fact (m);
375
376 switch (nargout)
377 {
378 case 0:
379 case 1:
380 retval(0) = fact.Y ();
381 break;
382
383 case 2:
384 {
385 PermMatrix P = fact.P ();
386 FloatMatrix L = P.transpose () * fact.L ();
387 retval(1) = get_lu_u (fact);
388 retval(0) = L;
389 }
390 break;
391
392 case 3:
393 default:
394 {
395 if (vecout)
396 retval(2) = fact.P_vec ();
397 else
398 retval(2) = fact.P ();
399 retval(1) = get_lu_u (fact);
400 retval(0) = get_lu_l (fact);
401 }
402 break;
403 }
404 }
405 }
406 else
407 {
408 Matrix m = arg.matrix_value ();
409
410 if (! error_state)
411 {
412 LU fact (m);
413
414 switch (nargout)
415 {
416 case 0:
417 case 1:
418 retval(0) = fact.Y ();
419 break;
420
421 case 2:
422 {
423 PermMatrix P = fact.P ();
424 Matrix L = P.transpose () * fact.L ();
425 retval(1) = get_lu_u (fact);
426 retval(0) = L;
427 }
428 break;
429
430 case 3:
431 default:
432 {
433 if (vecout)
434 retval(2) = fact.P_vec ();
435 else
436 retval(2) = fact.P ();
437 retval(1) = get_lu_u (fact);
438 retval(0) = get_lu_l (fact);
439 }
440 break;
441 }
442 }
443 }
444 }
445 else if (arg.is_complex_type ())
446 {
447 if (arg.is_single_type ())
448 {
449 FloatComplexMatrix m = arg.float_complex_matrix_value ();
450
451 if (! error_state)
452 {
453 FloatComplexLU fact (m);
454
455 switch (nargout)
456 {
457 case 0:
458 case 1:
459 retval(0) = fact.Y ();
460 break;
461
462 case 2:
463 {
464 PermMatrix P = fact.P ();
465 FloatComplexMatrix L = P.transpose () * fact.L ();
466 retval(1) = get_lu_u (fact);
467 retval(0) = L;
468 }
469 break;
470
471 case 3:
472 default:
473 {
474 if (vecout)
475 retval(2) = fact.P_vec ();
476 else
477 retval(2) = fact.P ();
478 retval(1) = get_lu_u (fact);
479 retval(0) = get_lu_l (fact);
480 }
481 break;
482 }
483 }
484 }
485 else
486 {
487 ComplexMatrix m = arg.complex_matrix_value ();
488
489 if (! error_state)
490 {
491 ComplexLU fact (m);
492
493 switch (nargout)
494 {
495 case 0:
496 case 1:
497 retval(0) = fact.Y ();
498 break;
499
500 case 2:
501 {
502 PermMatrix P = fact.P ();
503 ComplexMatrix L = P.transpose () * fact.L ();
504 retval(1) = get_lu_u (fact);
505 retval(0) = L;
506 }
507 break;
508
509 case 3:
510 default:
511 {
512 if (vecout)
513 retval(2) = fact.P_vec ();
514 else
515 retval(2) = fact.P ();
516 retval(1) = get_lu_u (fact);
517 retval(0) = get_lu_l (fact);
518 }
519 break;
520 }
521 }
522 }
523 }
524 else
525 gripe_wrong_type_arg ("lu", arg);
526 }
527
528 return retval;
529 }
530
531 /*
532 %!assert(lu ([1, 2; 3, 4]), [3, 4; 1/3, 2/3], eps);
533
534 %!test
535 %! [l, u] = lu ([1, 2; 3, 4]);
536 %! assert (l, [1/3, 1; 1, 0], sqrt (eps));
537 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
538
539 %!test
540 %! [l, u, p] = lu ([1, 2; 3, 4]);
541 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
542 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
543 %! assert (p(:,:), [0, 1; 1, 0], sqrt (eps));
544
545 %!test
546 %! [l, u, p] = lu ([1, 2; 3, 4], "vector");
547 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
548 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
549 %! assert (p, [2;1], sqrt (eps));
550
551 %!test
552 %! [l, u, p] = lu ([1, 2; 3, 4; 5, 6]);
553 %! assert (l, [1, 0; 1/5, 1; 3/5, 1/2], sqrt (eps));
554 %! assert (u, [5, 6; 0, 4/5], sqrt (eps));
555 %! assert (p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
556
557 %!assert (lu (single ([1, 2; 3, 4])), single ([3, 4; 1/3, 2/3]), eps ("single"))
558
559 %!test
560 %! [l, u] = lu (single ([1, 2; 3, 4]));
561 %! assert (l, single ([1/3, 1; 1, 0]), sqrt (eps ("single")));
562 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
563
564 %!test
565 %! [l, u, p] = lu (single ([1, 2; 3, 4]));
566 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
567 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
568 %! assert (p(:,:), single ([0, 1; 1, 0]), sqrt (eps ("single")));
569
570 %!test
571 %! [l, u, p] = lu (single ([1, 2; 3, 4]), "vector");
572 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
573 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
574 %! assert (p, single ([2;1]), sqrt (eps ("single")));
575
576 %!test
577 %! [l u p] = lu (single ([1, 2; 3, 4; 5, 6]));
578 %! assert (l, single ([1, 0; 1/5, 1; 3/5, 1/2]), sqrt (eps ("single")));
579 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single")));
580 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single")));
581
582 %!error lu ()
583 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2)
584 */
585
586 static
587 bool check_lu_dims (const octave_value& l, const octave_value& u,
588 const octave_value& p)
589 {
590 octave_idx_type m = l.rows (), k = u.rows (), n = u.columns ();
591 return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
592 && k == std::min (m, n) &&
593 (p.is_undefined () || p.rows () == m));
594 }
595
596 DEFUN (luupdate, args, ,
597 "-*- texinfo -*-\n\
598 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} luupdate (@var{L}, @var{U}, @var{x}, @var{y})\n\
599 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
600 Given an LU@tie{}factorization of a real or complex matrix\n\
601 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and\n\
602 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization\n\
603 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are\n\
604 column vectors (rank-1 update) or matrices with equal number of columns\n\
605 (rank-k update).\n\
606 Optionally, row-pivoted updating can be used by supplying\n\
607 a row permutation (pivoting) matrix @var{P};\n\
608 in that case, an updated permutation matrix is returned.\n\
609 Note that if @var{L}, @var{U}, @var{P} is a pivoted LU@tie{}factorization\n\
610 as obtained by @code{lu}:\n\
611 \n\
612 @example\n\
613 [@var{L}, @var{U}, @var{P}] = lu (@var{A});\n\
614 @end example\n\
615 \n\
616 @noindent\n\
617 then a factorization of @xcode{@var{A}+@var{x}*@var{y}.'} can be obtained\n\
618 either as\n\
619 \n\
620 @example\n\
621 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})\n\
622 @end example\n\
623 \n\
624 @noindent\n\
625 or\n\
626 \n\
627 @example\n\
628 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
629 @end example\n\
630 \n\
631 The first form uses the unpivoted algorithm, which is faster, but less\n\
632 stable. The second form uses a slower pivoted algorithm, which is more\n\
633 stable.\n\
634 \n\
635 The matrix case is done as a sequence of rank-1 updates;\n\
636 thus, for large enough k, it will be both faster and more accurate to\n\
637 recompute the factorization from scratch.\n\
638 @seealso{lu, qrupdate, cholupdate}\n\
639 @end deftypefn")
640 {
641 octave_idx_type nargin = args.length ();
642 octave_value_list retval;
643
644 bool pivoted = nargin == 5;
645
646 if (nargin != 4 && nargin != 5)
647 {
648 print_usage ();
649 return retval;
650 }
651
652 octave_value argl = args(0);
653 octave_value argu = args(1);
654 octave_value argp = pivoted ? args(2) : octave_value ();
655 octave_value argx = args(2 + pivoted);
656 octave_value argy = args(3 + pivoted);
657
658 if (argl.is_numeric_type () && argu.is_numeric_type ()
659 && argx.is_numeric_type () && argy.is_numeric_type ()
660 && (! pivoted || argp.is_perm_matrix ()))
661 {
662 if (check_lu_dims (argl, argu, argp))
663 {
664 PermMatrix P = (pivoted
665 ? argp.perm_matrix_value ()
666 : PermMatrix::eye (argl.rows ()));
667
668 if (argl.is_real_type ()
669 && argu.is_real_type ()
670 && argx.is_real_type ()
671 && argy.is_real_type ())
672 {
673 // all real case
674 if (argl.is_single_type ()
675 || argu.is_single_type ()
676 || argx.is_single_type ()
677 || argy.is_single_type ())
678 {
679 FloatMatrix L = argl.float_matrix_value ();
680 FloatMatrix U = argu.float_matrix_value ();
681 FloatMatrix x = argx.float_matrix_value ();
682 FloatMatrix y = argy.float_matrix_value ();
683
684 FloatLU fact (L, U, P);
685 if (pivoted)
686 fact.update_piv (x, y);
687 else
688 fact.update (x, y);
689
690 if (pivoted)
691 retval(2) = fact.P ();
692 retval(1) = get_lu_u (fact);
693 retval(0) = get_lu_l (fact);
694 }
695 else
696 {
697 Matrix L = argl.matrix_value ();
698 Matrix U = argu.matrix_value ();
699 Matrix x = argx.matrix_value ();
700 Matrix y = argy.matrix_value ();
701
702 LU fact (L, U, P);
703 if (pivoted)
704 fact.update_piv (x, y);
705 else
706 fact.update (x, y);
707
708 if (pivoted)
709 retval(2) = fact.P ();
710 retval(1) = get_lu_u (fact);
711 retval(0) = get_lu_l (fact);
712 }
713 }
714 else
715 {
716 // complex case
717 if (argl.is_single_type ()
718 || argu.is_single_type ()
719 || argx.is_single_type ()
720 || argy.is_single_type ())
721 {
722 FloatComplexMatrix L = argl.float_complex_matrix_value ();
723 FloatComplexMatrix U = argu.float_complex_matrix_value ();
724 FloatComplexMatrix x = argx.float_complex_matrix_value ();
725 FloatComplexMatrix y = argy.float_complex_matrix_value ();
726
727 FloatComplexLU fact (L, U, P);
728 if (pivoted)
729 fact.update_piv (x, y);
730 else
731 fact.update (x, y);
732
733 if (pivoted)
734 retval(2) = fact.P ();
735 retval(1) = get_lu_u (fact);
736 retval(0) = get_lu_l (fact);
737 }
738 else
739 {
740 ComplexMatrix L = argl.complex_matrix_value ();
741 ComplexMatrix U = argu.complex_matrix_value ();
742 ComplexMatrix x = argx.complex_matrix_value ();
743 ComplexMatrix y = argy.complex_matrix_value ();
744
745 ComplexLU fact (L, U, P);
746 if (pivoted)
747 fact.update_piv (x, y);
748 else
749 fact.update (x, y);
750
751 if (pivoted)
752 retval(2) = fact.P ();
753 retval(1) = get_lu_u (fact);
754 retval(0) = get_lu_l (fact);
755 }
756 }
757 }
758 else
759 error ("luupdate: dimension mismatch");
760 }
761 else
762 error ("luupdate: L, U, X, and Y must be numeric");
763
764 return retval;
765 }
766
767 /*
768 %!shared A, u, v, Ac, uc, vc
769 %! A = [0.091364 0.613038 0.999083;
770 %! 0.594638 0.425302 0.603537;
771 %! 0.383594 0.291238 0.085574;
772 %! 0.265712 0.268003 0.238409;
773 %! 0.669966 0.743851 0.445057 ];
774 %!
775 %! u = [0.85082;
776 %! 0.76426;
777 %! 0.42883;
778 %! 0.53010;
779 %! 0.80683 ];
780 %!
781 %! v = [0.98810;
782 %! 0.24295;
783 %! 0.43167 ];
784 %!
785 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
786 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
787 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
788 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
789 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
790 %!
791 %! uc = [0.20351 + 0.05401i;
792 %! 0.13141 + 0.43708i;
793 %! 0.29808 + 0.08789i;
794 %! 0.69821 + 0.38844i;
795 %! 0.74871 + 0.25821i ];
796 %!
797 %! vc = [0.85839 + 0.29468i;
798 %! 0.20820 + 0.93090i;
799 %! 0.86184 + 0.34689i ];
800 %!
801
802 %!testif HAVE_QRUPDATE_LUU
803 %! [L,U,P] = lu (A);
804 %! [L,U] = luupdate (L,U,P*u,v);
805 %! assert (norm (vec (tril (L)-L), Inf) == 0);
806 %! assert (norm (vec (triu (U)-U), Inf) == 0);
807 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
808 %!
809 %!testif HAVE_QRUPDATE_LUU
810 %! [L,U,P] = lu (Ac);
811 %! [L,U] = luupdate (L,U,P*uc,vc);
812 %! assert (norm (vec (tril (L)-L), Inf) == 0);
813 %! assert (norm (vec (triu (U)-U), Inf) == 0);
814 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
815
816 %!testif HAVE_QRUPDATE_LUU
817 %! [L,U,P] = lu (single (A));
818 %! [L,U] = luupdate (L,U,P*single (u), single (v));
819 %! assert (norm (vec (tril (L)-L), Inf) == 0);
820 %! assert (norm (vec (triu (U)-U), Inf) == 0);
821 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
822 %!
823 %!testif HAVE_QRUPDATE_LUU
824 %! [L,U,P] = lu (single (Ac));
825 %! [L,U] = luupdate (L,U,P*single (uc),single (vc));
826 %! assert (norm (vec (tril (L)-L), Inf) == 0);
827 %! assert (norm (vec (triu (U)-U), Inf) == 0);
828 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
829
830 %!testif HAVE_QRUPDATE_LUU
831 %! [L,U,P] = lu (A);
832 %! [L,U,P] = luupdate (L,U,P,u,v);
833 %! assert (norm (vec (tril (L)-L), Inf) == 0);
834 %! assert (norm (vec (triu (U)-U), Inf) == 0);
835 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
836 %!
837 %!testif HAVE_QRUPDATE_LUU
838 %! [L,U,P] = lu (Ac);
839 %! [L,U,P] = luupdate (L,U,P,uc,vc);
840 %! assert (norm (vec (tril (L)-L), Inf) == 0);
841 %! assert (norm (vec (triu (U)-U), Inf) == 0);
842 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
843
844 %!testif HAVE_QRUPDATE_LUU
845 %! [L,U,P] = lu (single (A));
846 %! [L,U,P] = luupdate (L,U,P,single (u),single (v));
847 %! assert (norm (vec (tril (L)-L), Inf) == 0);
848 %! assert (norm (vec (triu (U)-U), Inf) == 0);
849 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
850 %!
851 %!testif HAVE_QRUPDATE_LUU
852 %! [L,U,P] = lu (single (Ac));
853 %! [L,U,P] = luupdate (L,U,P,single (uc),single (vc));
854 %! assert (norm (vec (tril (L)-L), Inf) == 0);
855 %! assert (norm (vec (triu (U)-U), Inf) == 0);
856 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
857 */