GNU Octave  4.0.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-2015 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 ())
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 ())
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}.\n\
75 \n\
76 If @var{A} is full subroutines from @sc{lapack} are used and if @var{A} is\n\
77 sparse then @sc{umfpack} is used.\n\
78 \n\
79 The result is returned in a permuted form, according to the optional return\n\
80 value @var{P}. For example, given the matrix @code{a = [1, 2; 3, 4]},\n\
81 \n\
82 @example\n\
83 [l, u, p] = lu (@var{a})\n\
84 @end example\n\
85 \n\
86 @noindent\n\
87 returns\n\
88 \n\
89 @example\n\
90 @group\n\
91 l =\n\
92 \n\
93  1.00000 0.00000\n\
94  0.33333 1.00000\n\
95 \n\
96 u =\n\
97 \n\
98  3.00000 4.00000\n\
99  0.00000 0.66667\n\
100 \n\
101 p =\n\
102 \n\
103  0 1\n\
104  1 0\n\
105 @end group\n\
106 @end example\n\
107 \n\
108 The matrix is not required to be square.\n\
109 \n\
110 When called with two or three output arguments and a spare input matrix,\n\
111 @code{lu} does not attempt to perform sparsity preserving column\n\
112 permutations. Called with a fourth output argument, the sparsity\n\
113 preserving column transformation @var{Q} is returned, such that\n\
114 @code{@var{P} * @var{A} * @var{Q} = @var{L} * @var{U}}.\n\
115 \n\
116 Called with a fifth output argument and a sparse input matrix,\n\
117 @code{lu} attempts to use a scaling factor @var{R} on the input matrix\n\
118 such that\n\
119 @code{@var{P} * (@var{R} \\ @var{A}) * @var{Q} = @var{L} * @var{U}}.\n\
120 This typically leads to a sparser and more stable factorization.\n\
121 \n\
122 An additional input argument @var{thres}, that defines the pivoting\n\
123 threshold can be given. @var{thres} can be a scalar, in which case\n\
124 it defines the @sc{umfpack} pivoting tolerance for both symmetric and\n\
125 unsymmetric cases. If @var{thres} is a 2-element vector, then the first\n\
126 element defines the pivoting tolerance for the unsymmetric @sc{umfpack}\n\
127 pivoting strategy and the second for the symmetric strategy. By default,\n\
128 the values defined by @code{spparms} are used ([0.1, 0.001]).\n\
129 \n\
130 Given the string argument @qcode{\"vector\"}, @code{lu} returns the values\n\
131 of @var{P} and @var{Q} as vector values, such that for full matrix,\n\
132 @code{@var{A} (@var{P},:) = @var{L} * @var{U}}, and @code{@var{R}(@var{P},:)\n\
133 * @var{A} (:, @var{Q}) = @var{L} * @var{U}}.\n\
134 \n\
135 With two output arguments, returns the permuted forms of the upper and\n\
136 lower triangular matrices, such that @code{@var{A} = @var{L} * @var{U}}.\n\
137 With one output argument @var{y}, then the matrix returned by the @sc{lapack}\n\
138 routines is returned. If the input matrix is sparse then the matrix @var{L}\n\
139 is embedded into @var{U} to give a return value similar to the full case.\n\
140 For both full and sparse matrices, @code{lu} loses the permutation\n\
141 information.\n\
142 @seealso{luupdate, ilu, chol, hess, qr, qz, schur, svd}\n\
143 @end deftypefn")
144 {
145  octave_value_list retval;
146  int nargin = args.length ();
147  bool issparse = (nargin > 0 && args(0).is_sparse_type ());
148  bool scale = (nargout == 5);
149 
150  if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5))
151  || (!issparse && (nargin > 2 || nargout > 3)))
152  {
153  print_usage ();
154  return retval;
155  }
156 
157  bool vecout = false;
158  Matrix thres;
159 
160  int n = 1;
161  while (n < nargin && ! error_state)
162  {
163  if (args(n).is_string ())
164  {
165  std::string tmp = args(n++).string_value ();
166 
167  if (tmp.compare ("vector") == 0)
168  vecout = true;
169  else
170  error ("lu: unrecognized string argument");
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  if (arg.is_real_type ())
209  {
210 
212 
213  if (nargout < 4)
214  {
215 
216  ColumnVector Qinit;
217  Qinit.resize (nc);
218  for (octave_idx_type i = 0; i < nc; i++)
219  Qinit (i) = i;
220  SparseLU fact (m, Qinit, thres, false, true);
221 
222  if (nargout < 2)
223  retval(0) = fact.Y ();
224  else
225  {
226 
227  retval(1)
228  = octave_value (
229  fact.U () * fact.Pc_mat ().transpose (),
231  nc, fact.col_perm ()));
232 
233  PermMatrix P = fact.Pr_mat ();
234  SparseMatrix L = fact.L ();
235  if (nargout < 3)
236  retval(0)
237  = octave_value (P.transpose () * L,
239  nr, fact.row_perm ()));
240  else
241  {
242  retval(0) = L;
243  if (vecout)
244  retval(2) = fact.Pr_vec();
245  else
246  retval(2) = P;
247  }
248 
249  }
250 
251  }
252  else
253  {
254 
255  SparseLU fact (m, thres, scale);
256 
257  if (scale)
258  retval(4) = fact.R ();
259 
260  if (vecout)
261  {
262  retval(3) = fact.Pc_vec ();
263  retval(2) = fact.Pr_vec ();
264  }
265  else
266  {
267  retval(3) = fact.Pc_mat ();
268  retval(2) = fact.Pr_mat ();
269  }
270  retval(1) = octave_value (fact.U (),
272  retval(0) = octave_value (fact.L (),
274  }
275 
276  }
277  else if (arg.is_complex_type ())
278  {
280 
281  if (nargout < 4)
282  {
283 
284  ColumnVector Qinit;
285  Qinit.resize (nc);
286  for (octave_idx_type i = 0; i < nc; i++)
287  Qinit (i) = i;
288  SparseComplexLU fact (m, Qinit, thres, false, true);
289 
290  if (nargout < 2)
291 
292  retval(0) = fact.Y ();
293 
294  else
295  {
296 
297  retval(1)
298  = octave_value (
299  fact.U () * fact.Pc_mat ().transpose (),
301  nc, fact.col_perm ()));
302 
303  PermMatrix P = fact.Pr_mat ();
304  SparseComplexMatrix L = fact.L ();
305  if (nargout < 3)
306  retval(0)
307  = octave_value (P.transpose () * L,
309  nr, fact.row_perm ()));
310  else
311  {
312  retval(0) = L;
313  if (vecout)
314  retval(2) = fact.Pr_vec();
315  else
316  retval(2) = P;
317  }
318 
319  }
320 
321  }
322  else
323  {
324 
325  SparseComplexLU fact (m, thres, scale);
326 
327  if (scale)
328  retval(4) = fact.R ();
329 
330  if (vecout)
331  {
332  retval(3) = fact.Pc_vec ();
333  retval(2) = fact.Pr_vec ();
334  }
335  else
336  {
337  retval(3) = fact.Pc_mat ();
338  retval(2) = fact.Pr_mat ();
339  }
340  retval(1) = octave_value (fact.U (),
342  retval(0) = octave_value (fact.L (),
344  }
345 
346  }
347  else
348  gripe_wrong_type_arg ("lu", arg);
349  }
350  else
351  {
352  if (arg_is_empty < 0)
353  return retval;
354  else if (arg_is_empty > 0)
355  return octave_value_list (3, Matrix ());
356 
357  if (arg.is_real_type ())
358  {
359  if (arg.is_single_type ())
360  {
361  FloatMatrix m = arg.float_matrix_value ();
362 
363  if (! error_state)
364  {
365  FloatLU fact (m);
366 
367  switch (nargout)
368  {
369  case 0:
370  case 1:
371  retval(0) = fact.Y ();
372  break;
373 
374  case 2:
375  {
376  PermMatrix P = fact.P ();
377  FloatMatrix L = P.transpose () * fact.L ();
378  retval(1) = get_lu_u (fact);
379  retval(0) = L;
380  }
381  break;
382 
383  case 3:
384  default:
385  {
386  if (vecout)
387  retval(2) = fact.P_vec ();
388  else
389  retval(2) = fact.P ();
390  retval(1) = get_lu_u (fact);
391  retval(0) = get_lu_l (fact);
392  }
393  break;
394  }
395  }
396  }
397  else
398  {
399  Matrix m = arg.matrix_value ();
400 
401  if (! error_state)
402  {
403  LU fact (m);
404 
405  switch (nargout)
406  {
407  case 0:
408  case 1:
409  retval(0) = fact.Y ();
410  break;
411 
412  case 2:
413  {
414  PermMatrix P = fact.P ();
415  Matrix L = P.transpose () * fact.L ();
416  retval(1) = get_lu_u (fact);
417  retval(0) = L;
418  }
419  break;
420 
421  case 3:
422  default:
423  {
424  if (vecout)
425  retval(2) = fact.P_vec ();
426  else
427  retval(2) = fact.P ();
428  retval(1) = get_lu_u (fact);
429  retval(0) = get_lu_l (fact);
430  }
431  break;
432  }
433  }
434  }
435  }
436  else if (arg.is_complex_type ())
437  {
438  if (arg.is_single_type ())
439  {
441 
442  if (! error_state)
443  {
444  FloatComplexLU fact (m);
445 
446  switch (nargout)
447  {
448  case 0:
449  case 1:
450  retval(0) = fact.Y ();
451  break;
452 
453  case 2:
454  {
455  PermMatrix P = fact.P ();
456  FloatComplexMatrix L = P.transpose () * fact.L ();
457  retval(1) = get_lu_u (fact);
458  retval(0) = L;
459  }
460  break;
461 
462  case 3:
463  default:
464  {
465  if (vecout)
466  retval(2) = fact.P_vec ();
467  else
468  retval(2) = fact.P ();
469  retval(1) = get_lu_u (fact);
470  retval(0) = get_lu_l (fact);
471  }
472  break;
473  }
474  }
475  }
476  else
477  {
479 
480  if (! error_state)
481  {
482  ComplexLU fact (m);
483 
484  switch (nargout)
485  {
486  case 0:
487  case 1:
488  retval(0) = fact.Y ();
489  break;
490 
491  case 2:
492  {
493  PermMatrix P = fact.P ();
494  ComplexMatrix L = P.transpose () * fact.L ();
495  retval(1) = get_lu_u (fact);
496  retval(0) = L;
497  }
498  break;
499 
500  case 3:
501  default:
502  {
503  if (vecout)
504  retval(2) = fact.P_vec ();
505  else
506  retval(2) = fact.P ();
507  retval(1) = get_lu_u (fact);
508  retval(0) = get_lu_l (fact);
509  }
510  break;
511  }
512  }
513  }
514  }
515  else
516  gripe_wrong_type_arg ("lu", arg);
517  }
518 
519  return retval;
520 }
521 
522 /*
523 %!assert(lu ([1, 2; 3, 4]), [3, 4; 1/3, 2/3], eps);
524 
525 %!test
526 %! [l, u] = lu ([1, 2; 3, 4]);
527 %! assert (l, [1/3, 1; 1, 0], sqrt (eps));
528 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
529 
530 %!test
531 %! [l, u, p] = lu ([1, 2; 3, 4]);
532 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
533 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
534 %! assert (p(:,:), [0, 1; 1, 0], sqrt (eps));
535 
536 %!test
537 %! [l, u, p] = lu ([1, 2; 3, 4], "vector");
538 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
539 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
540 %! assert (p, [2;1], sqrt (eps));
541 
542 %!test
543 %! [l, u, p] = lu ([1, 2; 3, 4; 5, 6]);
544 %! assert (l, [1, 0; 1/5, 1; 3/5, 1/2], sqrt (eps));
545 %! assert (u, [5, 6; 0, 4/5], sqrt (eps));
546 %! assert (p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
547 
548 %!assert (lu (single ([1, 2; 3, 4])), single ([3, 4; 1/3, 2/3]), eps ("single"))
549 
550 %!test
551 %! [l, u] = lu (single ([1, 2; 3, 4]));
552 %! assert (l, single ([1/3, 1; 1, 0]), sqrt (eps ("single")));
553 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
554 
555 %!test
556 %! [l, u, p] = lu (single ([1, 2; 3, 4]));
557 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
558 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
559 %! assert (p(:,:), single ([0, 1; 1, 0]), sqrt (eps ("single")));
560 
561 %!test
562 %! [l, u, p] = lu (single ([1, 2; 3, 4]), "vector");
563 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
564 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
565 %! assert (p, single ([2;1]), sqrt (eps ("single")));
566 
567 %!test
568 %! [l u p] = lu (single ([1, 2; 3, 4; 5, 6]));
569 %! assert (l, single ([1, 0; 1/5, 1; 3/5, 1/2]), sqrt (eps ("single")));
570 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single")));
571 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single")));
572 
573 %!error lu ()
574 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2)
575 
576 %!testif HAVE_UMFPACK
577 %! Bi = [1 2 3 4 5 2 3 6 7 8 4 5 7 8 9];
578 %! Bj = [1 3 4 5 6 7 8 9 11 12 13 14 15 16 17];
579 %! Bv = [1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1];
580 %! B = sparse (Bi, Bj, Bv);
581 %! [L, U] = lu (B);
582 %! assert (L*U, B);
583 %! [L, U, P] = lu(B);
584 %! assert (P'*L*U, B);
585 %! [L, U, P, Q] = lu (B);
586 %! assert (P'*L*U*Q', B);
587 
588 */
589 
590 static
591 bool check_lu_dims (const octave_value& l, const octave_value& u,
592  const octave_value& p)
593 {
594  octave_idx_type m = l.rows ();
595  octave_idx_type k = u.rows ();
596  octave_idx_type n = u.columns ();
597  return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
598  && k == std::min (m, n)
599  && (p.is_undefined () || p.rows () == m));
600 }
601 
602 DEFUN (luupdate, args, ,
603  "-*- texinfo -*-\n\
604 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} luupdate (@var{L}, @var{U}, @var{x}, @var{y})\n\
605 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
606 Given an LU@tie{}factorization of a real or complex matrix\n\
607 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and\n\
608 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization\n\
609 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are\n\
610 column vectors (rank-1 update) or matrices with equal number of columns\n\
611 (rank-k update).\n\
612 \n\
613 Optionally, row-pivoted updating can be used by supplying a row permutation\n\
614 (pivoting) matrix @var{P}; in that case, an updated permutation matrix is\n\
615 returned. Note that if @var{L}, @var{U}, @var{P} is a pivoted\n\
616 LU@tie{}factorization as obtained by @code{lu}:\n\
617 \n\
618 @example\n\
619 [@var{L}, @var{U}, @var{P}] = lu (@var{A});\n\
620 @end example\n\
621 \n\
622 @noindent\n\
623 then a factorization of @tcode{@var{A}+@var{x}*@var{y}.'} can be obtained\n\
624 either as\n\
625 \n\
626 @example\n\
627 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})\n\
628 @end example\n\
629 \n\
630 @noindent\n\
631 or\n\
632 \n\
633 @example\n\
634 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
635 @end example\n\
636 \n\
637 The first form uses the unpivoted algorithm, which is faster, but less\n\
638 stable. The second form uses a slower pivoted algorithm, which is more\n\
639 stable.\n\
640 \n\
641 The matrix case is done as a sequence of rank-1 updates; thus, for large\n\
642 enough k, it will be both faster and more accurate to recompute the\n\
643 factorization from scratch.\n\
644 @seealso{lu, cholupdate, qrupdate}\n\
645 @end deftypefn")
646 {
647  octave_idx_type nargin = args.length ();
648  octave_value_list retval;
649 
650  bool pivoted = nargin == 5;
651 
652  if (nargin != 4 && nargin != 5)
653  {
654  print_usage ();
655  return retval;
656  }
657 
658  octave_value argl = args(0);
659  octave_value argu = args(1);
660  octave_value argp = pivoted ? args(2) : octave_value ();
661  octave_value argx = args(2 + pivoted);
662  octave_value argy = args(3 + pivoted);
663 
664  if (argl.is_numeric_type () && argu.is_numeric_type ()
665  && argx.is_numeric_type () && argy.is_numeric_type ()
666  && (! pivoted || argp.is_perm_matrix ()))
667  {
668  if (check_lu_dims (argl, argu, argp))
669  {
670  PermMatrix P = (pivoted
671  ? argp.perm_matrix_value ()
672  : PermMatrix::eye (argl.rows ()));
673 
674  if (argl.is_real_type ()
675  && argu.is_real_type ()
676  && argx.is_real_type ()
677  && argy.is_real_type ())
678  {
679  // all real case
680  if (argl.is_single_type ()
681  || argu.is_single_type ()
682  || argx.is_single_type ()
683  || argy.is_single_type ())
684  {
685  FloatMatrix L = argl.float_matrix_value ();
686  FloatMatrix U = argu.float_matrix_value ();
687  FloatMatrix x = argx.float_matrix_value ();
688  FloatMatrix y = argy.float_matrix_value ();
689 
690  FloatLU fact (L, U, P);
691  if (pivoted)
692  fact.update_piv (x, y);
693  else
694  fact.update (x, y);
695 
696  if (pivoted)
697  retval(2) = fact.P ();
698  retval(1) = get_lu_u (fact);
699  retval(0) = get_lu_l (fact);
700  }
701  else
702  {
703  Matrix L = argl.matrix_value ();
704  Matrix U = argu.matrix_value ();
705  Matrix x = argx.matrix_value ();
706  Matrix y = argy.matrix_value ();
707 
708  LU fact (L, U, P);
709  if (pivoted)
710  fact.update_piv (x, y);
711  else
712  fact.update (x, y);
713 
714  if (pivoted)
715  retval(2) = fact.P ();
716  retval(1) = get_lu_u (fact);
717  retval(0) = get_lu_l (fact);
718  }
719  }
720  else
721  {
722  // complex case
723  if (argl.is_single_type ()
724  || argu.is_single_type ()
725  || argx.is_single_type ()
726  || argy.is_single_type ())
727  {
732 
733  FloatComplexLU fact (L, U, P);
734  if (pivoted)
735  fact.update_piv (x, y);
736  else
737  fact.update (x, y);
738 
739  if (pivoted)
740  retval(2) = fact.P ();
741  retval(1) = get_lu_u (fact);
742  retval(0) = get_lu_l (fact);
743  }
744  else
745  {
750 
751  ComplexLU fact (L, U, P);
752  if (pivoted)
753  fact.update_piv (x, y);
754  else
755  fact.update (x, y);
756 
757  if (pivoted)
758  retval(2) = fact.P ();
759  retval(1) = get_lu_u (fact);
760  retval(0) = get_lu_l (fact);
761  }
762  }
763  }
764  else
765  error ("luupdate: dimension mismatch");
766  }
767  else
768  error ("luupdate: L, U, X, and Y must be numeric");
769 
770  return retval;
771 }
772 
773 /*
774 %!shared A, u, v, Ac, uc, vc
775 %! A = [0.091364 0.613038 0.999083;
776 %! 0.594638 0.425302 0.603537;
777 %! 0.383594 0.291238 0.085574;
778 %! 0.265712 0.268003 0.238409;
779 %! 0.669966 0.743851 0.445057 ];
780 %!
781 %! u = [0.85082;
782 %! 0.76426;
783 %! 0.42883;
784 %! 0.53010;
785 %! 0.80683 ];
786 %!
787 %! v = [0.98810;
788 %! 0.24295;
789 %! 0.43167 ];
790 %!
791 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
792 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
793 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
794 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
795 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
796 %!
797 %! uc = [0.20351 + 0.05401i;
798 %! 0.13141 + 0.43708i;
799 %! 0.29808 + 0.08789i;
800 %! 0.69821 + 0.38844i;
801 %! 0.74871 + 0.25821i ];
802 %!
803 %! vc = [0.85839 + 0.29468i;
804 %! 0.20820 + 0.93090i;
805 %! 0.86184 + 0.34689i ];
806 %!
807 
808 %!testif HAVE_QRUPDATE_LUU
809 %! [L,U,P] = lu (A);
810 %! [L,U] = luupdate (L,U,P*u,v);
811 %! assert (norm (vec (tril (L)-L), Inf) == 0);
812 %! assert (norm (vec (triu (U)-U), Inf) == 0);
813 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
814 %!
815 %!testif HAVE_QRUPDATE_LUU
816 %! [L,U,P] = lu (Ac);
817 %! [L,U] = luupdate (L,U,P*uc,vc);
818 %! assert (norm (vec (tril (L)-L), Inf) == 0);
819 %! assert (norm (vec (triu (U)-U), Inf) == 0);
820 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
821 
822 %!testif HAVE_QRUPDATE_LUU
823 %! [L,U,P] = lu (single (A));
824 %! [L,U] = luupdate (L,U,P*single (u), single (v));
825 %! assert (norm (vec (tril (L)-L), Inf) == 0);
826 %! assert (norm (vec (triu (U)-U), Inf) == 0);
827 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
828 %!
829 %!testif HAVE_QRUPDATE_LUU
830 %! [L,U,P] = lu (single (Ac));
831 %! [L,U] = luupdate (L,U,P*single (uc),single (vc));
832 %! assert (norm (vec (tril (L)-L), Inf) == 0);
833 %! assert (norm (vec (triu (U)-U), Inf) == 0);
834 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
835 
836 %!testif HAVE_QRUPDATE_LUU
837 %! [L,U,P] = lu (A);
838 %! [L,U,P] = luupdate (L,U,P,u,v);
839 %! assert (norm (vec (tril (L)-L), Inf) == 0);
840 %! assert (norm (vec (triu (U)-U), Inf) == 0);
841 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
842 %!
843 %!testif HAVE_QRUPDATE_LUU
844 %! [L,U,P] = lu (A);
845 %! [~,ordcols] = max (P,[],1);
846 %! [~,ordrows] = max (P,[],2);
847 %! P1 = eye (size(P))(:,ordcols);
848 %! P2 = eye (size(P))(ordrows,:);
849 %! assert(P1 == P);
850 %! assert(P2 == P);
851 %! [L,U,P] = luupdate (L,U,P,u,v);
852 %! [L,U,P1] = luupdate (L,U,P1,u,v);
853 %! [L,U,P2] = luupdate (L,U,P2,u,v);
854 %! assert(P1 == P);
855 %! assert(P2 == P);
856 %!
857 %!testif HAVE_QRUPDATE_LUU
858 %! [L,U,P] = lu (Ac);
859 %! [L,U,P] = luupdate (L,U,P,uc,vc);
860 %! assert (norm (vec (tril (L)-L), Inf) == 0);
861 %! assert (norm (vec (triu (U)-U), Inf) == 0);
862 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
863 
864 %!testif HAVE_QRUPDATE_LUU
865 %! [L,U,P] = lu (single (A));
866 %! [L,U,P] = luupdate (L,U,P,single (u),single (v));
867 %! assert (norm (vec (tril (L)-L), Inf) == 0);
868 %! assert (norm (vec (triu (U)-U), Inf) == 0);
869 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
870 %!
871 %!testif HAVE_QRUPDATE_LUU
872 %! [L,U,P] = lu (single (Ac));
873 %! [L,U,P] = luupdate (L,U,P,single (uc),single (vc));
874 %! assert (norm (vec (tril (L)-L), Inf) == 0);
875 %! assert (norm (vec (triu (U)-U), Inf) == 0);
876 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
877 */
void update(const ComplexColumnVector &u, const ComplexColumnVector &v)
Definition: CmplxLU.cc:88
bool regular(void) const
Definition: base-lu.cc:175
bool is_real_type(void) const
Definition: ov.h:651
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
Definition: gripes.cc:135
const octave_idx_type * col_perm(void) const
int ndims(void) const
Definition: ov.h:479
octave_idx_type rows(void) const
Definition: ov.h:473
lu_type U(void) const
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:130
const octave_idx_type * row_perm(void) const
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
bool is_numeric_type(void) const
Definition: ov.h:663
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
PermMatrix transpose(void) const
Definition: PermMatrix.cc:119
bool is_perm_matrix(void) const
Definition: ov.h:559
PermMatrix Pr_mat(void) const
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:246
void update(const ColumnVector &u, const ColumnVector &v)
Definition: dbleLU.cc:88
static bool check_lu_dims(const octave_value &l, const octave_value &u, const octave_value &p)
Definition: lu.cc:591
void update_piv(const ComplexColumnVector &u, const ComplexColumnVector &v)
Definition: CmplxLU.cc:139
PermMatrix P(void) const
Definition: base-lu.cc:152
SparseMatrix R(void) const
Definition: dbleLU.h:29
octave_idx_type nelem(void) const
Number of elements in the array.
Definition: Array.h:271
static octave_value get_lu_l(const base_lu< MT > &fact)
Definition: lu.cc:44
octave_idx_type columns(void) const
Definition: ov.h:475
ColumnVector Pr_vec(void) const
lu_type L(void) const
lu_type L(void) const
Definition: base-lu.cc:60
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
void update_piv(const FloatColumnVector &u, const FloatColumnVector &v)
Definition: floatLU.cc:139
Definition: dMatrix.h:35
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
void update(const FloatColumnVector &u, const FloatColumnVector &v)
Definition: floatLU.cc:88
static octave_value get_lu_u(const base_lu< MT > &fact)
Definition: lu.cc:55
double arg(double x)
Definition: lo-mappers.h:37
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
ColumnVector P_vec(void) const
Definition: base-lu.cc:159
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
void update_piv(const ColumnVector &u, const ColumnVector &v)
Definition: dbleLU.cc:138
PermMatrix Pc_mat(void) const
void update_piv(const FloatComplexColumnVector &u, const FloatComplexColumnVector &v)
Definition: fCmplxLU.cc:140
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5281
ColumnVector Pc_vec(void) const
bool is_undefined(void) const
Definition: ov.h:523
PermMatrix perm_matrix_value(void) const
Definition: ov.h:843
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:820
void update(const FloatComplexColumnVector &u, const FloatComplexColumnVector &v)
Definition: fCmplxLU.cc:88
bool is_single_type(void) const
Definition: ov.h:611
static PermMatrix eye(octave_idx_type n)
Definition: PermMatrix.cc:223
lu_type U(void) const
Definition: base-lu.cc:87
lu_type Y(void) const
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
F77_RET_T const double * x
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:210
void resize(octave_idx_type n, const double &rfv=0)
Definition: dColVector.h:106
lu_type Y(void) const
Definition: base-lu.cc:111