GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
lu.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1996-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include "lu.h"
31 #include "sparse-lu.h"
32 
33 #include "defun.h"
34 #include "error.h"
35 #include "errwarn.h"
36 #include "ovl.h"
37 #include "ov-re-sparse.h"
38 #include "ov-cx-sparse.h"
39 
40 template <typename MT>
41 static octave_value
43 {
44  MT L = fact.L ();
45  if (L.issquare ())
47  else
48  return L;
49 }
50 
51 template <typename MT>
52 static octave_value
54 {
55  MT U = fact.U ();
56  if (U.issquare () && fact.regular ())
58  else
59  return U;
60 }
61 
62 DEFUN (lu, args, nargout,
63  doc: /* -*- texinfo -*-
64 @deftypefn {} {[@var{L}, @var{U}] =} 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})
67 @deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}, @var{R}] =} lu (@var{S})
68 @deftypefnx {} {[@dots{}] =} lu (@var{S}, @var{thresh})
69 @deftypefnx {} {@var{y} =} lu (@dots{})
70 @deftypefnx {} {[@dots{}] =} lu (@dots{}, "vector")
71 @cindex LU decomposition
72 Compute the LU@tie{}decomposition of @var{A}.
73 
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.
76 
77 The result is returned in a permuted form, according to the optional return
78 value @var{P}. For example, given the matrix @code{@var{A} = [1, 2; 3, 4]},
79 
80 @example
81 [@var{L}, @var{U}, @var{P}] = lu (@var{A})
82 @end example
83 
84 @noindent
85 returns
86 
87 @example
88 @group
89 L =
90 
91  1.00000 0.00000
92  0.33333 1.00000
93 
94 U =
95 
96  3.00000 4.00000
97  0.00000 0.66667
98 
99 P =
100 
101  0 1
102  1 0
103 @end group
104 @end example
105 
106 The matrix is not required to be square.
107 
108 When called with two or three output arguments and a sparse input matrix,
109 @code{lu} does not attempt to perform sparsity preserving column permutations.
110 Called with a fourth output argument, the sparsity preserving column
111 transformation @var{Q} is returned, such that
112 @code{@var{P} * @var{A} * @var{Q} = @var{L} * @var{U}}. This is the
113 @strong{preferred} way to call @code{lu} with sparse input matrices.
114 
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
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.
119 
120 An additional input argument @var{thresh} that defines the pivoting
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
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}
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]).
127 
128 Given the string argument @qcode{"vector"}, @code{lu} returns the values
129 of @var{P} and @var{Q} as vector values, such that for full matrix,
130 @code{@var{A}(@var{P},:) = @var{L} * @var{U}}, and @code{@var{R}(@var{P},:)
131 * @var{A}(:,@var{Q}) = @var{L} * @var{U}}.
132 
133 With two output arguments, returns the permuted forms of the upper and
134 lower triangular matrices, such that @code{@var{A} = @var{L} * @var{U}}.
135 With one output argument @var{y}, then the matrix returned by the
136 @sc{lapack} routines is returned. If the input matrix is sparse then the
137 matrix @var{L} is embedded into @var{U} to give a return value similar to
138 the full case. For both full and sparse matrices, @code{lu} loses the
139 permutation information.
140 @seealso{luupdate, ilu, chol, hess, qr, qz, schur, svd}
141 @end deftypefn */)
142 {
143  int nargin = args.length ();
144  bool issparse = (nargin > 0 && args(0).issparse ());
145 
146  if (nargin < 1 || (issparse && nargin > 3) || (! issparse && nargin > 2))
147  print_usage ();
148 
149  bool vecout = false;
150  Matrix thresh;
151  int n = 1;
152 
153  while (n < nargin)
154  {
155  if (args(n).is_string ())
156  {
157  std::string tmp = args(n++).string_value ();
158 
159  if (tmp == "vector")
160  vecout = true;
161  else
162  error ("lu: unrecognized string argument");
163  }
164  else
165  {
166  if (! issparse)
167  error ("lu: can not define pivoting threshold THRESH for full matrices");
168 
169  Matrix tmp = args(n++).matrix_value ();
170  if (tmp.numel () == 1)
171  {
172  thresh.resize (1,2);
173  thresh(0) = tmp(0);
174  thresh(1) = tmp(0);
175  }
176  else if (tmp.numel () == 2)
177  thresh = tmp;
178  else
179  error ("lu: THRESH must be a 1- or 2-element vector");
180  }
181  }
182 
184  bool scale = (nargout == 5);
185 
186  octave_value arg = args(0);
187 
188  octave_idx_type nr = arg.rows ();
189  octave_idx_type nc = arg.columns ();
190 
191  if (issparse)
192  {
193  if (arg.isempty ())
194  return octave_value_list (5, SparseMatrix ());
195 
196  if (arg.isreal ())
197  {
199 
200  if (nargout < 4)
201  {
202  warning_with_id ("Octave:lu:sparse_input",
203  "lu: function may fail when called with less than 4 output arguments and a sparse input");
204 
205  ColumnVector Qinit (nc);
206  for (octave_idx_type i = 0; i < nc; i++)
207  Qinit(i) = i;
208  octave::math::sparse_lu<SparseMatrix> fact (m, Qinit, thresh,
209  false, true);
210 
211  if (nargout < 2)
212  retval(0) = fact.Y ();
213  else
214  {
215  retval.resize (nargout == 3 ? 3 : 2);
216  retval(1)
217  = octave_value (fact.U () * fact.Pc_mat ().transpose (),
219  nc, fact.col_perm ()));
220 
221  PermMatrix P = fact.Pr_mat ();
222  SparseMatrix L = fact.L ();
223 
224  if (nargout == 2)
225  retval(0)
226  = octave_value (P.transpose () * L,
228  nr, fact.row_perm ()));
229  else
230  {
231  retval(0) = L;
232  if (vecout)
233  retval(2) = fact.Pr_vec();
234  else
235  retval(2) = P;
236  }
237  }
238  }
239  else
240  {
241  retval.resize (scale ? 5 : 4);
243 
244  retval(0) = octave_value (fact.L (),
246  retval(1) = octave_value (fact.U (),
248 
249  if (vecout)
250  {
251  retval(2) = fact.Pr_vec ();
252  retval(3) = fact.Pc_vec ();
253  }
254  else
255  {
256  retval(2) = fact.Pr_mat ();
257  retval(3) = fact.Pc_mat ();
258  }
259 
260  if (scale)
261  retval(4) = fact.R ();
262  }
263  }
264  else if (arg.iscomplex ())
265  {
267 
268  if (nargout < 4)
269  {
270  warning_with_id ("Octave:lu:sparse_input",
271  "lu: function may fail when called with less than 4 output arguments and a sparse input");
272 
273  ColumnVector Qinit (nc);
274  for (octave_idx_type i = 0; i < nc; i++)
275  Qinit(i) = i;
277  thresh, false,
278  true);
279 
280  if (nargout < 2)
281  retval(0) = fact.Y ();
282  else
283  {
284  retval.resize (nargout == 3 ? 3 : 2);
285  retval(1)
286  = octave_value (fact.U () * fact.Pc_mat ().transpose (),
288  nc, fact.col_perm ()));
289 
290  PermMatrix P = fact.Pr_mat ();
291  SparseComplexMatrix L = fact.L ();
292  if (nargout == 2)
293  retval(0)
294  = octave_value (P.transpose () * L,
296  nr, fact.row_perm ()));
297  else
298  {
299  retval(0) = L;
300  if (vecout)
301  retval(2) = fact.Pr_vec();
302  else
303  retval(2) = P;
304  }
305  }
306  }
307  else
308  {
309  retval.resize (scale ? 5 : 4);
311  scale);
312 
313  retval(0) = octave_value (fact.L (),
315  retval(1) = octave_value (fact.U (),
317 
318  if (vecout)
319  {
320  retval(2) = fact.Pr_vec ();
321  retval(3) = fact.Pc_vec ();
322  }
323  else
324  {
325  retval(2) = fact.Pr_mat ();
326  retval(3) = fact.Pc_mat ();
327  }
328 
329  if (scale)
330  retval(4) = fact.R ();
331  }
332 
333  }
334  else
335  err_wrong_type_arg ("lu", arg);
336  }
337  else
338  {
339  if (arg.isempty ())
340  return octave_value_list (3, Matrix ());
341 
342  if (arg.isreal ())
343  {
344  if (arg.is_single_type ())
345  {
347 
349 
350  switch (nargout)
351  {
352  case 0:
353  case 1:
354  retval = ovl (fact.Y ());
355  break;
356 
357  case 2:
358  {
359  PermMatrix P = fact.P ();
360  FloatMatrix L = P.transpose () * fact.L ();
361  retval = ovl (L, get_lu_u (fact));
362  }
363  break;
364 
365  case 3:
366  default:
367  {
368  if (vecout)
369  retval = ovl (get_lu_l (fact), get_lu_u (fact),
370  fact.P_vec ());
371  else
372  retval = ovl (get_lu_l (fact), get_lu_u (fact),
373  fact.P ());
374  }
375  break;
376  }
377  }
378  else
379  {
380  Matrix m = arg.matrix_value ();
381 
383 
384  switch (nargout)
385  {
386  case 0:
387  case 1:
388  retval = ovl (fact.Y ());
389  break;
390 
391  case 2:
392  {
393  PermMatrix P = fact.P ();
394  Matrix L = P.transpose () * fact.L ();
395  retval = ovl (L, get_lu_u (fact));
396  }
397  break;
398 
399  case 3:
400  default:
401  {
402  if (vecout)
403  retval = ovl (get_lu_l (fact), get_lu_u (fact),
404  fact.P_vec ());
405  else
406  retval = ovl (get_lu_l (fact), get_lu_u (fact),
407  fact.P ());
408  }
409  break;
410  }
411  }
412  }
413  else if (arg.iscomplex ())
414  {
415  if (arg.is_single_type ())
416  {
418 
420 
421  switch (nargout)
422  {
423  case 0:
424  case 1:
425  retval = ovl (fact.Y ());
426  break;
427 
428  case 2:
429  {
430  PermMatrix P = fact.P ();
431  FloatComplexMatrix L = P.transpose () * fact.L ();
432  retval = ovl (L, get_lu_u (fact));
433  }
434  break;
435 
436  case 3:
437  default:
438  {
439  if (vecout)
440  retval = ovl (get_lu_l (fact), get_lu_u (fact),
441  fact.P_vec ());
442  else
443  retval = ovl (get_lu_l (fact), get_lu_u (fact),
444  fact.P ());
445  }
446  break;
447  }
448  }
449  else
450  {
452 
454 
455  switch (nargout)
456  {
457  case 0:
458  case 1:
459  retval = ovl (fact.Y ());
460  break;
461 
462  case 2:
463  {
464  PermMatrix P = fact.P ();
465  ComplexMatrix L = P.transpose () * fact.L ();
466  retval = ovl (L, get_lu_u (fact));
467  }
468  break;
469 
470  case 3:
471  default:
472  {
473  if (vecout)
474  retval = ovl (get_lu_l (fact), get_lu_u (fact),
475  fact.P_vec ());
476  else
477  retval = ovl (get_lu_l (fact), get_lu_u (fact),
478  fact.P ());
479  }
480  break;
481  }
482  }
483  }
484  else
485  err_wrong_type_arg ("lu", arg);
486  }
487 
488  return retval;
489 }
490 
491 /*
492 %!assert(lu ([1, 2; 3, 4]), [3, 4; 1/3, 2/3], eps)
493 
494 %!test
495 %! [l, u] = lu ([1, 2; 3, 4]);
496 %! assert (l, [1/3, 1; 1, 0], sqrt (eps));
497 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
498 
499 %!test
500 %! [l, u, p] = lu ([1, 2; 3, 4]);
501 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
502 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
503 %! assert (p(:,:), [0, 1; 1, 0], sqrt (eps));
504 
505 %!test
506 %! [l, u, p] = lu ([1, 2; 3, 4], "vector");
507 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
508 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
509 %! assert (p, [2;1], sqrt (eps));
510 
511 %!test
512 %! [l, u, p] = lu ([1, 2; 3, 4; 5, 6]);
513 %! assert (l, [1, 0; 1/5, 1; 3/5, 1/2], sqrt (eps));
514 %! assert (u, [5, 6; 0, 4/5], sqrt (eps));
515 %! assert (p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
516 
517 %!assert (lu (single ([1, 2; 3, 4])), single ([3, 4; 1/3, 2/3]), eps ("single"))
518 
519 %!test
520 %! [l, u] = lu (single ([1, 2; 3, 4]));
521 %! assert (l, single ([1/3, 1; 1, 0]), sqrt (eps ("single")));
522 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
523 
524 %!test
525 %! [l, u, p] = lu (single ([1, 2; 3, 4]));
526 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
527 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
528 %! assert (p(:,:), single ([0, 1; 1, 0]), sqrt (eps ("single")));
529 
530 %!test
531 %! [l, u, p] = lu (single ([1, 2; 3, 4]), "vector");
532 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
533 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
534 %! assert (p, single ([2;1]), sqrt (eps ("single")));
535 
536 %!test
537 %! [l u p] = lu (single ([1, 2; 3, 4; 5, 6]));
538 %! assert (l, single ([1, 0; 1/5, 1; 3/5, 1/2]), sqrt (eps ("single")));
539 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single")));
540 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single")));
541 
542 %!testif HAVE_UMFPACK
543 %! Bi = [1 2 3 4 5 2 3 6 7 8 4 5 7 8 9];
544 %! Bj = [1 3 4 5 6 7 8 9 11 12 13 14 15 16 17];
545 %! Bv = [1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1];
546 %! B = sparse (Bi, Bj, Bv);
547 %! warning ("off", "Octave:lu:sparse_input", "local");
548 %! [L, U] = lu (B);
549 %! assert (L*U, B);
550 %! [L, U, P] = lu(B);
551 %! assert (P'*L*U, B);
552 %! [L, U, P, Q] = lu (B);
553 %! assert (P'*L*U*Q', B);
554 
555 %!error lu ()
556 %!testif HAVE_UMFPACK
557 %! fail ("[l,u] = lu (sparse (magic (3)))", "warning", "function may fail");
558 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2)
559 
560 */
561 
562 static
563 bool check_lu_dims (const octave_value& l, const octave_value& u,
564  const octave_value& p)
565 {
566  octave_idx_type m = l.rows ();
567  octave_idx_type k = u.rows ();
568  octave_idx_type n = u.columns ();
569 
570  return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
571  && k == std::min (m, n)
572  && (p.is_undefined () || p.rows () == m));
573 }
574 
575 DEFUN (luupdate, args, ,
576  doc: /* -*- texinfo -*-
577 @deftypefn {} {[@var{L}, @var{U}] =} luupdate (@var{L}, @var{U}, @var{x}, @var{y})
578 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})
579 Given an LU@tie{}factorization of a real or complex matrix
580 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and
581 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization
582 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are
583 column vectors (rank-1 update) or matrices with equal number of columns
584 (rank-k update).
585 
586 Optionally, row-pivoted updating can be used by supplying a row permutation
587 (pivoting) matrix @var{P}; in that case, an updated permutation matrix is
588 returned. Note that if @var{L}, @var{U}, @var{P} is a pivoted
589 LU@tie{}factorization as obtained by @code{lu}:
590 
591 @example
592 [@var{L}, @var{U}, @var{P}] = lu (@var{A});
593 @end example
594 
595 @noindent
596 then a factorization of @tcode{@var{A}+@var{x}*@var{y}.'} can be obtained
597 either as
598 
599 @example
600 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})
601 @end example
602 
603 @noindent
604 or
605 
606 @example
607 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})
608 @end example
609 
610 The first form uses the unpivoted algorithm, which is faster, but less
611 stable. The second form uses a slower pivoted algorithm, which is more
612 stable.
613 
614 The matrix case is done as a sequence of rank-1 updates; thus, for large
615 enough k, it will be both faster and more accurate to recompute the
616 factorization from scratch.
617 @seealso{lu, cholupdate, qrupdate}
618 @end deftypefn */)
619 {
620  int nargin = args.length ();
621 
622  if (nargin != 4 && nargin != 5)
623  print_usage ();
624 
625  bool pivoted = (nargin == 5);
626 
627  octave_value argl = args(0);
628  octave_value argu = args(1);
629  octave_value argp = (pivoted ? args(2) : octave_value ());
630  octave_value argx = args(2 + pivoted);
631  octave_value argy = args(3 + pivoted);
632 
633  if (! (argl.isnumeric () && argu.isnumeric ()
634  && argx.isnumeric () && argy.isnumeric ()
635  && (! pivoted || argp.is_perm_matrix ())))
636  error ("luupdate: L, U, X, and Y must be numeric");
637 
638  if (! check_lu_dims (argl, argu, argp))
639  error ("luupdate: dimension mismatch");
640 
641  PermMatrix P = (pivoted
642  ? argp.perm_matrix_value ()
643  : PermMatrix::eye (argl.rows ()));
644 
645  if (argl.isreal () && argu.isreal ()
646  && argx.isreal () && argy.isreal ())
647  {
648  // all real case
649  if (argl.is_single_type () || argu.is_single_type ()
650  || argx.is_single_type () || argy.is_single_type ())
651  {
652  FloatMatrix L = argl.float_matrix_value ();
653  FloatMatrix U = argu.float_matrix_value ();
654  FloatMatrix x = argx.float_matrix_value ();
655  FloatMatrix y = argy.float_matrix_value ();
656 
657  octave::math::lu<FloatMatrix> fact (L, U, P);
658  if (pivoted)
659  fact.update_piv (x, y);
660  else
661  fact.update (x, y);
662 
663  if (pivoted)
664  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
665  else
666  return ovl (get_lu_l (fact), get_lu_u (fact));
667  }
668  else
669  {
670  Matrix L = argl.matrix_value ();
671  Matrix U = argu.matrix_value ();
672  Matrix x = argx.matrix_value ();
673  Matrix y = argy.matrix_value ();
674 
675  octave::math::lu<Matrix> fact (L, U, P);
676  if (pivoted)
677  fact.update_piv (x, y);
678  else
679  fact.update (x, y);
680 
681  if (pivoted)
682  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
683  else
684  return ovl (get_lu_l (fact), get_lu_u (fact));
685  }
686  }
687  else
688  {
689  // complex case
690  if (argl.is_single_type () || argu.is_single_type ()
691  || argx.is_single_type () || argy.is_single_type ())
692  {
697 
699  if (pivoted)
700  fact.update_piv (x, y);
701  else
702  fact.update (x, y);
703 
704  if (pivoted)
705  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
706  else
707  return ovl (get_lu_l (fact), get_lu_u (fact));
708  }
709  else
710  {
715 
716  octave::math::lu<ComplexMatrix> fact (L, U, P);
717  if (pivoted)
718  fact.update_piv (x, y);
719  else
720  fact.update (x, y);
721 
722  if (pivoted)
723  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
724  else
725  return ovl (get_lu_l (fact), get_lu_u (fact));
726  }
727  }
728 }
729 
730 /*
731 %!shared A, u, v, Ac, uc, vc
732 %! A = [0.091364 0.613038 0.999083;
733 %! 0.594638 0.425302 0.603537;
734 %! 0.383594 0.291238 0.085574;
735 %! 0.265712 0.268003 0.238409;
736 %! 0.669966 0.743851 0.445057 ];
737 %!
738 %! u = [0.85082;
739 %! 0.76426;
740 %! 0.42883;
741 %! 0.53010;
742 %! 0.80683 ];
743 %!
744 %! v = [0.98810;
745 %! 0.24295;
746 %! 0.43167 ];
747 %!
748 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
749 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
750 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
751 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
752 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
753 %!
754 %! uc = [0.20351 + 0.05401i;
755 %! 0.13141 + 0.43708i;
756 %! 0.29808 + 0.08789i;
757 %! 0.69821 + 0.38844i;
758 %! 0.74871 + 0.25821i ];
759 %!
760 %! vc = [0.85839 + 0.29468i;
761 %! 0.20820 + 0.93090i;
762 %! 0.86184 + 0.34689i ];
763 %!
764 
765 %!testif HAVE_QRUPDATE_LUU
766 %! [L,U,P] = lu (A);
767 %! [L,U] = luupdate (L,U,P*u,v);
768 %! assert (norm (vec (tril (L)-L), Inf) == 0);
769 %! assert (norm (vec (triu (U)-U), Inf) == 0);
770 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
771 %!
772 %!testif HAVE_QRUPDATE_LUU
773 %! [L,U,P] = lu (Ac);
774 %! [L,U] = luupdate (L,U,P*uc,vc);
775 %! assert (norm (vec (tril (L)-L), Inf) == 0);
776 %! assert (norm (vec (triu (U)-U), Inf) == 0);
777 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
778 
779 %!testif HAVE_QRUPDATE_LUU
780 %! [L,U,P] = lu (single (A));
781 %! [L,U] = luupdate (L,U,P*single (u), single (v));
782 %! assert (norm (vec (tril (L)-L), Inf) == 0);
783 %! assert (norm (vec (triu (U)-U), Inf) == 0);
784 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
785 %!
786 %!testif HAVE_QRUPDATE_LUU
787 %! [L,U,P] = lu (single (Ac));
788 %! [L,U] = luupdate (L,U,P*single (uc),single (vc));
789 %! assert (norm (vec (tril (L)-L), Inf) == 0);
790 %! assert (norm (vec (triu (U)-U), Inf) == 0);
791 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
792 
793 %!testif HAVE_QRUPDATE_LUU
794 %! [L,U,P] = lu (A);
795 %! [L,U,P] = luupdate (L,U,P,u,v);
796 %! assert (norm (vec (tril (L)-L), Inf) == 0);
797 %! assert (norm (vec (triu (U)-U), Inf) == 0);
798 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
799 %!
800 %!testif HAVE_QRUPDATE_LUU
801 %! [L,U,P] = lu (A);
802 %! [~,ordcols] = max (P,[],1);
803 %! [~,ordrows] = max (P,[],2);
804 %! P1 = eye (size(P))(:,ordcols);
805 %! P2 = eye (size(P))(ordrows,:);
806 %! assert(P1 == P);
807 %! assert(P2 == P);
808 %! [L,U,P] = luupdate (L,U,P,u,v);
809 %! [L,U,P1] = luupdate (L,U,P1,u,v);
810 %! [L,U,P2] = luupdate (L,U,P2,u,v);
811 %! assert(P1 == P);
812 %! assert(P2 == P);
813 %!
814 %!testif HAVE_QRUPDATE_LUU
815 %! [L,U,P] = lu (Ac);
816 %! [L,U,P] = luupdate (L,U,P,uc,vc);
817 %! assert (norm (vec (tril (L)-L), Inf) == 0);
818 %! assert (norm (vec (triu (U)-U), Inf) == 0);
819 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
820 
821 %!testif HAVE_QRUPDATE_LUU
822 %! [L,U,P] = lu (single (A));
823 %! [L,U,P] = luupdate (L,U,P,single (u),single (v));
824 %! assert (norm (vec (tril (L)-L), Inf) == 0);
825 %! assert (norm (vec (triu (U)-U), Inf) == 0);
826 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
827 %!
828 %!testif HAVE_QRUPDATE_LUU
829 %! [L,U,P] = lu (single (Ac));
830 %! [L,U,P] = luupdate (L,U,P,single (uc),single (vc));
831 %! assert (norm (vec (tril (L)-L), Inf) == 0);
832 %! assert (norm (vec (triu (U)-U), Inf) == 0);
833 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
834 */
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:207
void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array.cc:1011
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:377
@ Permuted_Lower
Definition: MatrixType.h:54
@ Permuted_Upper
Definition: MatrixType.h:53
Definition: dMatrix.h:42
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:151
static PermMatrix eye(octave_idx_type n)
Definition: PermMatrix.cc:205
PermMatrix transpose(void) const
Definition: PermMatrix.cc:101
Definition: mx-defs.h:80
T U(void) const
Definition: lu.cc:123
bool regular(void) const
Definition: lu.cc:212
PermMatrix P(void) const
Definition: lu.cc:189
void update(const VT &u, const VT &v)
T Y(void) const
Definition: lu.cc:147
void update_piv(const VT &u, const VT &v)
ColumnVector P_vec(void) const
Definition: lu.cc:196
T L(void) const
Definition: lu.cc:96
lu_type U(void) const
Definition: sparse-lu.h:89
PermMatrix Pc_mat(void) const
Definition: sparse-lu.cc:978
ColumnVector Pc_vec(void) const
Definition: sparse-lu.cc:964
PermMatrix Pr_mat(void) const
Definition: sparse-lu.cc:937
const octave_idx_type * row_perm(void) const
Definition: sparse-lu.h:107
ColumnVector Pr_vec(void) const
Definition: sparse-lu.cc:923
const octave_idx_type * col_perm(void) const
Definition: sparse-lu.h:109
lu_type L(void) const
Definition: sparse-lu.h:87
SparseMatrix R(void) const
Definition: sparse-lu.h:91
lu_type Y(void) const
Definition: sparse-lu.cc:866
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:853
bool isreal(void) const
Definition: ov.h:691
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:824
octave_idx_type rows(void) const
Definition: ov.h:504
bool isnumeric(void) const
Definition: ov.h:703
octave_idx_type columns(void) const
Definition: ov.h:506
int ndims(void) const
Definition: ov.h:510
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:809
PermMatrix perm_matrix_value(void) const
Definition: ov.h:876
bool isempty(void) const
Definition: ov.h:557
bool is_single_type(void) const
Definition: ov.h:651
bool is_undefined(void) const
Definition: ov.h:554
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:806
bool is_perm_matrix(void) const
Definition: ov.h:590
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:828
bool iscomplex(void) const
Definition: ov.h:694
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:857
OCTINTERP_API void print_usage(void)
Definition: defun.cc:53
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
void warning_with_id(const char *id, const char *fmt,...)
Definition: error.cc:1065
void error(const char *fmt,...)
Definition: error.cc:968
void err_wrong_type_arg(const char *name, const char *s)
Definition: errwarn.cc:166
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5846
static bool check_lu_dims(const octave_value &l, const octave_value &u, const octave_value &p)
Definition: lu.cc:563
static octave_value get_lu_u(const octave::math::lu< MT > &fact)
Definition: lu.cc:53
static octave_value get_lu_l(const octave::math::lu< MT > &fact)
Definition: lu.cc:42
F77_RET_T const F77_DBLE * x
T octave_idx_type m
Definition: mx-inlines.cc:773
octave_idx_type n
Definition: mx-inlines.cc:753
return octave_value(v1.char_array_value() . concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string()) ? '\'' :'"'))
octave_value::octave_value(const Array< char > &chm, char type) return retval
Definition: ov.cc:811
octave_value_list ovl(const OV_Args &... args)
Construct an octave_value_list with less typing.
Definition: ovl.h:211