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