GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
lu.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2013 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 @qcode{\"vector\"}, @code{lu} returns the values\n\
129 of @var{P} and @var{Q} as vector values, such that for full matrix,\n\
130 @code{@var{A} (@var{P},:) = @var{L} * @var{U}}, and @code{@var{R}(@var{P},:)\n\
131 * @var{A} (:, @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 @seealso{luupdate, chol, hess, qr, qz, schur, svd}\n\
141 @end deftypefn")
142 {
143  octave_value_list retval;
144  int nargin = args.length ();
145  bool issparse = (nargin > 0 && args(0).is_sparse_type ());
146  bool scale = (nargout == 5);
147 
148  if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5))
149  || (!issparse && (nargin > 2 || nargout > 3)))
150  {
151  print_usage ();
152  return retval;
153  }
154 
155  bool vecout = false;
156  Matrix thres;
157 
158  int n = 1;
159  while (n < nargin && ! error_state)
160  {
161  if (args (n).is_string ())
162  {
163  std::string tmp = args(n++).string_value ();
164 
165  if (! error_state )
166  {
167  if (tmp.compare ("vector") == 0)
168  vecout = true;
169  else
170  error ("lu: unrecognized string argument");
171  }
172  }
173  else
174  {
175  Matrix tmp = args(n++).matrix_value ();
176 
177  if (! error_state )
178  {
179  if (!issparse)
180  error ("lu: can not define pivoting threshold THRES for full matrices");
181  else if (tmp.nelem () == 1)
182  {
183  thres.resize (1,2);
184  thres(0) = tmp(0);
185  thres(1) = tmp(0);
186  }
187  else if (tmp.nelem () == 2)
188  thres = tmp;
189  else
190  error ("lu: expecting 2-element vector for THRES");
191  }
192  }
193  }
194 
195  octave_value arg = args(0);
196 
197  octave_idx_type nr = arg.rows ();
198  octave_idx_type nc = arg.columns ();
199 
200  int arg_is_empty = empty_arg ("lu", nr, nc);
201 
202  if (issparse)
203  {
204  if (arg_is_empty < 0)
205  return retval;
206  else if (arg_is_empty > 0)
207  return octave_value_list (5, SparseMatrix ());
208 
209  ColumnVector Qinit;
210  if (nargout < 4)
211  {
212  Qinit.resize (nc);
213  for (octave_idx_type i = 0; i < nc; i++)
214  Qinit (i) = i;
215  }
216 
217  if (arg.is_real_type ())
218  {
219  SparseMatrix m = arg.sparse_matrix_value ();
220 
221  switch (nargout)
222  {
223  case 0:
224  case 1:
225  case 2:
226  {
227  SparseLU fact (m, Qinit, thres, false, true);
228 
229  if (nargout < 2)
230  retval(0) = fact.Y ();
231  else
232  {
233  PermMatrix P = fact.Pr_mat ();
234  SparseMatrix L = P.transpose () * fact.L ();
235  retval(1) = octave_value (fact.U (),
236  MatrixType (MatrixType::Upper));
237 
238  retval(0)
239  = octave_value (L, MatrixType (MatrixType::Permuted_Lower,
240  nr, fact.row_perm ()));
241  }
242  }
243  break;
244 
245  case 3:
246  {
247  SparseLU fact (m, Qinit, thres, false, true);
248 
249  if (vecout)
250  retval(2) = fact.Pr_vec ();
251  else
252  retval(2) = fact.Pr_mat ();
253 
254  retval(1) = octave_value (fact.U (),
255  MatrixType (MatrixType::Upper));
256  retval(0) = octave_value (fact.L (),
257  MatrixType (MatrixType::Lower));
258  }
259  break;
260 
261  case 4:
262  default:
263  {
264  SparseLU fact (m, thres, scale);
265 
266  if (scale)
267  retval(4) = fact.R ();
268 
269  if (vecout)
270  {
271  retval(3) = fact.Pc_vec ();
272  retval(2) = fact.Pr_vec ();
273  }
274  else
275  {
276  retval(3) = fact.Pc_mat ();
277  retval(2) = fact.Pr_mat ();
278  }
279  retval(1) = octave_value (fact.U (),
280  MatrixType (MatrixType::Upper));
281  retval(0) = octave_value (fact.L (),
282  MatrixType (MatrixType::Lower));
283  }
284  break;
285  }
286  }
287  else if (arg.is_complex_type ())
288  {
289  SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
290 
291  switch (nargout)
292  {
293  case 0:
294  case 1:
295  case 2:
296  {
297  SparseComplexLU fact (m, Qinit, thres, false, true);
298 
299  if (nargout < 2)
300  retval(0) = fact.Y ();
301  else
302  {
303  PermMatrix P = fact.Pr_mat ();
304  SparseComplexMatrix L = P.transpose () * fact.L ();
305  retval(1) = octave_value (fact.U (),
306  MatrixType (MatrixType::Upper));
307 
308  retval(0)
309  = octave_value (L, MatrixType (MatrixType::Permuted_Lower,
310  nr, fact.row_perm ()));
311  }
312  }
313  break;
314 
315  case 3:
316  {
317  SparseComplexLU fact (m, Qinit, thres, false, true);
318 
319  if (vecout)
320  retval(2) = fact.Pr_vec ();
321  else
322  retval(2) = fact.Pr_mat ();
323 
324  retval(1) = octave_value (fact.U (),
325  MatrixType (MatrixType::Upper));
326  retval(0) = octave_value (fact.L (),
327  MatrixType (MatrixType::Lower));
328  }
329  break;
330 
331  case 4:
332  default:
333  {
334  SparseComplexLU fact (m, thres, scale);
335 
336  if (scale)
337  retval(4) = fact.R ();
338 
339  if (vecout)
340  {
341  retval(3) = fact.Pc_vec ();
342  retval(2) = fact.Pr_vec ();
343  }
344  else
345  {
346  retval(3) = fact.Pc_mat ();
347  retval(2) = fact.Pr_mat ();
348  }
349  retval(1) = octave_value (fact.U (),
350  MatrixType (MatrixType::Upper));
351  retval(0) = octave_value (fact.L (),
352  MatrixType (MatrixType::Lower));
353  }
354  break;
355  }
356  }
357  else
358  gripe_wrong_type_arg ("lu", arg);
359  }
360  else
361  {
362  if (arg_is_empty < 0)
363  return retval;
364  else if (arg_is_empty > 0)
365  return octave_value_list (3, Matrix ());
366 
367  if (arg.is_real_type ())
368  {
369  if (arg.is_single_type ())
370  {
371  FloatMatrix m = arg.float_matrix_value ();
372 
373  if (! error_state)
374  {
375  FloatLU fact (m);
376 
377  switch (nargout)
378  {
379  case 0:
380  case 1:
381  retval(0) = fact.Y ();
382  break;
383 
384  case 2:
385  {
386  PermMatrix P = fact.P ();
387  FloatMatrix L = P.transpose () * fact.L ();
388  retval(1) = get_lu_u (fact);
389  retval(0) = L;
390  }
391  break;
392 
393  case 3:
394  default:
395  {
396  if (vecout)
397  retval(2) = fact.P_vec ();
398  else
399  retval(2) = fact.P ();
400  retval(1) = get_lu_u (fact);
401  retval(0) = get_lu_l (fact);
402  }
403  break;
404  }
405  }
406  }
407  else
408  {
409  Matrix m = arg.matrix_value ();
410 
411  if (! error_state)
412  {
413  LU fact (m);
414 
415  switch (nargout)
416  {
417  case 0:
418  case 1:
419  retval(0) = fact.Y ();
420  break;
421 
422  case 2:
423  {
424  PermMatrix P = fact.P ();
425  Matrix L = P.transpose () * fact.L ();
426  retval(1) = get_lu_u (fact);
427  retval(0) = L;
428  }
429  break;
430 
431  case 3:
432  default:
433  {
434  if (vecout)
435  retval(2) = fact.P_vec ();
436  else
437  retval(2) = fact.P ();
438  retval(1) = get_lu_u (fact);
439  retval(0) = get_lu_l (fact);
440  }
441  break;
442  }
443  }
444  }
445  }
446  else if (arg.is_complex_type ())
447  {
448  if (arg.is_single_type ())
449  {
450  FloatComplexMatrix m = arg.float_complex_matrix_value ();
451 
452  if (! error_state)
453  {
454  FloatComplexLU fact (m);
455 
456  switch (nargout)
457  {
458  case 0:
459  case 1:
460  retval(0) = fact.Y ();
461  break;
462 
463  case 2:
464  {
465  PermMatrix P = fact.P ();
466  FloatComplexMatrix L = P.transpose () * fact.L ();
467  retval(1) = get_lu_u (fact);
468  retval(0) = L;
469  }
470  break;
471 
472  case 3:
473  default:
474  {
475  if (vecout)
476  retval(2) = fact.P_vec ();
477  else
478  retval(2) = fact.P ();
479  retval(1) = get_lu_u (fact);
480  retval(0) = get_lu_l (fact);
481  }
482  break;
483  }
484  }
485  }
486  else
487  {
488  ComplexMatrix m = arg.complex_matrix_value ();
489 
490  if (! error_state)
491  {
492  ComplexLU fact (m);
493 
494  switch (nargout)
495  {
496  case 0:
497  case 1:
498  retval(0) = fact.Y ();
499  break;
500 
501  case 2:
502  {
503  PermMatrix P = fact.P ();
504  ComplexMatrix L = P.transpose () * fact.L ();
505  retval(1) = get_lu_u (fact);
506  retval(0) = L;
507  }
508  break;
509 
510  case 3:
511  default:
512  {
513  if (vecout)
514  retval(2) = fact.P_vec ();
515  else
516  retval(2) = fact.P ();
517  retval(1) = get_lu_u (fact);
518  retval(0) = get_lu_l (fact);
519  }
520  break;
521  }
522  }
523  }
524  }
525  else
526  gripe_wrong_type_arg ("lu", arg);
527  }
528 
529  return retval;
530 }
531 
532 /*
533 %!assert(lu ([1, 2; 3, 4]), [3, 4; 1/3, 2/3], eps);
534 
535 %!test
536 %! [l, u] = lu ([1, 2; 3, 4]);
537 %! assert (l, [1/3, 1; 1, 0], sqrt (eps));
538 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
539 
540 %!test
541 %! [l, u, p] = lu ([1, 2; 3, 4]);
542 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
543 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
544 %! assert (p(:,:), [0, 1; 1, 0], sqrt (eps));
545 
546 %!test
547 %! [l, u, p] = lu ([1, 2; 3, 4], "vector");
548 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
549 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
550 %! assert (p, [2;1], sqrt (eps));
551 
552 %!test
553 %! [l, u, p] = lu ([1, 2; 3, 4; 5, 6]);
554 %! assert (l, [1, 0; 1/5, 1; 3/5, 1/2], sqrt (eps));
555 %! assert (u, [5, 6; 0, 4/5], sqrt (eps));
556 %! assert (p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
557 
558 %!assert (lu (single ([1, 2; 3, 4])), single ([3, 4; 1/3, 2/3]), eps ("single"))
559 
560 %!test
561 %! [l, u] = lu (single ([1, 2; 3, 4]));
562 %! assert (l, single ([1/3, 1; 1, 0]), sqrt (eps ("single")));
563 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
564 
565 %!test
566 %! [l, u, p] = lu (single ([1, 2; 3, 4]));
567 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
568 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
569 %! assert (p(:,:), single ([0, 1; 1, 0]), sqrt (eps ("single")));
570 
571 %!test
572 %! [l, u, p] = lu (single ([1, 2; 3, 4]), "vector");
573 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
574 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
575 %! assert (p, single ([2;1]), sqrt (eps ("single")));
576 
577 %!test
578 %! [l u p] = lu (single ([1, 2; 3, 4; 5, 6]));
579 %! assert (l, single ([1, 0; 1/5, 1; 3/5, 1/2]), sqrt (eps ("single")));
580 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single")));
581 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single")));
582 
583 %!error lu ()
584 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2)
585 */
586 
587 static
588 bool check_lu_dims (const octave_value& l, const octave_value& u,
589  const octave_value& p)
590 {
591  octave_idx_type m = l.rows (), k = u.rows (), n = u.columns ();
592  return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
593  && k == std::min (m, n) &&
594  (p.is_undefined () || p.rows () == m));
595 }
596 
597 DEFUN (luupdate, args, ,
598  "-*- texinfo -*-\n\
599 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} luupdate (@var{L}, @var{U}, @var{x}, @var{y})\n\
600 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
601 Given an LU@tie{}factorization of a real or complex matrix\n\
602 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and\n\
603 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization\n\
604 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are\n\
605 column vectors (rank-1 update) or matrices with equal number of columns\n\
606 (rank-k update).\n\
607 Optionally, row-pivoted updating can be used by supplying\n\
608 a row permutation (pivoting) matrix @var{P};\n\
609 in that case, an updated permutation matrix is returned.\n\
610 Note that if @var{L}, @var{U}, @var{P} is a pivoted LU@tie{}factorization\n\
611 as obtained by @code{lu}:\n\
612 \n\
613 @example\n\
614 [@var{L}, @var{U}, @var{P}] = lu (@var{A});\n\
615 @end example\n\
616 \n\
617 @noindent\n\
618 then a factorization of @tcode{@var{A}+@var{x}*@var{y}.'} can be obtained\n\
619 either as\n\
620 \n\
621 @example\n\
622 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})\n\
623 @end example\n\
624 \n\
625 @noindent\n\
626 or\n\
627 \n\
628 @example\n\
629 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
630 @end example\n\
631 \n\
632 The first form uses the unpivoted algorithm, which is faster, but less\n\
633 stable. The second form uses a slower pivoted algorithm, which is more\n\
634 stable.\n\
635 \n\
636 The matrix case is done as a sequence of rank-1 updates;\n\
637 thus, for large enough k, it will be both faster and more accurate to\n\
638 recompute the factorization from scratch.\n\
639 @seealso{lu, cholupdate, qrupdate}\n\
640 @end deftypefn")
641 {
642  octave_idx_type nargin = args.length ();
643  octave_value_list retval;
644 
645  bool pivoted = nargin == 5;
646 
647  if (nargin != 4 && nargin != 5)
648  {
649  print_usage ();
650  return retval;
651  }
652 
653  octave_value argl = args(0);
654  octave_value argu = args(1);
655  octave_value argp = pivoted ? args(2) : octave_value ();
656  octave_value argx = args(2 + pivoted);
657  octave_value argy = args(3 + pivoted);
658 
659  if (argl.is_numeric_type () && argu.is_numeric_type ()
660  && argx.is_numeric_type () && argy.is_numeric_type ()
661  && (! pivoted || argp.is_perm_matrix ()))
662  {
663  if (check_lu_dims (argl, argu, argp))
664  {
665  PermMatrix P = (pivoted
666  ? argp.perm_matrix_value ()
667  : PermMatrix::eye (argl.rows ()));
668 
669  if (argl.is_real_type ()
670  && argu.is_real_type ()
671  && argx.is_real_type ()
672  && argy.is_real_type ())
673  {
674  // all real case
675  if (argl.is_single_type ()
676  || argu.is_single_type ()
677  || argx.is_single_type ()
678  || argy.is_single_type ())
679  {
680  FloatMatrix L = argl.float_matrix_value ();
681  FloatMatrix U = argu.float_matrix_value ();
682  FloatMatrix x = argx.float_matrix_value ();
683  FloatMatrix y = argy.float_matrix_value ();
684 
685  FloatLU fact (L, U, P);
686  if (pivoted)
687  fact.update_piv (x, y);
688  else
689  fact.update (x, y);
690 
691  if (pivoted)
692  retval(2) = fact.P ();
693  retval(1) = get_lu_u (fact);
694  retval(0) = get_lu_l (fact);
695  }
696  else
697  {
698  Matrix L = argl.matrix_value ();
699  Matrix U = argu.matrix_value ();
700  Matrix x = argx.matrix_value ();
701  Matrix y = argy.matrix_value ();
702 
703  LU fact (L, U, P);
704  if (pivoted)
705  fact.update_piv (x, y);
706  else
707  fact.update (x, y);
708 
709  if (pivoted)
710  retval(2) = fact.P ();
711  retval(1) = get_lu_u (fact);
712  retval(0) = get_lu_l (fact);
713  }
714  }
715  else
716  {
717  // complex case
718  if (argl.is_single_type ()
719  || argu.is_single_type ()
720  || argx.is_single_type ()
721  || argy.is_single_type ())
722  {
723  FloatComplexMatrix L = argl.float_complex_matrix_value ();
724  FloatComplexMatrix U = argu.float_complex_matrix_value ();
725  FloatComplexMatrix x = argx.float_complex_matrix_value ();
726  FloatComplexMatrix y = argy.float_complex_matrix_value ();
727 
728  FloatComplexLU fact (L, U, P);
729  if (pivoted)
730  fact.update_piv (x, y);
731  else
732  fact.update (x, y);
733 
734  if (pivoted)
735  retval(2) = fact.P ();
736  retval(1) = get_lu_u (fact);
737  retval(0) = get_lu_l (fact);
738  }
739  else
740  {
741  ComplexMatrix L = argl.complex_matrix_value ();
742  ComplexMatrix U = argu.complex_matrix_value ();
743  ComplexMatrix x = argx.complex_matrix_value ();
744  ComplexMatrix y = argy.complex_matrix_value ();
745 
746  ComplexLU fact (L, U, P);
747  if (pivoted)
748  fact.update_piv (x, y);
749  else
750  fact.update (x, y);
751 
752  if (pivoted)
753  retval(2) = fact.P ();
754  retval(1) = get_lu_u (fact);
755  retval(0) = get_lu_l (fact);
756  }
757  }
758  }
759  else
760  error ("luupdate: dimension mismatch");
761  }
762  else
763  error ("luupdate: L, U, X, and Y must be numeric");
764 
765  return retval;
766 }
767 
768 /*
769 %!shared A, u, v, Ac, uc, vc
770 %! A = [0.091364 0.613038 0.999083;
771 %! 0.594638 0.425302 0.603537;
772 %! 0.383594 0.291238 0.085574;
773 %! 0.265712 0.268003 0.238409;
774 %! 0.669966 0.743851 0.445057 ];
775 %!
776 %! u = [0.85082;
777 %! 0.76426;
778 %! 0.42883;
779 %! 0.53010;
780 %! 0.80683 ];
781 %!
782 %! v = [0.98810;
783 %! 0.24295;
784 %! 0.43167 ];
785 %!
786 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
787 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
788 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
789 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
790 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
791 %!
792 %! uc = [0.20351 + 0.05401i;
793 %! 0.13141 + 0.43708i;
794 %! 0.29808 + 0.08789i;
795 %! 0.69821 + 0.38844i;
796 %! 0.74871 + 0.25821i ];
797 %!
798 %! vc = [0.85839 + 0.29468i;
799 %! 0.20820 + 0.93090i;
800 %! 0.86184 + 0.34689i ];
801 %!
802 
803 %!testif HAVE_QRUPDATE_LUU
804 %! [L,U,P] = lu (A);
805 %! [L,U] = luupdate (L,U,P*u,v);
806 %! assert (norm (vec (tril (L)-L), Inf) == 0);
807 %! assert (norm (vec (triu (U)-U), Inf) == 0);
808 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
809 %!
810 %!testif HAVE_QRUPDATE_LUU
811 %! [L,U,P] = lu (Ac);
812 %! [L,U] = luupdate (L,U,P*uc,vc);
813 %! assert (norm (vec (tril (L)-L), Inf) == 0);
814 %! assert (norm (vec (triu (U)-U), Inf) == 0);
815 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
816 
817 %!testif HAVE_QRUPDATE_LUU
818 %! [L,U,P] = lu (single (A));
819 %! [L,U] = luupdate (L,U,P*single (u), single (v));
820 %! assert (norm (vec (tril (L)-L), Inf) == 0);
821 %! assert (norm (vec (triu (U)-U), Inf) == 0);
822 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
823 %!
824 %!testif HAVE_QRUPDATE_LUU
825 %! [L,U,P] = lu (single (Ac));
826 %! [L,U] = luupdate (L,U,P*single (uc),single (vc));
827 %! assert (norm (vec (tril (L)-L), Inf) == 0);
828 %! assert (norm (vec (triu (U)-U), Inf) == 0);
829 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
830 
831 %!testif HAVE_QRUPDATE_LUU
832 %! [L,U,P] = lu (A);
833 %! [L,U,P] = luupdate (L,U,P,u,v);
834 %! assert (norm (vec (tril (L)-L), Inf) == 0);
835 %! assert (norm (vec (triu (U)-U), Inf) == 0);
836 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
837 %!
838 %!testif HAVE_QRUPDATE_LUU
839 %! [L,U,P] = lu (Ac);
840 %! [L,U,P] = luupdate (L,U,P,uc,vc);
841 %! assert (norm (vec (tril (L)-L), Inf) == 0);
842 %! assert (norm (vec (triu (U)-U), Inf) == 0);
843 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
844 
845 %!testif HAVE_QRUPDATE_LUU
846 %! [L,U,P] = lu (single (A));
847 %! [L,U,P] = luupdate (L,U,P,single (u),single (v));
848 %! assert (norm (vec (tril (L)-L), Inf) == 0);
849 %! assert (norm (vec (triu (U)-U), Inf) == 0);
850 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
851 %!
852 %!testif HAVE_QRUPDATE_LUU
853 %! [L,U,P] = lu (single (Ac));
854 %! [L,U,P] = luupdate (L,U,P,single (uc),single (vc));
855 %! assert (norm (vec (tril (L)-L), Inf) == 0);
856 %! assert (norm (vec (triu (U)-U), Inf) == 0);
857 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
858 */
Generated on Mon Dec 30 2013 03:04:27 for GNU Octave by   doxygen 1.8.1.2