GNU Octave  4.2.1 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qr.cc
Go to the documentation of this file.
1 /*
2
3 Copyright (C) 1996-2017 John W. Eaton
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 Copyright (C) 2008-2009 VZLU Prague
6
7 This file is part of Octave.
8
9 Octave is free software; you can redistribute it and/or modify it
11 Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 Octave is distributed in the hope that it will be useful, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
17 for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with Octave; see the file COPYING. If not, see
22
23 */
24
25 #if defined (HAVE_CONFIG_H)
26 # include "config.h"
27 #endif
28
29 #include "qr.h"
30 #include "qrp.h"
31 #include "sparse-qr.h"
32
33 #include "defun-dld.h"
34 #include "error.h"
35 #include "errwarn.h"
36 #include "ovl.h"
37
38 template <typename MT>
39 static octave_value
41 {
42  MT R = fact.R ();
43  if (R.is_square () && fact.regular ())
45  else
46  return R;
47 }
48
49 template <typename T>
50 static typename octave::math::qr<T>::type
52 {
53  return ((nargout == 0 || nargout == 1)
56 }
57
58 // [Q, R] = qr (X): form Q unitary and R upper triangular such
59 // that Q * R = X
60 //
61 // [Q, R] = qr (X, 0): form the economy decomposition such that if X is
62 // m by n then only the first n columns of Q are
63 // computed.
64 //
65 // [Q, R, P] = qr (X): form QRP factorization of X where
66 // P is a permutation matrix such that
67 // A * P = Q * R
68 //
69 // [Q, R, P] = qr (X, 0): form the economy decomposition with
70 // permutation vector P such that Q * R = X (:, P)
71 //
72 // qr (X) alone returns the output of the LAPACK routine dgeqrf, such
73 // that R = triu (qr (X))
74
76  doc: /* -*- texinfo -*-
77 @deftypefn {} {[@var{Q}, @var{R}] =} qr (@var{A})
78 @deftypefnx {} {[@var{Q}, @var{R}, @var{P}] =} qr (@var{A}) # non-sparse A
79 @deftypefnx {} {@var{X} =} qr (@var{A})
80 @deftypefnx {} {@var{R} =} qr (@var{A}) # sparse A
81 @deftypefnx {} {[@var{C}, @var{R}] =} qr (@var{A}, @var{B})
82 @deftypefnx {} {[@dots{}] =} qr (@dots{}, 0)
83 @deftypefnx {} {[@dots{}] =} qr (@dots{}, 'vector')
84 @deftypefnx {} {[@dots{}] =} qr (@dots{}, 'matrix')
85 @cindex QR factorization
86 Compute the QR@tie{}factorization of @var{A}, using standard @sc{lapack}
87 subroutines.
88 The QR@tie{}factorization is
89 @tex
90 $QR = A$ where $Q$ is an orthogonal matrix and $R$ is upper triangular.
91 @end tex
92 @ifnottex
93 @code{@var{Q} * @var{R} = @var{A}} where @var{Q} is an orthogonal matrix and
94 @var{R} is upper triangular.
95 @end ifnottex
96
97 For example, given the matrix @code{@var{A} = [1, 2; 3, 4]},
98
99 @example
100 [@var{Q}, @var{R}] = qr (@var{A})
101 @end example
102
103 @noindent
104 returns
105
106 @example
107 @group
108 @var{Q} =
109
110  -0.31623 -0.94868
111  -0.94868 0.31623
112
113 @var{R} =
114
115  -3.16228 -4.42719
116  0.00000 -0.63246
117 @end group
118 @end example
119
120 The @code{qr} factorization has applications in the solution of least
121 squares problems
122 @tex
123 $$124 \min_x \left\Vert A x - b \right\Vert_2 125$$
126 @end tex
127 @ifnottex
128
129 @example
130 min norm(A x - b)
131 @end example
132
133 @end ifnottex
134 for overdetermined systems of equations (i.e.,
135 @tex
136 $A$
137 @end tex
138 @ifnottex
139 @var{A}
140 @end ifnottex
141 is a tall, thin matrix).
142
143 If only a single return value is requested, then it is either @var{R} if
144 @var{A} is sparse, or @var{X} such that @code{@var{R} = triu (@var{X})} if
145 @var{A} is full.
146 (Note: Unlike most commands, the single return value is not
147 the first return value when multiple are requested.)
148
149
150 If the matrix @var{A} is full, the permuted QR@tie{}factorization
151 @code{[@var{Q}, @var{R}, @var{P}] = qr (@var{A})} forms the
152 QR@tie{}factorization such that the diagonal entries of @var{R} are
153 decreasing in magnitude order. For example, given the matrix
154 @code{a = [1, 2; 3, 4]},
155
156 @example
157 [@var{Q}, @var{R}, @var{P}] = qr (@var{A})
158 @end example
159
160 @noindent
161 returns
162
163 @example
164 @group
165 @var{Q} =
166
167  -0.44721 -0.89443
168  -0.89443 0.44721
169
170 @var{R} =
171
172  -4.47214 -3.13050
173  0.00000 0.44721
174
175 @var{P} =
176
177  0 1
178  1 0
179 @end group
180 @end example
181
182 The permuted @code{qr} factorization
183 @code{[@var{Q}, @var{R}, @var{P}] = qr (@var{A})} factorization allows the
184 construction of an orthogonal basis of @code{span (A)}.
185
186 If the matrix @var{A} is sparse, then the sparse QR@tie{}factorization
187 of @var{A} is computed using @sc{CSparse}.
188 As the matrix @var{Q} is in general a full matrix,
189 it is recommended to request only one return value,
190 which is the @var{Q}-less factorization @var{R} of @var{A}, such that
191 @code{@var{R} = chol (@var{A}' * @var{A})}.
192
193 If an additional matrix @var{B} is supplied and two return values are
194 requested, then @code{qr} returns
195 @var{C}, where @code{@var{C} = @var{Q}' * @var{B}}. This allows the
196 least squares approximation of @code{@var{A} \ @var{B}} to be calculated
197 as
198
199 @example
200 @group
201 [@var{C}, @var{R}] = qr (@var{A}, @var{B})
202 x = @var{R} \ @var{C}
203 @end group
204 @end example
205
206 If the final argument is the scalar 0 and the number of rows is larger
207 than the number of columns, then an @qcode{"economy"} factorization is
208 returned, omitting zeroes of @var{R} and the corresponding columns of @var{Q}.
209 That is, @var{R} will have only @code{size (@var{A},1)} rows. In this case,
210 @var{P} is a vector rather than a matrix.
211
212 If the final argument is the string @qcode{"vector"} then @var{P} is a
213 permutation vector instead of a permutation matrix.
214
215 @seealso{chol, hess, lu, qz, schur, svd, qrupdate, qrinsert, qrdelete, qrshift}
216 @end deftypefn */)
217 {
218  int nargin = args.length ();
219
220  if (nargin < 1 || nargin > 3)
221  print_usage ();
222
224
225  octave_value arg = args(0);
226
227  bool economy = false;
228  bool is_cmplx = false;
229  bool have_b = 0;
230  bool vector_p = 0;
231
232  if (arg.is_complex_type ())
233  is_cmplx = true;
234  if (nargin > 1)
235  {
236  have_b = true;
237  if (args(nargin-1).is_scalar_type ())
238  {
239  int val = args(nargin-1).int_value ();
240  if (val == 0)
241  {
242  economy = true;
243  have_b = (nargin > 2);
244  }
245  else if (nargin == 3) // argument 3 should be 0 or a string
246  print_usage ();
247  }
248  else if (args(nargin-1).is_string ())
249  {
250  if (args(nargin-1).string_value () == "vector")
251  vector_p = true;
252  else if (args(nargin-1).string_value () != "matrix")
253  error ("qr: type for P must be 'matrix' or 'vector', not %s",
254  args(nargin-1).string_value ().c_str ());
255  have_b = (nargin > 2);
256  }
257  else if (! args(nargin-1).is_matrix_type ())
258  err_wrong_type_arg ("qr", args(nargin-1));
259  else if (nargin == 3) // should be caught by is_scalar_type or is_string
260  print_usage ();
261
262  if (have_b && args(1).is_complex_type ())
263  is_cmplx = true;
264  }
265
266  if (arg.is_sparse_type ())
267  {
268  if (nargout > 2)
269  error ("qr: Permutation output is not supported for sparse input");
270
271  if (is_cmplx)
272  {
274
275  if (have_b)
276  {
277  retval = ovl (q.C (args(1).complex_matrix_value ()),
278  q.R (economy));
279  if (arg.rows () < arg.columns ())
280  warning ("qr: non minimum norm solution for under-determined "
281  "problem %dx%d", arg.rows (), arg.columns ());
282  }
283  else if (nargout > 1)
284  retval = ovl (q.Q (), q.R (economy));
285  else
286  retval = ovl (q.R (economy));
287  }
288  else
289  {
291
292  if (have_b)
293  {
294  retval = ovl (q.C (args(1).matrix_value ()), q.R (economy));
295  if (arg.rows () < arg.columns ())
296  warning ("qr: non minimum norm solution for under-determined "
297  "problem %dx%d", arg.rows (), arg.columns ());
298  }
299  else if (nargout > 1)
300  retval = ovl (q.Q (), q.R (economy));
301  else
302  retval = ovl (q.R (economy));
303  }
304  }
305  else
306  {
307  if (arg.is_single_type ())
308  {
309  if (arg.is_real_type ())
310  {
312  = qr_type<FloatMatrix> (nargin, nargout);
313
315
316  switch (nargout)
317  {
318  case 0:
319  case 1:
320  {
321  octave::math::qr<FloatMatrix> fact (m, type);
322  retval = ovl (fact.R ());
323  }
324  break;
325
326  case 2:
327  {
328  octave::math::qr<FloatMatrix> fact (m, type);
329  retval = ovl (fact.Q (), get_qr_r (fact));
330  if (have_b)
331  {
332  if (is_cmplx)
333  retval(0) = fact.Q ().transpose ()
334  * args(1).float_complex_matrix_value ();
335  else
336  retval(0) = fact.Q ().transpose ()
337  * args(1).float_matrix_value ();
338  }
339  }
340  break;
341
342  default:
343  {
344  octave::math::qrp<FloatMatrix> fact (m, type);
345
346  if (type == octave::math::qr<FloatMatrix>::economy || vector_p)
347  retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
348  else
349  retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
350  }
351  break;
352  }
353  }
354  else if (arg.is_complex_type ())
355  {
357  = qr_type<FloatComplexMatrix> (nargin, nargout);
358
360
361  switch (nargout)
362  {
363  case 0:
364  case 1:
365  {
367  retval = ovl (fact.R ());
368  }
369  break;
370
371  case 2:
372  {
374  retval = ovl (fact.Q (), get_qr_r (fact));
375  if (have_b)
376  retval (0) = conj (fact.Q ().transpose ())
377  * args(1).float_complex_matrix_value ();
378  }
379  break;
380
381  default:
382  {
384  if (type == octave::math::qr<FloatComplexMatrix>::economy || vector_p)
385  retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
386  else
387  retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
388  }
389  break;
390  }
391  }
392  }
393  else
394  {
395  if (arg.is_real_type ())
396  {
397  octave::math::qr<Matrix>::type type = qr_type<Matrix> (nargin, nargout);
398
399  Matrix m = arg.matrix_value ();
400
401  switch (nargout)
402  {
403  case 0:
404  case 1:
405  {
406  octave::math::qr<Matrix> fact (m, type);
407  retval = ovl (fact.R ());
408  }
409  break;
410
411  case 2:
412  {
413  octave::math::qr<Matrix> fact (m, type);
414  retval = ovl (fact.Q (), get_qr_r (fact));
415  if (have_b)
416  {
417  if (is_cmplx)
418  retval(0) = fact.Q ().transpose ()
419  * args(1).complex_matrix_value ();
420  else
421  retval(0) = fact.Q ().transpose ()
422  * args(1).matrix_value ();
423  }
424  }
425  break;
426
427  default:
428  {
429  octave::math::qrp<Matrix> fact (m, type);
430  if (type == octave::math::qr<Matrix>::economy || vector_p)
431  retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
432  else
433  retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
434  }
435  break;
436  }
437  }
438  else if (arg.is_complex_type ())
439  {
441  = qr_type<ComplexMatrix> (nargin, nargout);
442
444
445  switch (nargout)
446  {
447  case 0:
448  case 1:
449  {
450  octave::math::qr<ComplexMatrix> fact (m, type);
451  retval = ovl (fact.R ());
452  }
453  break;
454
455  case 2:
456  {
457  octave::math::qr<ComplexMatrix> fact (m, type);
458  retval = ovl (fact.Q (), get_qr_r (fact));
459  if (have_b)
460  retval (0) = conj (fact.Q ().transpose ())
461  * args(1).complex_matrix_value ();
462  }
463  break;
464
465  default:
466  {
467  octave::math::qrp<ComplexMatrix> fact (m, type);
468  if (type == octave::math::qr<ComplexMatrix>::economy || vector_p)
469  retval = ovl (fact.Q (), get_qr_r (fact), fact.Pvec ());
470  else
471  retval = ovl (fact.Q (), get_qr_r (fact), fact.P ());
472  }
473  break;
474  }
475  }
476  else
477  err_wrong_type_arg ("qr", arg);
478  }
479  }
480
481  return retval;
482 }
483
484 /*
485 %!test
486 %! a = [0, 2, 1; 2, 1, 2];
487 %!
488 %! [q, r] = qr (a);
489 %! [qe, re] = qr (a, 0);
490 %!
491 %! assert (q * r, a, sqrt (eps));
492 %! assert (qe * re, a, sqrt (eps));
493
494 %!test
495 %! a = [0, 2, 1; 2, 1, 2];
496 %!
497 %! [q, r] = qr (a);
498 %! [qe, re] = qr (a, 0);
499 %!
500 %! assert (q * r, a, sqrt (eps));
501 %! assert (qe * re, a, sqrt (eps));
502
503 %!test
504 %! a = [0, 2, 1; 2, 1, 2];
505 %!
506 %! [q, r, p] = qr (a); # FIXME: not giving right dimensions.
507 %! [qe, re, pe] = qr (a, 0);
508 %!
509 %! assert (q * r, a * p, sqrt (eps));
510 %! assert (qe * re, a(:, pe), sqrt (eps));
511
512 %!test
513 %! a = [0, 2; 2, 1; 1, 2];
514 %!
515 %! [q, r] = qr (a);
516 %! [qe, re] = qr (a, 0);
517 %!
518 %! assert (q * r, a, sqrt (eps));
519 %! assert (qe * re, a, sqrt (eps));
520
521 %!test
522 %! a = [0, 2; 2, 1; 1, 2];
523 %!
524 %! [q, r, p] = qr (a);
525 %! [qe, re, pe] = qr (a, 0);
526 %!
527 %! assert (q * r, a * p, sqrt (eps));
528 %! assert (qe * re, a(:, pe), sqrt (eps));
529
530 %!test
531 %! a = [0, 2, 1; 2, 1, 2; 3, 1, 2];
532 %! b = [1, 3, 2; 1, 1, 0; 3, 0, 2];
533 %!
534 %! [q, r] = qr (a);
535 %! [c, re] = qr (a, b);
536 %!
537 %! assert (r, re, sqrt (eps));
538 %! assert (q'*b, c, sqrt (eps));
539
540 %!test
541 %! a = [0, 2, i; 2, 1, 2; 3, 1, 2];
542 %! b = [1, 3, 2; 1, i, 0; 3, 0, 2];
543 %!
544 %! [q, r] = qr (a);
545 %! [c, re] = qr (a, b);
546 %!
547 %! assert (r, re, sqrt (eps));
548 %! assert (q'*b, c, sqrt (eps));
549
550 %!test
551 %! a = [0, 2, i; 2, 1, 2; 3, 1, 2];
552 %! b = [1, 3, 2; 1, 1, 0; 3, 0, 2];
553 %!
554 %! [q, r] = qr (a);
555 %! [c, re] = qr (a, b);
556 %!
557 %! assert (r, re, sqrt (eps));
558 %! assert (q'*b, c, sqrt (eps));
559
560 %!test
561 %! a = [0, 2, 1; 2, 1, 2; 3, 1, 2];
562 %! b = [1, 3, 2; 1, i, 0; 3, 0, 2];
563 %!
564 %! [q, r] = qr (a);
565 %! [c, re] = qr (a, b);
566 %!
567 %! assert (r, re, sqrt (eps));
568 %! assert (q'*b, c, sqrt (eps));
569
570 %!test
571 %! assert (qr (zeros (0, 0)), zeros (0, 0))
572 %! assert (qr (zeros (1, 0)), zeros (1, 0))
573 %! assert (qr (zeros (0, 1)), zeros (0, 1))
574
575 %!error qr ()
576 %!error qr ([1, 2; 3, 4], 0, 2)
577
578 %!function retval = __testqr (q, r, a, p)
579 %! tol = 100*eps (class (q));
580 %! retval = 0;
581 %! if (nargin == 3)
582 %! n1 = norm (q*r - a);
583 %! n2 = norm (q'*q - eye (columns (q)));
584 %! retval = (n1 < tol && n2 < tol);
585 %! else
586 %! n1 = norm (q'*q - eye (columns (q)));
587 %! retval = (n1 < tol);
588 %! if (isvector (p))
589 %! n2 = norm (q*r - a(:,p));
590 %! retval = (retval && n2 < tol);
591 %! else
592 %! n2 = norm (q*r - a*p);
593 %! retval = (retval && n2 < tol);
594 %! endif
595 %! endif
596 %!endfunction
597
598 %!test
599 %! t = ones (24, 1);
600 %! j = 1;
601 %!
602 %! if (false) # eliminate big matrix tests
603 %! a = rand (5000, 20);
604 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
605 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
606 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
607 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
608 %!
609 %! a = a+1i*eps;
610 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
611 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
612 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
613 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
614 %! endif
615 %!
616 %! a = [ ones(1,15); sqrt(eps)*eye(15) ];
617 %! [q,r] = qr (a); t(j++) = __testqr (q, r, a);
618 %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a');
619 %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p);
620 %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p);
621 %!
622 %! a = a+1i*eps;
623 %! [q,r] = qr (a); t(j++) = __testqr (q, r, a);
624 %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a');
625 %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p);
626 %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p);
627 %!
628 %! a = [ ones(1,15); sqrt(eps)*eye(15) ];
629 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
630 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
631 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
632 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
633 %!
634 %! a = a+1i*eps;
635 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
636 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
637 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
638 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
639 %!
640 %! a = [ 611 196 -192 407 -8 -52 -49 29
641 %! 196 899 113 -192 -71 -43 -8 -44
642 %! -192 113 899 196 61 49 8 52
643 %! 407 -192 196 611 8 44 59 -23
644 %! -8 -71 61 8 411 -599 208 208
645 %! -52 -43 49 44 -599 411 208 208
646 %! -49 -8 8 59 208 208 99 -911
647 %! 29 -44 52 -23 208 208 -911 99 ];
648 %! [q,r] = qr (a);
649 %!
650 %! assert (all (t) && norm (q*r - a) < 5000*eps);
651
652 %!test
653 %! a = single ([0, 2, 1; 2, 1, 2]);
654 %!
655 %! [q, r] = qr (a);
656 %! [qe, re] = qr (a, 0);
657 %!
658 %! assert (q * r, a, sqrt (eps ("single")));
659 %! assert (qe * re, a, sqrt (eps ("single")));
660
661 %!test
662 %! a = single ([0, 2, 1; 2, 1, 2]);
663 %!
664 %! [q, r, p] = qr (a); # FIXME: not giving right dimensions.
665 %! [qe, re, pe] = qr (a, 0);
666 %!
667 %! assert (q * r, a * p, sqrt (eps ("single")));
668 %! assert (qe * re, a(:, pe), sqrt (eps ("single")));
669
670 %!test
671 %! a = single ([0, 2; 2, 1; 1, 2]);
672 %!
673 %! [q, r] = qr (a);
674 %! [qe, re] = qr (a, 0);
675 %!
676 %! assert (q * r, a, sqrt (eps ("single")));
677 %! assert (qe * re, a, sqrt (eps ("single")));
678
679 %!test
680 %! a = single ([0, 2; 2, 1; 1, 2]);
681 %!
682 %! [q, r, p] = qr (a);
683 %! [qe, re, pe] = qr (a, 0);
684 %!
685 %! assert (q * r, a * p, sqrt (eps ("single")));
686 %! assert (qe * re, a(:, pe), sqrt (eps ("single")));
687
688 %!test
689 %! a = single([0, 2, 1; 2, 1, 2; 3, 1, 2]);
690 %! b = single([1, 3, 2; 1, 1, 0; 3, 0, 2]);
691 %!
692 %! [q, r] = qr (a);
693 %! [c, re] = qr (a, b);
694 %!
695 %! assert (r, re, sqrt (eps ("single")));
696 %! assert (q'*b, c, sqrt (eps ("single")));
697
698 %!test
699 %! a = single([0, 2, i; 2, 1, 2; 3, 1, 2]);
700 %! b = single([1, 3, 2; 1, i, 0; 3, 0, 2]);
701 %!
702 %! [q, r] = qr (a);
703 %! [c, re] = qr (a, b);
704 %!
705 %! assert (r, re, sqrt (eps ("single")));
706 %! assert (q'*b, c, sqrt (eps ("single")));
707
708 %!test
709 %! a = single([0, 2, i; 2, 1, 2; 3, 1, 2]);
710 %! b = single([1, 3, 2; 1, 1, 0; 3, 0, 2]);
711 %!
712 %! [q, r] = qr (a);
713 %! [c, re] = qr (a, b);
714 %!
715 %! assert (r, re, sqrt (eps));
716 %! assert (q'*b, c, sqrt (eps));
717
718 %!test
719 %! a = single([0, 2, 1; 2, 1, 2; 3, 1, 2]);
720 %! b = single([1, 3, 2; 1, i, 0; 3, 0, 2]);
721 %!
722 %! [q, r] = qr (a);
723 %! [c, re] = qr (a, b);
724 %!
725 %! assert (r, re, sqrt (eps ("single")));
726 %! assert (q'*b, c, sqrt (eps ("single")));
727
728 %!error qr ()
729 %!error qr ([1, 2; 3, 4], 0, 2)
730
731 %!test
732 %! t = ones (24, 1);
733 %! j = 1;
734 %!
735 %! if (false) # eliminate big matrix tests
736 %! a = rand (5000,20);
737 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
738 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
739 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
740 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
741 %!
742 %! a = a+1i*eps ("single");
743 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
744 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
745 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
746 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
747 %! endif
748 %!
749 %! a = [ ones(1,15); sqrt(eps("single"))*eye(15) ];
750 %! [q,r] = qr (a); t(j++) = __testqr (q, r, a);
751 %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a');
752 %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p);
753 %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p);
754 %!
755 %! a = a+1i*eps ("single");
756 %! [q,r] = qr (a); t(j++) = __testqr (q, r, a);
757 %! [q,r] = qr (a'); t(j++) = __testqr (q, r, a');
758 %! [q,r,p] = qr (a); t(j++) = __testqr (q, r, a, p);
759 %! [q,r,p] = qr (a'); t(j++) = __testqr (q, r, a', p);
760 %!
761 %! a = [ ones(1,15); sqrt(eps("single"))*eye(15) ];
762 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
763 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
764 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
765 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a', p);
766 %!
767 %! a = a+1i*eps ("single");
768 %! [q,r] = qr (a, 0); t(j++) = __testqr (q, r, a);
769 %! [q,r] = qr (a',0); t(j++) = __testqr (q, r, a');
770 %! [q,r,p] = qr (a, 0); t(j++) = __testqr (q, r, a, p);
771 %! [q,r,p] = qr (a',0); t(j++) = __testqr (q, r, a',p);
772 %!
773 %! a = [ 611 196 -192 407 -8 -52 -49 29
774 %! 196 899 113 -192 -71 -43 -8 -44
775 %! -192 113 899 196 61 49 8 52
776 %! 407 -192 196 611 8 44 59 -23
777 %! -8 -71 61 8 411 -599 208 208
778 %! -52 -43 49 44 -599 411 208 208
779 %! -49 -8 8 59 208 208 99 -911
780 %! 29 -44 52 -23 208 208 -911 99 ];
781 %! [q,r] = qr (a);
782 %!
783 %! assert (all (t) && norm (q*r-a) < 5000*eps ("single"));
784
785 ## The deactivated tests below can't be tested till rectangular back-subs is
786 ## implemented for sparse matrices.
787
788 %!testif HAVE_CXSPARSE
789 %! n = 20; d = 0.2;
790 %! a = sprandn (n,n,d) + speye (n,n);
791 %! r = qr (a);
792 %! assert (r'*r, a'*a, 1e-10);
793
794 %!testif HAVE_COLAMD
795 %! n = 20; d = 0.2;
796 %! a = sprandn (n,n,d) + speye (n,n);
797 %! q = symamd (a);
798 %! a = a(q,q);
799 %! r = qr (a);
800 %! assert (r'*r, a'*a, 1e-10);
801
802 %!testif HAVE_CXSPARSE
803 %! n = 20; d = 0.2;
804 %! a = sprandn (n,n,d) + speye (n,n);
805 %! [c,r] = qr (a, ones (n,1));
806 %! assert (r\c, full (a)\ones (n,1), 10e-10);
807
808 %!testif HAVE_CXSPARSE
809 %! n = 20; d = 0.2;
810 %! a = sprandn (n,n,d) + speye (n,n);
811 %! b = randn (n,2);
812 %! [c,r] = qr (a, b);
813 %! assert (r\c, full (a)\b, 10e-10);
814
815 %% Test under-determined systems!!
816 %!#testif HAVE_CXSPARSE
817 %! n = 20; d = 0.2;
818 %! a = sprandn (n,n+1,d) + speye (n,n+1);
819 %! b = randn (n,2);
820 %! [c,r] = qr (a, b);
821 %! assert (r\c, full (a)\b, 10e-10);
822
823 %!testif HAVE_CXSPARSE
824 %! n = 20; d = 0.2;
825 %! a = 1i*sprandn (n,n,d) + speye (n,n);
826 %! r = qr (a);
827 %! assert (r'*r,a'*a,1e-10);
828
829 %!testif HAVE_COLAMD
830 %! n = 20; d = 0.2;
831 %! a = 1i*sprandn (n,n,d) + speye (n,n);
832 %! q = symamd (a);
833 %! a = a(q,q);
834 %! r = qr (a);
835 %! assert (r'*r, a'*a, 1e-10);
836
837 %!testif HAVE_CXSPARSE
838 %! n = 20; d = 0.2;
839 %! a = 1i*sprandn (n,n,d) + speye (n,n);
840 %! [c,r] = qr (a, ones (n,1));
841 %! assert (r\c, full (a)\ones (n,1), 10e-10);
842
843 %!testif HAVE_CXSPARSE
844 %! n = 20; d = 0.2;
845 %! a = 1i*sprandn (n,n,d) + speye (n,n);
846 %! b = randn (n,2);
847 %! [c,r] = qr (a, b);
848 %! assert (r\c, full (a)\b, 10e-10);
849
850 %% Test under-determined systems!!
851 %!#testif HAVE_CXSPARSE
852 %! n = 20; d = 0.2;
853 %! a = 1i*sprandn (n,n+1,d) + speye (n,n+1);
854 %! b = randn (n,2);
855 %! [c,r] = qr (a, b);
856 %! assert (r\c, full (a)\b, 10e-10);
857
858 */
859
860 static
861 bool check_qr_dims (const octave_value& q, const octave_value& r,
862  bool allow_ecf = false)
863 {
864  octave_idx_type m = q.rows ();
865  octave_idx_type k = r.rows ();
866  octave_idx_type n = r.columns ();
867  return ((q.ndims () == 2 && r.ndims () == 2 && k == q.columns ())
868  && (m == k || (allow_ecf && k == n && k < m)));
869 }
870
871 static
872 bool check_index (const octave_value& i, bool vector_allowed = false)
873 {
874  return ((i.is_real_type () || i.is_integer_type ())
875  && (i.is_scalar_type () || vector_allowed));
876 }
877
878 DEFUN_DLD (qrupdate, args, ,
879  doc: /* -*- texinfo -*-
880 @deftypefn {} {[@var{Q1}, @var{R1}] =} qrupdate (@var{Q}, @var{R}, @var{u}, @var{v})
881 Given a QR@tie{}factorization of a real or complex matrix
882 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and
883 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of
884 @w{@var{A} + @var{u}*@var{v}'}, where @var{u} and @var{v} are column vectors
885 (rank-1 update) or matrices with equal number of columns
886 (rank-k update). Notice that the latter case is done as a sequence of
887 rank-1 updates; thus, for k large enough, it will be both faster and more
888 accurate to recompute the factorization from scratch.
889
890 The QR@tie{}factorization supplied may be either full (Q is square) or
891 economized (R is square).
892
893 @seealso{qr, qrinsert, qrdelete, qrshift}
894 @end deftypefn */)
895 {
897
898  if (args.length () != 4)
899  print_usage ();
900
901  octave_value argq = args(0);
902  octave_value argr = args(1);
903  octave_value argu = args(2);
904  octave_value argv = args(3);
905
906  if (! argq.is_numeric_type () || ! argr.is_numeric_type ()
907  || ! argu.is_numeric_type () || ! argv.is_numeric_type ())
908  print_usage ();
909
910  if (! check_qr_dims (argq, argr, true))
911  error ("qrupdate: Q and R dimensions don't match");
912
913  if (argq.is_real_type () && argr.is_real_type () && argu.is_real_type ()
914  && argv.is_real_type ())
915  {
916  // all real case
917  if (argq.is_single_type () || argr.is_single_type ()
918  || argu.is_single_type () || argv.is_single_type ())
919  {
920  FloatMatrix Q = argq.float_matrix_value ();
921  FloatMatrix R = argr.float_matrix_value ();
922  FloatMatrix u = argu.float_matrix_value ();
923  FloatMatrix v = argv.float_matrix_value ();
924
925  octave::math::qr<FloatMatrix> fact (Q, R);
926  fact.update (u, v);
927
928  retval = ovl (fact.Q (), get_qr_r (fact));
929  }
930  else
931  {
932  Matrix Q = argq.matrix_value ();
933  Matrix R = argr.matrix_value ();
934  Matrix u = argu.matrix_value ();
935  Matrix v = argv.matrix_value ();
936
937  octave::math::qr<Matrix> fact (Q, R);
938  fact.update (u, v);
939
940  retval = ovl (fact.Q (), get_qr_r (fact));
941  }
942  }
943  else
944  {
945  // complex case
946  if (argq.is_single_type () || argr.is_single_type ()
947  || argu.is_single_type () || argv.is_single_type ())
948  {
953
955  fact.update (u, v);
956
957  retval = ovl (fact.Q (), get_qr_r (fact));
958  }
959  else
960  {
965
967  fact.update (u, v);
968
969  retval = ovl (fact.Q (), get_qr_r (fact));
970  }
971  }
972
973  return retval;
974 }
975
976 /*
977 %!shared A, u, v, Ac, uc, vc
978 %! A = [0.091364 0.613038 0.999083;
979 %! 0.594638 0.425302 0.603537;
980 %! 0.383594 0.291238 0.085574;
981 %! 0.265712 0.268003 0.238409;
982 %! 0.669966 0.743851 0.445057 ];
983 %!
984 %! u = [0.85082;
985 %! 0.76426;
986 %! 0.42883;
987 %! 0.53010;
988 %! 0.80683 ];
989 %!
990 %! v = [0.98810;
991 %! 0.24295;
992 %! 0.43167 ];
993 %!
994 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
995 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
996 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
997 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
998 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
999 %!
1000 %! uc = [0.20351 + 0.05401i;
1001 %! 0.13141 + 0.43708i;
1002 %! 0.29808 + 0.08789i;
1003 %! 0.69821 + 0.38844i;
1004 %! 0.74871 + 0.25821i ];
1005 %!
1006 %! vc = [0.85839 + 0.29468i;
1007 %! 0.20820 + 0.93090i;
1008 %! 0.86184 + 0.34689i ];
1009 %!
1010
1011 %!test
1012 %! [Q,R] = qr (A);
1013 %! [Q,R] = qrupdate (Q, R, u, v);
1014 %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps);
1015 %! assert (norm (vec (triu (R)-R), Inf) == 0);
1016 %! assert (norm (vec (Q*R - A - u*v'), Inf) < norm (A)*1e1*eps);
1017 %!
1018 %!test
1019 %! [Q,R] = qr (Ac);
1020 %! [Q,R] = qrupdate (Q, R, uc, vc);
1021 %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps);
1022 %! assert (norm (vec (triu (R)-R), Inf) == 0);
1023 %! assert (norm (vec (Q*R - Ac - uc*vc'), Inf) < norm (Ac)*1e1*eps);
1024
1025 %!test
1026 %! [Q,R] = qr (single (A));
1027 %! [Q,R] = qrupdate (Q, R, single (u), single (v));
1028 %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single"));
1029 %! assert (norm (vec (triu (R)-R), Inf) == 0);
1030 %! assert (norm (vec (Q*R - single (A) - single (u)*single (v)'), Inf) < norm (single (A))*1e1*eps ("single"));
1031 %!
1032 %!test
1033 %! [Q,R] = qr (single (Ac));
1034 %! [Q,R] = qrupdate (Q, R, single (uc), single (vc));
1035 %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single"));
1036 %! assert (norm (vec (triu (R)-R), Inf) == 0);
1037 %! assert (norm (vec (Q*R - single (Ac) - single (uc)*single (vc)'), Inf) < norm (single (Ac))*1e1*eps ("single"));
1038 */
1039
1040 DEFUN_DLD (qrinsert, args, ,
1041  doc: /* -*- texinfo -*-
1042 @deftypefn {} {[@var{Q1}, @var{R1}] =} qrinsert (@var{Q}, @var{R}, @var{j}, @var{x}, @var{orient})
1043 Given a QR@tie{}factorization of a real or complex matrix
1044 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and
1045 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of
1046 @w{[A(:,1:j-1) x A(:,j:n)]}, where @var{u} is a column vector to be inserted
1047 into @var{A} (if @var{orient} is @qcode{"col"}), or the
1048 QR@tie{}factorization of @w{[A(1:j-1,:);x;A(:,j:n)]}, where @var{x} is a row
1049 vector to be inserted into @var{A} (if @var{orient} is @qcode{"row"}).
1050
1051 The default value of @var{orient} is @qcode{"col"}. If @var{orient} is
1052 @qcode{"col"}, @var{u} may be a matrix and @var{j} an index vector
1053 resulting in the QR@tie{}factorization of a matrix @var{B} such that
1054 @w{B(:,@var{j})} gives @var{u} and @w{B(:,@var{j}) = []} gives @var{A}.
1055 Notice that the latter case is done as a sequence of k insertions;
1056 thus, for k large enough, it will be both faster and more accurate to
1057 recompute the factorization from scratch.
1058
1059 If @var{orient} is @qcode{"col"}, the QR@tie{}factorization supplied may
1060 be either full (Q is square) or economized (R is square).
1061
1062 If @var{orient} is @qcode{"row"}, full factorization is needed.
1063 @seealso{qr, qrupdate, qrdelete, qrshift}
1064 @end deftypefn */)
1065 {
1066  int nargin = args.length ();
1067
1068  if (nargin < 4 || nargin > 5)
1069  print_usage ();
1070
1071  octave_value argq = args(0);
1072  octave_value argr = args(1);
1073  octave_value argj = args(2);
1074  octave_value argx = args(3);
1075
1076  if (! argq.is_numeric_type () || ! argr.is_numeric_type ()
1077  || ! argx.is_numeric_type ()
1078  || (nargin > 4 && ! args(4).is_string ()))
1079  print_usage ();
1080
1081  std::string orient = (nargin < 5) ? "col" : args(4).string_value ();
1082  bool col = (orient == "col");
1083
1084  if (! col && orient != "row")
1085  error ("qrinsert: ORIENT must be \"col\" or \"row\"");
1086
1087  if (! check_qr_dims (argq, argr, col) || (! col && argx.rows () != 1))
1088  error ("qrinsert: dimension mismatch");
1089
1090  if (! check_index (argj, col))
1091  error ("qrinsert: invalid index J");
1092
1094
1096
1097  octave_idx_type one = 1;
1098
1099  if (argq.is_real_type () && argr.is_real_type () && argx.is_real_type ())
1100  {
1101  // real case
1102  if (argq.is_single_type () || argr.is_single_type ()
1103  || argx.is_single_type ())
1104  {
1105  FloatMatrix Q = argq.float_matrix_value ();
1106  FloatMatrix R = argr.float_matrix_value ();
1107  FloatMatrix x = argx.float_matrix_value ();
1108
1109  octave::math::qr<FloatMatrix> fact (Q, R);
1110
1111  if (col)
1112  fact.insert_col (x, j-one);
1113  else
1114  fact.insert_row (x.row (0), j(0)-one);
1115
1116  retval = ovl (fact.Q (), get_qr_r (fact));
1117  }
1118  else
1119  {
1120  Matrix Q = argq.matrix_value ();
1121  Matrix R = argr.matrix_value ();
1122  Matrix x = argx.matrix_value ();
1123
1124  octave::math::qr<Matrix> fact (Q, R);
1125
1126  if (col)
1127  fact.insert_col (x, j-one);
1128  else
1129  fact.insert_row (x.row (0), j(0)-one);
1130
1131  retval = ovl (fact.Q (), get_qr_r (fact));
1132  }
1133  }
1134  else
1135  {
1136  // complex case
1137  if (argq.is_single_type () || argr.is_single_type ()
1138  || argx.is_single_type ())
1139  {
1140  FloatComplexMatrix Q =
1142  FloatComplexMatrix R =
1144  FloatComplexMatrix x =
1146
1148
1149  if (col)
1150  fact.insert_col (x, j-one);
1151  else
1152  fact.insert_row (x.row (0), j(0)-one);
1153
1154  retval = ovl (fact.Q (), get_qr_r (fact));
1155  }
1156  else
1157  {
1158  ComplexMatrix Q = argq.complex_matrix_value ();
1159  ComplexMatrix R = argr.complex_matrix_value ();
1160  ComplexMatrix x = argx.complex_matrix_value ();
1161
1162  octave::math::qr<ComplexMatrix> fact (Q, R);
1163
1164  if (col)
1165  fact.insert_col (x, j-one);
1166  else
1167  fact.insert_row (x.row (0), j(0)-one);
1168
1169  retval = ovl (fact.Q (), get_qr_r (fact));
1170  }
1171  }
1172
1173  return retval;
1174 }
1175
1176 /*
1177 %!test
1178 %! [Q,R] = qr (A);
1179 %! [Q,R] = qrinsert (Q, R, 3, u);
1180 %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps);
1181 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1182 %! assert (norm (vec (Q*R - [A(:,1:2) u A(:,3)]), Inf) < norm (A)*1e1*eps);
1183 %!test
1184 %! [Q,R] = qr (Ac);
1185 %! [Q,R] = qrinsert (Q, R, 3, uc);
1186 %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 1e1*eps);
1187 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1188 %! assert (norm (vec (Q*R - [Ac(:,1:2) uc Ac(:,3)]), Inf) < norm (Ac)*1e1*eps);
1189 %!test
1190 %! x = [0.85082 0.76426 0.42883 ];
1191 %!
1192 %! [Q,R] = qr (A);
1193 %! [Q,R] = qrinsert (Q, R, 3, x, "row");
1194 %! assert (norm (vec (Q'*Q - eye (6)), Inf) < 1e1*eps);
1195 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1196 %! assert (norm (vec (Q*R - [A(1:2,:);x;A(3:5,:)]), Inf) < norm (A)*1e1*eps);
1197 %!test
1198 %! x = [0.20351 + 0.05401i 0.13141 + 0.43708i 0.29808 + 0.08789i ];
1199 %!
1200 %! [Q,R] = qr (Ac);
1201 %! [Q,R] = qrinsert (Q, R, 3, x, "row");
1202 %! assert (norm (vec (Q'*Q - eye (6)), Inf) < 1e1*eps);
1203 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1204 %! assert (norm (vec (Q*R - [Ac(1:2,:);x;Ac(3:5,:)]), Inf) < norm (Ac)*1e1*eps);
1205
1206 %!test
1207 %! [Q,R] = qr (single (A));
1208 %! [Q,R] = qrinsert (Q, R, 3, single (u));
1209 %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single"));
1210 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1211 %! assert (norm (vec (Q*R - single ([A(:,1:2) u A(:,3)])), Inf) < norm (single (A))*1e1*eps ("single"));
1212 %!test
1213 %! [Q,R] = qr (single (Ac));
1214 %! [Q,R] = qrinsert (Q, R, 3, single (uc));
1215 %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single"));
1216 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1217 %! assert (norm (vec (Q*R - single ([Ac(:,1:2) uc Ac(:,3)])), Inf) < norm (single (Ac))*1e1*eps ("single"));
1218 %!test
1219 %! x = single ([0.85082 0.76426 0.42883 ]);
1220 %!
1221 %! [Q,R] = qr (single (A));
1222 %! [Q,R] = qrinsert (Q, R, 3, x, "row");
1223 %! assert (norm (vec (Q'*Q - eye (6,"single")), Inf) < 1e1*eps ("single"));
1224 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1225 %! assert (norm (vec (Q*R - single ([A(1:2,:);x;A(3:5,:)])), Inf) < norm (single (A))*1e1*eps ("single"));
1226 %!test
1227 %! x = single ([0.20351 + 0.05401i 0.13141 + 0.43708i 0.29808 + 0.08789i ]);
1228 %!
1229 %! [Q,R] = qr (single (Ac));
1230 %! [Q,R] = qrinsert (Q, R, 3, x, "row");
1231 %! assert (norm (vec (Q'*Q - eye (6,"single")), Inf) < 1e1*eps ("single"));
1232 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1233 %! assert (norm (vec (Q*R - single ([Ac(1:2,:);x;Ac(3:5,:)])), Inf) < norm (single (Ac))*1e1*eps ("single"));
1234 */
1235
1236 DEFUN_DLD (qrdelete, args, ,
1237  doc: /* -*- texinfo -*-
1238 @deftypefn {} {[@var{Q1}, @var{R1}] =} qrdelete (@var{Q}, @var{R}, @var{j}, @var{orient})
1239 Given a QR@tie{}factorization of a real or complex matrix
1240 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and
1241 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization of
1242 @w{[A(:,1:j-1), U, A(:,j:n)]},
1243 where @var{u} is a column vector to be inserted into @var{A}
1244 (if @var{orient} is @qcode{"col"}),
1245 or the QR@tie{}factorization of @w{[A(1:j-1,:);X;A(:,j:n)]},
1246 where @var{x} is a row @var{orient} is @qcode{"row"}).
1247 The default value of @var{orient} is @qcode{"col"}.
1248
1249 If @var{orient} is @qcode{"col"}, @var{j} may be an index vector
1250 resulting in the QR@tie{}factorization of a matrix @var{B} such that
1251 @w{A(:,@var{j}) = []} gives @var{B}. Notice that the latter case is done as
1252 a sequence of k deletions; thus, for k large enough, it will be both faster
1253 and more accurate to recompute the factorization from scratch.
1254
1255 If @var{orient} is @qcode{"col"}, the QR@tie{}factorization supplied may
1256 be either full (Q is square) or economized (R is square).
1257
1258 If @var{orient} is @qcode{"row"}, full factorization is needed.
1259 @seealso{qr, qrupdate, qrinsert, qrshift}
1260 @end deftypefn */)
1261 {
1262  int nargin = args.length ();
1263
1264  if (nargin < 3 || nargin > 4)
1265  print_usage ();
1266
1267  octave_value argq = args(0);
1268  octave_value argr = args(1);
1269  octave_value argj = args(2);
1270
1271  if (! argq.is_numeric_type () || ! argr.is_numeric_type ()
1272  || (nargin > 3 && ! args(3).is_string ()))
1273  print_usage ();
1274
1275  std::string orient = (nargin < 4) ? "col" : args(3).string_value ();
1276  bool col = orient == "col";
1277
1278  if (! col && orient != "row")
1279  error ("qrdelete: ORIENT must be \"col\" or \"row\"");
1280
1281  if (! check_qr_dims (argq, argr, col))
1282  error ("qrdelete: dimension mismatch");
1283
1285  if (! check_index (argj, col))
1286  error ("qrdelete: invalid index J");
1287
1289
1290  octave_idx_type one = 1;
1291
1292  if (argq.is_real_type () && argr.is_real_type ())
1293  {
1294  // real case
1295  if (argq.is_single_type () || argr.is_single_type ())
1296  {
1297  FloatMatrix Q = argq.float_matrix_value ();
1298  FloatMatrix R = argr.float_matrix_value ();
1299
1300  octave::math::qr<FloatMatrix> fact (Q, R);
1301
1302  if (col)
1303  fact.delete_col (j-one);
1304  else
1305  fact.delete_row (j(0)-one);
1306
1307  retval = ovl (fact.Q (), get_qr_r (fact));
1308  }
1309  else
1310  {
1311  Matrix Q = argq.matrix_value ();
1312  Matrix R = argr.matrix_value ();
1313
1314  octave::math::qr<Matrix> fact (Q, R);
1315
1316  if (col)
1317  fact.delete_col (j-one);
1318  else
1319  fact.delete_row (j(0)-one);
1320
1321  retval = ovl (fact.Q (), get_qr_r (fact));
1322  }
1323  }
1324  else
1325  {
1326  // complex case
1327  if (argq.is_single_type () || argr.is_single_type ())
1328  {
1329  FloatComplexMatrix Q =
1331  FloatComplexMatrix R =
1333
1335
1336  if (col)
1337  fact.delete_col (j-one);
1338  else
1339  fact.delete_row (j(0)-one);
1340
1341  retval = ovl (fact.Q (), get_qr_r (fact));
1342  }
1343  else
1344  {
1345  ComplexMatrix Q = argq.complex_matrix_value ();
1346  ComplexMatrix R = argr.complex_matrix_value ();
1347
1348  octave::math::qr<ComplexMatrix> fact (Q, R);
1349
1350  if (col)
1351  fact.delete_col (j-one);
1352  else
1353  fact.delete_row (j(0)-one);
1354
1355  retval = ovl (fact.Q (), get_qr_r (fact));
1356  }
1357  }
1358
1359  return retval;
1360 }
1361
1362 /*
1363 %!test
1364 %! AA = [0.091364 0.613038 0.027504 0.999083;
1365 %! 0.594638 0.425302 0.562834 0.603537;
1366 %! 0.383594 0.291238 0.742073 0.085574;
1367 %! 0.265712 0.268003 0.783553 0.238409;
1368 %! 0.669966 0.743851 0.457255 0.445057 ];
1369 %!
1370 %! [Q,R] = qr (AA);
1371 %! [Q,R] = qrdelete (Q, R, 3);
1372 %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 16*eps);
1373 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1374 %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps);
1375 %!
1376 %!test
1377 %! AA = [0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i;
1378 %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i;
1379 %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i;
1380 %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i;
1381 %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ] * I;
1382 %!
1383 %! [Q,R] = qr (AA);
1384 %! [Q,R] = qrdelete (Q, R, 3);
1385 %! assert (norm (vec (Q'*Q - eye (5)), Inf) < 16*eps);
1386 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1387 %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps);
1388 %!
1389 %!test
1390 %! AA = [0.091364 0.613038 0.027504 0.999083;
1391 %! 0.594638 0.425302 0.562834 0.603537;
1392 %! 0.383594 0.291238 0.742073 0.085574;
1393 %! 0.265712 0.268003 0.783553 0.238409;
1394 %! 0.669966 0.743851 0.457255 0.445057 ];
1395 %!
1396 %! [Q,R] = qr (AA);
1397 %! [Q,R] = qrdelete (Q, R, 3, "row");
1398 %! assert (norm (vec (Q'*Q - eye (4)), Inf) < 1e1*eps);
1399 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1400 %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps);
1401 %!
1402 %!test
1403 %! AA = [0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i;
1404 %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i;
1405 %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i;
1406 %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i;
1407 %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ] * I;
1408 %!
1409 %! [Q,R] = qr (AA);
1410 %! [Q,R] = qrdelete (Q, R, 3, "row");
1411 %! assert (norm (vec (Q'*Q - eye (4)), Inf) < 1e1*eps);
1412 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1413 %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps);
1414
1415 %!test
1416 %! AA = single ([0.091364 0.613038 0.027504 0.999083;
1417 %! 0.594638 0.425302 0.562834 0.603537;
1418 %! 0.383594 0.291238 0.742073 0.085574;
1419 %! 0.265712 0.268003 0.783553 0.238409;
1420 %! 0.669966 0.743851 0.457255 0.445057 ]);
1421 %!
1422 %! [Q,R] = qr (AA);
1423 %! [Q,R] = qrdelete (Q, R, 3);
1424 %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single"));
1425 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1426 %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps ("single"));
1427 %!
1428 %!test
1429 %! AA = single ([0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i;
1430 %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i;
1431 %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i;
1432 %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i;
1433 %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ]) * I;
1434 %!
1435 %! [Q,R] = qr (AA);
1436 %! [Q,R] = qrdelete (Q, R, 3);
1437 %! assert (norm (vec (Q'*Q - eye (5,"single")), Inf) < 1e1*eps ("single"));
1438 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1439 %! assert (norm (vec (Q*R - [AA(:,1:2) AA(:,4)]), Inf) < norm (AA)*1e1*eps ("single"));
1440 %!
1441 %!test
1442 %! AA = single ([0.091364 0.613038 0.027504 0.999083;
1443 %! 0.594638 0.425302 0.562834 0.603537;
1444 %! 0.383594 0.291238 0.742073 0.085574;
1445 %! 0.265712 0.268003 0.783553 0.238409;
1446 %! 0.669966 0.743851 0.457255 0.445057 ]);
1447 %!
1448 %! [Q,R] = qr (AA);
1449 %! [Q,R] = qrdelete (Q, R, 3, "row");
1450 %! assert (norm (vec (Q'*Q - eye (4,"single")), Inf) < 1.5e1*eps ("single"));
1451 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1452 %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps ("single"));
1453 %!testif HAVE_QRUPDATE
1454 %! ## Same test as above but with more precicision
1455 %! AA = single ([0.091364 0.613038 0.027504 0.999083;
1456 %! 0.594638 0.425302 0.562834 0.603537;
1457 %! 0.383594 0.291238 0.742073 0.085574;
1458 %! 0.265712 0.268003 0.783553 0.238409;
1459 %! 0.669966 0.743851 0.457255 0.445057 ]);
1460 %!
1461 %! [Q,R] = qr (AA);
1462 %! [Q,R] = qrdelete (Q, R, 3, "row");
1463 %! assert (norm (vec (Q'*Q - eye (4,"single")), Inf) < 1e1*eps ("single"));
1464 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1465 %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps ("single"));
1466 %!
1467 %!test
1468 %! AA = single ([0.364554 + 0.993117i 0.669818 + 0.510234i 0.426568 + 0.041337i 0.847051 + 0.233291i;
1469 %! 0.049600 + 0.242783i 0.448946 + 0.484022i 0.141155 + 0.074420i 0.446746 + 0.392706i;
1470 %! 0.581922 + 0.657416i 0.581460 + 0.030016i 0.219909 + 0.447288i 0.201144 + 0.069132i;
1471 %! 0.694986 + 0.000571i 0.682327 + 0.841712i 0.807537 + 0.166086i 0.192767 + 0.358098i;
1472 %! 0.945002 + 0.066788i 0.350492 + 0.642638i 0.579629 + 0.048102i 0.600170 + 0.636938i ]) * I;
1473 %!
1474 %! [Q,R] = qr (AA);
1475 %! [Q,R] = qrdelete (Q, R, 3, "row");
1476 %! assert (norm (vec (Q'*Q - eye (4,"single")), Inf) < 1e1*eps ("single"));
1477 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1478 %! assert (norm (vec (Q*R - [AA(1:2,:);AA(4:5,:)]), Inf) < norm (AA)*1e1*eps ("single"));
1479 */
1480
1481 DEFUN_DLD (qrshift, args, ,
1482  doc: /* -*- texinfo -*-
1483 @deftypefn {} {[@var{Q1}, @var{R1}] =} qrshift (@var{Q}, @var{R}, @var{i}, @var{j})
1484 Given a QR@tie{}factorization of a real or complex matrix
1485 @w{@var{A} = @var{Q}*@var{R}}, @var{Q}@tie{}unitary and
1486 @var{R}@tie{}upper trapezoidal, return the QR@tie{}factorization
1487 of @w{@var{A}(:,p)}, where @w{p} is the permutation @*
1488 @code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*
1489  or @*
1490 @code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*
1491
1492 @seealso{qr, qrupdate, qrinsert, qrdelete}
1493 @end deftypefn */)
1494 {
1495  if (args.length () != 4)
1496  print_usage ();
1497
1498  octave_value argq = args(0);
1499  octave_value argr = args(1);
1500  octave_value argi = args(2);
1501  octave_value argj = args(3);
1502
1503  if (! argq.is_numeric_type () || ! argr.is_numeric_type ())
1504  print_usage ();
1505
1506  if (! check_qr_dims (argq, argr, true))
1507  error ("qrshift: dimensions mismatch");
1508
1509  octave_idx_type i = argi.idx_type_value ();
1510  octave_idx_type j = argj.idx_type_value ();
1511
1512  if (! check_index (argi) || ! check_index (argj))
1513  error ("qrshift: invalid index I or J");
1514
1516
1517  if (argq.is_real_type () && argr.is_real_type ())
1518  {
1519  // all real case
1520  if (argq.is_single_type ()
1521  && argr.is_single_type ())
1522  {
1523  FloatMatrix Q = argq.float_matrix_value ();
1524  FloatMatrix R = argr.float_matrix_value ();
1525
1526  octave::math::qr<FloatMatrix> fact (Q, R);
1527  fact.shift_cols (i-1, j-1);
1528
1529  retval = ovl (fact.Q (), get_qr_r (fact));
1530  }
1531  else
1532  {
1533  Matrix Q = argq.matrix_value ();
1534  Matrix R = argr.matrix_value ();
1535
1536  octave::math::qr<Matrix> fact (Q, R);
1537  fact.shift_cols (i-1, j-1);
1538
1539  retval = ovl (fact.Q (), get_qr_r (fact));
1540  }
1541  }
1542  else
1543  {
1544  // complex case
1545  if (argq.is_single_type ()
1546  && argr.is_single_type ())
1547  {
1550
1552  fact.shift_cols (i-1, j-1);
1553
1554  retval = ovl (fact.Q (), get_qr_r (fact));
1555  }
1556  else
1557  {
1558  ComplexMatrix Q = argq.complex_matrix_value ();
1559  ComplexMatrix R = argr.complex_matrix_value ();
1560
1561  octave::math::qr<ComplexMatrix> fact (Q, R);
1562  fact.shift_cols (i-1, j-1);
1563
1564  retval = ovl (fact.Q (), get_qr_r (fact));
1565  }
1566  }
1567
1568  return retval;
1569 }
1570
1571 /*
1572 %!test
1573 %! AA = A.';
1574 %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
1575 %!
1576 %! [Q,R] = qr (AA);
1577 %! [Q,R] = qrshift (Q, R, i, j);
1578 %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps);
1579 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1580 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps);
1581 %!
1582 %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
1583 %!
1584 %! [Q,R] = qr (AA);
1585 %! [Q,R] = qrshift (Q, R, i, j);
1586 %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps);
1587 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1588 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps);
1589 %!
1590 %!test
1591 %! AA = Ac.';
1592 %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
1593 %!
1594 %! [Q,R] = qr (AA);
1595 %! [Q,R] = qrshift (Q, R, i, j);
1596 %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps);
1597 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1598 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps);
1599 %!
1600 %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
1601 %!
1602 %! [Q,R] = qr (AA);
1603 %! [Q,R] = qrshift (Q, R, i, j);
1604 %! assert (norm (vec (Q'*Q - eye (3)), Inf) < 1e1*eps);
1605 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1606 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps);
1607
1608 %!test
1609 %! AA = single (A).';
1610 %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
1611 %!
1612 %! [Q,R] = qr (AA);
1613 %! [Q,R] = qrshift (Q, R, i, j);
1614 %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single"));
1615 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1616 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single"));
1617 %!
1618 %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
1619 %!
1620 %! [Q,R] = qr (AA);
1621 %! [Q,R] = qrshift (Q, R, i, j);
1622 %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single"));
1623 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1624 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single"));
1625 %!
1626 %!test
1627 %! AA = single (Ac).';
1628 %! i = 2; j = 4; p = [1:i-1, shift(i:j,-1), j+1:5];
1629 %!
1630 %! [Q,R] = qr (AA);
1631 %! [Q,R] = qrshift (Q, R, i, j);
1632 %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single"));
1633 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1634 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single"));
1635 %!
1636 %! j = 2; i = 4; p = [1:j-1, shift(j:i,+1), i+1:5];
1637 %!
1638 %! [Q,R] = qr (AA);
1639 %! [Q,R] = qrshift (Q, R, i, j);
1640 %! assert (norm (vec (Q'*Q - eye (3,"single")), Inf) < 1e1*eps ("single"));
1641 %! assert (norm (vec (triu (R) - R), Inf) == 0);
1642 %! assert (norm (vec (Q*R - AA(:,p)), Inf) < norm (AA)*1e1*eps ("single"));
1643 */
FloatComplexMatrix transpose(void) const
Definition: fCMatrix.h:170
bool is_real_type(void) const
Definition: ov.h:667
T Q(void) const
Definition: qr.h:78
int ndims(void) const
Definition: ov.h:495
octave_idx_type rows(void) const
Definition: ov.h:489
static octave::math::qr< T >::type qr_type(int nargin, int nargout)
Definition: qr.cc:51
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
identity matrix If supplied two scalar respectively For allows like xample val
Definition: data.cc:5068
bool is_scalar_type(void) const
Definition: ov.h:673
FloatMatrix transpose(void) const
Definition: fMatrix.h:133
bool is_numeric_type(void) const
Definition: ov.h:679
static octave_value get_qr_r(const octave::math::qr< MT > &fact)
Definition: qr.cc:40
for large enough k
Definition: lu.cc:606
bool regular(void) const
Definition: qr.cc:85
void error(const char *fmt,...)
Definition: error.cc:570
RowVector row(octave_idx_type i) const
Definition: dMatrix.cc:417
ComplexRowVector row(octave_idx_type i) const
Definition: CMatrix.cc:705
u
Definition: lu.cc:138
FloatRowVector row(octave_idx_type i) const
Definition: fMatrix.cc:423
octave_idx_type idx_type_value(bool req_int=false, bool frc_str_conv=false) const
Definition: ov.cc:1677
RV_T Pvec(void) const
octave_value arg
Definition: pr-output.cc:3440
string_vector argv
JNIEnv void * args
Definition: ov-java.cc:67
octave_idx_type columns(void) const
Definition: ov.h:491
ComplexColumnVector conj(const ComplexColumnVector &a)
Definition: CColVector.cc:216
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
bool is_sparse_type(void) const
Definition: ov.h:682
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
int nargin
Definition: graphics.cc:10115
Array< octave_idx_type > octave_idx_type_vector_value(bool req_int=false, bool frc_str_conv=false, bool frc_vec_conv=false) const
Definition: ov.cc:1891
bool is_complex_type(void) const
Definition: ov.h:670
octave_value retval
Definition: data.cc:6294
Matrix transpose(void) const
Definition: dMatrix.h:129
idx type
Definition: ov.cc:3129
Definition: dMatrix.h:37
ComplexMatrix transpose(void) const
Definition: CMatrix.h:165
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
void warning(const char *fmt,...)
Definition: error.cc:788
PermMatrix P(void) const
Definition: qrp.h:67
Definition: mx-defs.h:77
=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
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:790
#define DEFUN_DLD(name, args_name, nargout_name, doc)
Definition: defun-dld.h:45
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:834
T R(void) const
Definition: qr.h:80
bool is_single_type(void) const
Definition: ov.h:627
FloatComplexRowVector row(octave_idx_type i) const
Definition: fCMatrix.cc:708
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
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
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
bool is_integer_type(void) const
Definition: ov.h:664