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