GNU Octave  4.2.1
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-2017 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 #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.is_square ())
44  else
45  return L;
46 }
47 
48 template <typename MT>
49 static octave_value
51 {
52  MT U = fact.U ();
53  if (U.is_square () && 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
107 permutations. Called with a fourth output argument, the sparsity
108 preserving column transformation @var{Q} is returned, such that
109 @code{@var{P} * @var{A} * @var{Q} = @var{L} * @var{U}}.
110 
111 Called with a fifth output argument and a sparse input matrix,
112 @code{lu} attempts to use a scaling factor @var{R} on the input matrix
113 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).is_sparse_type ());
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.is_empty ())
191  return octave_value_list (5, SparseMatrix ());
192 
193  if (arg.is_real_type ())
194  {
196 
197  if (nargout < 4)
198  {
199  ColumnVector Qinit (nc);
200  for (octave_idx_type i = 0; i < nc; i++)
201  Qinit(i) = i;
202  octave::math::sparse_lu<SparseMatrix> fact (m, Qinit, thres, false, true);
203 
204  if (nargout < 2)
205  retval(0) = fact.Y ();
206  else
207  {
208  retval.resize (nargout == 3 ? 3 : 2);
209  retval(1)
210  = octave_value (
211  fact.U () * fact.Pc_mat ().transpose (),
213  nc, fact.col_perm ()));
214 
215  PermMatrix P = fact.Pr_mat ();
216  SparseMatrix L = fact.L ();
217 
218  if (nargout == 2)
219  retval(0)
220  = octave_value (P.transpose () * L,
222  nr, fact.row_perm ()));
223  else
224  {
225  retval(0) = L;
226  if (vecout)
227  retval(2) = fact.Pr_vec();
228  else
229  retval(2) = P;
230  }
231  }
232  }
233  else
234  {
235  retval.resize (scale ? 5 : 4);
236  octave::math::sparse_lu<SparseMatrix> fact (m, thres, scale);
237 
238  retval(0) = octave_value (fact.L (),
240  retval(1) = octave_value (fact.U (),
242 
243  if (vecout)
244  {
245  retval(2) = fact.Pr_vec ();
246  retval(3) = fact.Pc_vec ();
247  }
248  else
249  {
250  retval(2) = fact.Pr_mat ();
251  retval(3) = fact.Pc_mat ();
252  }
253 
254  if (scale)
255  retval(4) = fact.R ();
256  }
257  }
258  else if (arg.is_complex_type ())
259  {
261 
262  if (nargout < 4)
263  {
264  ColumnVector Qinit (nc);
265  for (octave_idx_type i = 0; i < nc; i++)
266  Qinit(i) = i;
268  thres, false,
269  true);
270 
271  if (nargout < 2)
272  retval(0) = fact.Y ();
273  else
274  {
275  retval.resize (nargout == 3 ? 3 : 2);
276  retval(1)
277  = octave_value (
278  fact.U () * fact.Pc_mat ().transpose (),
280  nc, fact.col_perm ()));
281 
282  PermMatrix P = fact.Pr_mat ();
283  SparseComplexMatrix L = fact.L ();
284  if (nargout == 2)
285  retval(0)
286  = octave_value (P.transpose () * L,
288  nr, fact.row_perm ()));
289  else
290  {
291  retval(0) = L;
292  if (vecout)
293  retval(2) = fact.Pr_vec();
294  else
295  retval(2) = P;
296  }
297  }
298  }
299  else
300  {
301  retval.resize (scale ? 5 : 4);
302  octave::math::sparse_lu<SparseComplexMatrix> fact (m, thres, scale);
303 
304  retval(0) = octave_value (fact.L (),
306  retval(1) = octave_value (fact.U (),
308 
309  if (vecout)
310  {
311  retval(2) = fact.Pr_vec ();
312  retval(3) = fact.Pc_vec ();
313  }
314  else
315  {
316  retval(2) = fact.Pr_mat ();
317  retval(3) = fact.Pc_mat ();
318  }
319 
320  if (scale)
321  retval(4) = fact.R ();
322  }
323 
324  }
325  else
326  err_wrong_type_arg ("lu", arg);
327  }
328  else
329  {
330  if (arg.is_empty ())
331  return octave_value_list (3, Matrix ());
332 
333  if (arg.is_real_type ())
334  {
335  if (arg.is_single_type ())
336  {
338 
340 
341  switch (nargout)
342  {
343  case 0:
344  case 1:
345  retval = ovl (fact.Y ());
346  break;
347 
348  case 2:
349  {
350  PermMatrix P = fact.P ();
351  FloatMatrix L = P.transpose () * fact.L ();
352  retval = ovl (L, get_lu_u (fact));
353  }
354  break;
355 
356  case 3:
357  default:
358  {
359  if (vecout)
360  retval = ovl (get_lu_l (fact), get_lu_u (fact),
361  fact.P_vec ());
362  else
363  retval = ovl (get_lu_l (fact), get_lu_u (fact),
364  fact.P ());
365  }
366  break;
367  }
368  }
369  else
370  {
371  Matrix m = arg.matrix_value ();
372 
373  octave::math::lu<Matrix> fact (m);
374 
375  switch (nargout)
376  {
377  case 0:
378  case 1:
379  retval = ovl (fact.Y ());
380  break;
381 
382  case 2:
383  {
384  PermMatrix P = fact.P ();
385  Matrix L = P.transpose () * fact.L ();
386  retval = ovl (L, get_lu_u (fact));
387  }
388  break;
389 
390  case 3:
391  default:
392  {
393  if (vecout)
394  retval = ovl (get_lu_l (fact), get_lu_u (fact),
395  fact.P_vec ());
396  else
397  retval = ovl (get_lu_l (fact), get_lu_u (fact),
398  fact.P ());
399  }
400  break;
401  }
402  }
403  }
404  else if (arg.is_complex_type ())
405  {
406  if (arg.is_single_type ())
407  {
409 
411 
412  switch (nargout)
413  {
414  case 0:
415  case 1:
416  retval = ovl (fact.Y ());
417  break;
418 
419  case 2:
420  {
421  PermMatrix P = fact.P ();
422  FloatComplexMatrix L = P.transpose () * fact.L ();
423  retval = ovl (L, get_lu_u (fact));
424  }
425  break;
426 
427  case 3:
428  default:
429  {
430  if (vecout)
431  retval = ovl (get_lu_l (fact), get_lu_u (fact),
432  fact.P_vec ());
433  else
434  retval = ovl (get_lu_l (fact), get_lu_u (fact),
435  fact.P ());
436  }
437  break;
438  }
439  }
440  else
441  {
443 
445 
446  switch (nargout)
447  {
448  case 0:
449  case 1:
450  retval = ovl (fact.Y ());
451  break;
452 
453  case 2:
454  {
455  PermMatrix P = fact.P ();
456  ComplexMatrix L = P.transpose () * fact.L ();
457  retval = ovl (L, get_lu_u (fact));
458  }
459  break;
460 
461  case 3:
462  default:
463  {
464  if (vecout)
465  retval = ovl (get_lu_l (fact), get_lu_u (fact),
466  fact.P_vec ());
467  else
468  retval = ovl (get_lu_l (fact), get_lu_u (fact),
469  fact.P ());
470  }
471  break;
472  }
473  }
474  }
475  else
476  err_wrong_type_arg ("lu", arg);
477  }
478 
479  return retval;
480 }
481 
482 /*
483 %!assert(lu ([1, 2; 3, 4]), [3, 4; 1/3, 2/3], eps)
484 
485 %!test
486 %! [l, u] = lu ([1, 2; 3, 4]);
487 %! assert (l, [1/3, 1; 1, 0], sqrt (eps));
488 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
489 
490 %!test
491 %! [l, u, p] = lu ([1, 2; 3, 4]);
492 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
493 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
494 %! assert (p(:,:), [0, 1; 1, 0], sqrt (eps));
495 
496 %!test
497 %! [l, u, p] = lu ([1, 2; 3, 4], "vector");
498 %! assert (l, [1, 0; 1/3, 1], sqrt (eps));
499 %! assert (u, [3, 4; 0, 2/3], sqrt (eps));
500 %! assert (p, [2;1], sqrt (eps));
501 
502 %!test
503 %! [l, u, p] = lu ([1, 2; 3, 4; 5, 6]);
504 %! assert (l, [1, 0; 1/5, 1; 3/5, 1/2], sqrt (eps));
505 %! assert (u, [5, 6; 0, 4/5], sqrt (eps));
506 %! assert (p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
507 
508 %!assert (lu (single ([1, 2; 3, 4])), single ([3, 4; 1/3, 2/3]), eps ("single"))
509 
510 %!test
511 %! [l, u] = lu (single ([1, 2; 3, 4]));
512 %! assert (l, single ([1/3, 1; 1, 0]), sqrt (eps ("single")));
513 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
514 
515 %!test
516 %! [l, u, p] = lu (single ([1, 2; 3, 4]));
517 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
518 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
519 %! assert (p(:,:), single ([0, 1; 1, 0]), sqrt (eps ("single")));
520 
521 %!test
522 %! [l, u, p] = lu (single ([1, 2; 3, 4]), "vector");
523 %! assert (l, single ([1, 0; 1/3, 1]), sqrt (eps ("single")));
524 %! assert (u, single ([3, 4; 0, 2/3]), sqrt (eps ("single")));
525 %! assert (p, single ([2;1]), sqrt (eps ("single")));
526 
527 %!test
528 %! [l u p] = lu (single ([1, 2; 3, 4; 5, 6]));
529 %! assert (l, single ([1, 0; 1/5, 1; 3/5, 1/2]), sqrt (eps ("single")));
530 %! assert (u, single ([5, 6; 0, 4/5]), sqrt (eps ("single")));
531 %! assert (p(:,:), single ([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps ("single")));
532 
533 %!error lu ()
534 %!error <can not define pivoting threshold> lu ([1, 2; 3, 4], 2)
535 
536 %!testif HAVE_UMFPACK
537 %! Bi = [1 2 3 4 5 2 3 6 7 8 4 5 7 8 9];
538 %! Bj = [1 3 4 5 6 7 8 9 11 12 13 14 15 16 17];
539 %! Bv = [1 1 1 1 1 1 -1 1 1 1 1 -1 1 -1 1];
540 %! B = sparse (Bi, Bj, Bv);
541 %! [L, U] = lu (B);
542 %! assert (L*U, B);
543 %! [L, U, P] = lu(B);
544 %! assert (P'*L*U, B);
545 %! [L, U, P, Q] = lu (B);
546 %! assert (P'*L*U*Q', B);
547 
548 */
549 
550 static
551 bool check_lu_dims (const octave_value& l, const octave_value& u,
552  const octave_value& p)
553 {
554  octave_idx_type m = l.rows ();
555  octave_idx_type k = u.rows ();
556  octave_idx_type n = u.columns ();
557 
558  return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
559  && k == std::min (m, n)
560  && (p.is_undefined () || p.rows () == m));
561 }
562 
563 DEFUN (luupdate, args, ,
564  doc: /* -*- texinfo -*-
565 @deftypefn {} {[@var{L}, @var{U}] =} luupdate (@var{L}, @var{U}, @var{x}, @var{y})
566 @deftypefnx {} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})
567 Given an LU@tie{}factorization of a real or complex matrix
568 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and
569 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization
570 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are
571 column vectors (rank-1 update) or matrices with equal number of columns
572 (rank-k update).
573 
574 Optionally, row-pivoted updating can be used by supplying a row permutation
575 (pivoting) matrix @var{P}; in that case, an updated permutation matrix is
576 returned. Note that if @var{L}, @var{U}, @var{P} is a pivoted
577 LU@tie{}factorization as obtained by @code{lu}:
578 
579 @example
580 [@var{L}, @var{U}, @var{P}] = lu (@var{A});
581 @end example
582 
583 @noindent
584 then a factorization of @tcode{@var{A}+@var{x}*@var{y}.'} can be obtained
585 either as
586 
587 @example
588 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})
589 @end example
590 
591 @noindent
592 or
593 
594 @example
595 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})
596 @end example
597 
598 The first form uses the unpivoted algorithm, which is faster, but less
599 stable. The second form uses a slower pivoted algorithm, which is more
600 stable.
601 
602 The matrix case is done as a sequence of rank-1 updates; thus, for large
603 enough k, it will be both faster and more accurate to recompute the
604 factorization from scratch.
605 @seealso{lu, cholupdate, qrupdate}
606 @end deftypefn */)
607 {
608  int nargin = args.length ();
609 
610  if (nargin != 4 && nargin != 5)
611  print_usage ();
612 
613  bool pivoted = (nargin == 5);
614 
615  octave_value argl = args(0);
616  octave_value argu = args(1);
617  octave_value argp = pivoted ? args(2) : octave_value ();
618  octave_value argx = args(2 + pivoted);
619  octave_value argy = args(3 + pivoted);
620 
621  if (! (argl.is_numeric_type () && argu.is_numeric_type ()
622  && argx.is_numeric_type () && argy.is_numeric_type ()
623  && (! pivoted || argp.is_perm_matrix ())))
624  error ("luupdate: L, U, X, and Y must be numeric");
625 
626  if (! check_lu_dims (argl, argu, argp))
627  error ("luupdate: dimension mismatch");
628 
629  PermMatrix P = (pivoted
630  ? argp.perm_matrix_value ()
631  : PermMatrix::eye (argl.rows ()));
632 
633  if (argl.is_real_type () && argu.is_real_type ()
634  && argx.is_real_type () && argy.is_real_type ())
635  {
636  // all real case
637  if (argl.is_single_type () || argu.is_single_type ()
638  || argx.is_single_type () || argy.is_single_type ())
639  {
640  FloatMatrix L = argl.float_matrix_value ();
641  FloatMatrix U = argu.float_matrix_value ();
642  FloatMatrix x = argx.float_matrix_value ();
643  FloatMatrix y = argy.float_matrix_value ();
644 
645  octave::math::lu<FloatMatrix> fact (L, U, P);
646  if (pivoted)
647  fact.update_piv (x, y);
648  else
649  fact.update (x, y);
650 
651  if (pivoted)
652  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
653  else
654  return ovl (get_lu_l (fact), get_lu_u (fact));
655  }
656  else
657  {
658  Matrix L = argl.matrix_value ();
659  Matrix U = argu.matrix_value ();
660  Matrix x = argx.matrix_value ();
661  Matrix y = argy.matrix_value ();
662 
663  octave::math::lu<Matrix> fact (L, U, P);
664  if (pivoted)
665  fact.update_piv (x, y);
666  else
667  fact.update (x, y);
668 
669  if (pivoted)
670  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
671  else
672  return ovl (get_lu_l (fact), get_lu_u (fact));
673  }
674  }
675  else
676  {
677  // complex case
678  if (argl.is_single_type () || argu.is_single_type ()
679  || argx.is_single_type () || argy.is_single_type ())
680  {
685 
687  if (pivoted)
688  fact.update_piv (x, y);
689  else
690  fact.update (x, y);
691 
692  if (pivoted)
693  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
694  else
695  return ovl (get_lu_l (fact), get_lu_u (fact));
696  }
697  else
698  {
703 
704  octave::math::lu<ComplexMatrix> fact (L, U, P);
705  if (pivoted)
706  fact.update_piv (x, y);
707  else
708  fact.update (x, y);
709 
710  if (pivoted)
711  return ovl (get_lu_l (fact), get_lu_u (fact), fact.P ());
712  else
713  return ovl (get_lu_l (fact), get_lu_u (fact));
714  }
715  }
716 }
717 
718 /*
719 %!shared A, u, v, Ac, uc, vc
720 %! A = [0.091364 0.613038 0.999083;
721 %! 0.594638 0.425302 0.603537;
722 %! 0.383594 0.291238 0.085574;
723 %! 0.265712 0.268003 0.238409;
724 %! 0.669966 0.743851 0.445057 ];
725 %!
726 %! u = [0.85082;
727 %! 0.76426;
728 %! 0.42883;
729 %! 0.53010;
730 %! 0.80683 ];
731 %!
732 %! v = [0.98810;
733 %! 0.24295;
734 %! 0.43167 ];
735 %!
736 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
737 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
738 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
739 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
740 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
741 %!
742 %! uc = [0.20351 + 0.05401i;
743 %! 0.13141 + 0.43708i;
744 %! 0.29808 + 0.08789i;
745 %! 0.69821 + 0.38844i;
746 %! 0.74871 + 0.25821i ];
747 %!
748 %! vc = [0.85839 + 0.29468i;
749 %! 0.20820 + 0.93090i;
750 %! 0.86184 + 0.34689i ];
751 %!
752 
753 %!testif HAVE_QRUPDATE_LUU
754 %! [L,U,P] = lu (A);
755 %! [L,U] = luupdate (L,U,P*u,v);
756 %! assert (norm (vec (tril (L)-L), Inf) == 0);
757 %! assert (norm (vec (triu (U)-U), Inf) == 0);
758 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
759 %!
760 %!testif HAVE_QRUPDATE_LUU
761 %! [L,U,P] = lu (Ac);
762 %! [L,U] = luupdate (L,U,P*uc,vc);
763 %! assert (norm (vec (tril (L)-L), Inf) == 0);
764 %! assert (norm (vec (triu (U)-U), Inf) == 0);
765 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
766 
767 %!testif HAVE_QRUPDATE_LUU
768 %! [L,U,P] = lu (single (A));
769 %! [L,U] = luupdate (L,U,P*single (u), single (v));
770 %! assert (norm (vec (tril (L)-L), Inf) == 0);
771 %! assert (norm (vec (triu (U)-U), Inf) == 0);
772 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
773 %!
774 %!testif HAVE_QRUPDATE_LUU
775 %! [L,U,P] = lu (single (Ac));
776 %! [L,U] = luupdate (L,U,P*single (uc),single (vc));
777 %! assert (norm (vec (tril (L)-L), Inf) == 0);
778 %! assert (norm (vec (triu (U)-U), Inf) == 0);
779 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
780 
781 %!testif HAVE_QRUPDATE_LUU
782 %! [L,U,P] = lu (A);
783 %! [L,U,P] = luupdate (L,U,P,u,v);
784 %! assert (norm (vec (tril (L)-L), Inf) == 0);
785 %! assert (norm (vec (triu (U)-U), Inf) == 0);
786 %! assert (norm (vec (P'*L*U - A - u*v.'), Inf) < norm (A)*1e1*eps);
787 %!
788 %!testif HAVE_QRUPDATE_LUU
789 %! [L,U,P] = lu (A);
790 %! [~,ordcols] = max (P,[],1);
791 %! [~,ordrows] = max (P,[],2);
792 %! P1 = eye (size(P))(:,ordcols);
793 %! P2 = eye (size(P))(ordrows,:);
794 %! assert(P1 == P);
795 %! assert(P2 == P);
796 %! [L,U,P] = luupdate (L,U,P,u,v);
797 %! [L,U,P1] = luupdate (L,U,P1,u,v);
798 %! [L,U,P2] = luupdate (L,U,P2,u,v);
799 %! assert(P1 == P);
800 %! assert(P2 == P);
801 %!
802 %!testif HAVE_QRUPDATE_LUU
803 %! [L,U,P] = lu (Ac);
804 %! [L,U,P] = luupdate (L,U,P,uc,vc);
805 %! assert (norm (vec (tril (L)-L), Inf) == 0);
806 %! assert (norm (vec (triu (U)-U), Inf) == 0);
807 %! assert (norm (vec (P'*L*U - Ac - uc*vc.'), Inf) < norm (Ac)*1e1*eps);
808 
809 %!testif HAVE_QRUPDATE_LUU
810 %! [L,U,P] = lu (single (A));
811 %! [L,U,P] = luupdate (L,U,P,single (u),single (v));
812 %! assert (norm (vec (tril (L)-L), Inf) == 0);
813 %! assert (norm (vec (triu (U)-U), Inf) == 0);
814 %! assert (norm (vec (P'*L*U - single (A) - single (u)*single (v).'), Inf) < norm (single (A))*1e1*eps ("single"));
815 %!
816 %!testif HAVE_QRUPDATE_LUU
817 %! [L,U,P] = lu (single (Ac));
818 %! [L,U,P] = luupdate (L,U,P,single (uc),single (vc));
819 %! assert (norm (vec (tril (L)-L), Inf) == 0);
820 %! assert (norm (vec (triu (U)-U), Inf) == 0);
821 %! assert (norm (vec (P'*L*U - single (Ac) - single (uc)*single (vc).'), Inf) < norm (single (Ac))*1e1*eps ("single"));
822 */
ColumnVector Pc_vec(void) const
Definition: sparse-lu.cc:938
bool is_real_type(void) const
Definition: ov.h:667
static octave_value get_lu_l(const octave::math::lu< MT > &fact)
Definition: lu.cc:39
int ndims(void) const
Definition: ov.h:495
octave_idx_type rows(void) const
Definition: ov.h:489
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:145
ar P
Definition: __luinc__.cc:49
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:809
OCTAVE_EXPORT octave_value_list isa nd deftypefn *return ovl(args(0).is_integer_type())
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
PermMatrix Pr_mat(void) const
Definition: sparse-lu.cc:911
lu_type L(void) const
Definition: sparse-lu.h:82
bool is_numeric_type(void) const
Definition: ov.h:679
for large enough k
Definition: lu.cc:606
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
PermMatrix transpose(void) const
Definition: PermMatrix.cc:110
bool is_perm_matrix(void) const
Definition: ov.h:575
ColumnVector Pr_vec(void) const
Definition: sparse-lu.cc:897
u
Definition: lu.cc:138
const octave_idx_type * col_perm(void) const
Definition: sparse-lu.h:104
octave_value arg
Definition: pr-output.cc:3440
lu_type Y(void) const
Definition: sparse-lu.cc:840
T L(void) const
Definition: lu.cc:76
JNIEnv void * args
Definition: ov-java.cc:67
bool regular(void) const
Definition: lu.cc:192
octave_idx_type columns(void) const
Definition: ov.h:491
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:935
const octave_idx_type * row_perm(void) const
Definition: sparse-lu.h:102
SparseMatrix R(void) const
Definition: sparse-lu.h:86
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
int nargin
Definition: graphics.cc:10115
bool is_complex_type(void) const
Definition: ov.h:670
double tmp
Definition: data.cc:6300
PermMatrix P(void) const
Definition: lu.cc:169
octave_value retval
Definition: data.cc:6294
void update_piv(const VT &u, const VT &v)
Definition: dMatrix.h:37
T U(void) const
Definition: lu.cc:103
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:838
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:787
void err_wrong_type_arg(const char *name, const char *s)
Definition: errwarn.cc:156
bool is_empty(void) const
Definition: ov.h:542
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:805
p
Definition: lu.cc:138
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:790
PermMatrix Pc_mat(void) const
Definition: sparse-lu.cc:952
the element is set to zero In other the statement xample y
Definition: data.cc:5342
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5150
bool is_undefined(void) const
Definition: ov.h:539
void update(const VT &u, const VT &v)
void resize(octave_idx_type n, const octave_value &rfv=octave_value())
Definition: ovl.h:100
ColumnVector P_vec(void) const
Definition: lu.cc:176
PermMatrix perm_matrix_value(void) const
Definition: ov.h:857
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:834
bool is_single_type(void) const
Definition: ov.h:627
static PermMatrix eye(octave_idx_type n)
Definition: PermMatrix.cc:214
lu_type U(void) const
Definition: sparse-lu.h:84
T Y(void) const
Definition: lu.cc:127
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:854
Definition: mx-defs.h:75
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &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 F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205