GNU Octave  4.2.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
qr.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2017 John W. Eaton
4 Copyright (C) 2008-2009 Jaroslav Hajek
5 Copyright (C) 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 #if defined (HAVE_CONFIG_H)
26 # include "config.h"
27 #endif
28 
29 #include "CColVector.h"
30 #include "CMatrix.h"
31 #include "CRowVector.h"
32 #include "dColVector.h"
33 #include "dMatrix.h"
34 #include "dRowVector.h"
35 #include "f77-fcn.h"
36 #include "fCColVector.h"
37 #include "fCMatrix.h"
38 #include "fCRowVector.h"
39 #include "fColVector.h"
40 #include "fMatrix.h"
41 #include "fRowVector.h"
42 #include "idx-vector.h"
43 #include "lo-error.h"
44 #include "lo-lapack-proto.h"
45 #include "lo-qrupdate-proto.h"
46 #include "oct-locbuf.h"
47 #include "qr.h"
48 
49 namespace octave
50 {
51  namespace math
52  {
53  template <typename T>
54  qr<T>::qr (const T& q_arg, const T& r_arg)
55  : q (q_arg), r (r_arg)
56  {
57  octave_idx_type q_nr = q.rows ();
58  octave_idx_type q_nc = q.columns ();
59 
60  octave_idx_type r_nr = r.rows ();
61  octave_idx_type r_nc = r.columns ();
62 
63  if (! (q_nc == r_nr && (q_nr == q_nc || (q_nr > q_nc && r_nr == r_nc))))
64  (*current_liboctave_error_handler) ("QR dimensions mismatch");
65  }
66 
67  template <typename T>
68  typename qr<T>::type
69  qr<T>::get_type (void) const
70  {
71  type retval;
72 
73  if (! q.is_empty () && q.is_square ())
74  retval = qr<T>::std;
75  else if (q.rows () > q.columns () && r.is_square ())
76  retval = qr<T>::economy;
77  else
78  retval = qr<T>::raw;
79 
80  return retval;
81  }
82 
83  template <typename T>
84  bool
85  qr<T>::regular (void) const
86  {
87  bool retval = true;
88 
89  octave_idx_type k = std::min (r.rows (), r.columns ());
90 
91  for (octave_idx_type i = 0; i < k; i++)
92  {
93  if (r(i, i) == ELT_T ())
94  {
95  retval = false;
96  break;
97  }
98  }
99 
100  return retval;
101  }
102 
103 #if ! defined (HAVE_QRUPDATE)
104 
105  // Replacement update methods.
106 
107  void
108  warn_qrupdate_once (void)
109  {
110  static bool warned = false;
111 
112  if (! warned)
113  {
114  (*current_liboctave_warning_with_id_handler)
115  ("Octave:missing-dependency",
116  "In this version of Octave, QR & Cholesky updating routines "
117  "simply update the matrix and recalculate factorizations. "
118  "To use fast algorithms, link Octave with the qrupdate library. "
119  "See <http://sourceforge.net/projects/qrupdate>.");
120 
121  warned = true;
122  }
123  }
124 
125  template <typename T>
126  void
127  qr<T>::update (const CV_T& u, const CV_T& v)
128  {
130 
131  octave_idx_type m = q.rows ();
132  octave_idx_type n = r.columns ();
133 
134  if (u.numel () != m || v.numel () != n)
135  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
136 
137  init (q*r + T (u) * T (v).hermitian (), get_type ());
138  }
139 
140  template <typename T>
141  void
142  qr<T>::update (const T& u, const T& v)
143  {
145 
146  octave_idx_type m = q.rows ();
147  octave_idx_type n = r.columns ();
148 
149  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
150  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
151 
152  init (q*r + u * v.hermitian (), get_type ());
153  }
154 
155  template <typename T, typename CV_T>
156  static
157  T
158  insert_col (const T& a, octave_idx_type i, const CV_T& x)
159  {
160  T retval (a.rows (), a.columns () + 1);
162  a.index (idx_vector::colon, idx_vector (0, i)));
165  a.index (idx_vector::colon, idx_vector (i, a.columns ())));
166  return retval;
167  }
168 
169  template <typename T, typename RV_T>
170  static
171  T
172  insert_row (const T& a, octave_idx_type i, const RV_T& x)
173  {
174  T retval (a.rows () + 1, a.columns ());
176  a.index (idx_vector (0, i), idx_vector::colon));
179  a.index (idx_vector (i, a.rows ()), idx_vector::colon));
180  return retval;
181  }
182 
183  template <typename T>
184  static
185  T
186  delete_col (const T& a, octave_idx_type i)
187  {
188  T retval = a;
189  retval.delete_elements (1, idx_vector (i));
190  return retval;
191  }
192 
193  template <typename T>
194  static
195  T
196  delete_row (const T& a, octave_idx_type i)
197  {
198  T retval = a;
199  retval.delete_elements (0, idx_vector (i));
200  return retval;
201  }
202 
203  template <typename T>
204  static
205  T
206  shift_cols (const T& a, octave_idx_type i, octave_idx_type j)
207  {
208  octave_idx_type n = a.columns ();
210  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
211  if (i < j)
212  {
213  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
214  p(j) = i;
215  }
216  else if (j < i)
217  {
218  p(j) = i;
219  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
220  }
221 
222  return a.index (idx_vector::colon, idx_vector (p));
223  }
224 
225  template <typename T>
226  void
227  qr<T>::insert_col (const CV_T& u, octave_idx_type j)
228  {
230 
231  octave_idx_type m = q.rows ();
232  octave_idx_type n = r.columns ();
233 
234  if (u.numel () != m)
235  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
236  if (j < 0 || j > n)
237  (*current_liboctave_error_handler) ("qrinsert: index out of range");
238 
239  init (octave::math::insert_col (q*r, j, u), get_type ());
240  }
241 
242  template <typename T>
243  void
244  qr<T>::insert_col (const T& u, const Array<octave_idx_type>& j)
245  {
247 
248  octave_idx_type m = q.rows ();
249  octave_idx_type n = r.columns ();
250 
252  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
253  octave_idx_type nj = js.numel ();
254  bool dups = false;
255  for (octave_idx_type i = 0; i < nj - 1; i++)
256  dups = dups && js(i) == js(i+1);
257 
258  if (dups)
259  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
260  if (u.numel () != m || u.columns () != nj)
261  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
262  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
263  (*current_liboctave_error_handler) ("qrinsert: index out of range");
264 
265  if (nj > 0)
266  {
267  T a = q*r;
268  for (octave_idx_type i = 0; i < js.numel (); i++)
269  a = octave::math::insert_col (a, js(i), u.column (i));
270  init (a, get_type ());
271  }
272  }
273 
274  template <typename T>
275  void
277  {
279 
280  octave_idx_type n = r.columns ();
281 
282  if (j < 0 || j > n-1)
283  (*current_liboctave_error_handler) ("qrdelete: index out of range");
284 
285  init (octave::math::delete_col (q*r, j), get_type ());
286  }
287 
288  template <typename T>
289  void
291  {
293 
294  octave_idx_type n = r.columns ();
295 
297  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
298  octave_idx_type nj = js.numel ();
299  bool dups = false;
300  for (octave_idx_type i = 0; i < nj - 1; i++)
301  dups = dups && js(i) == js(i+1);
302 
303  if (dups)
304  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
305  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
306  (*current_liboctave_error_handler) ("qrinsert: index out of range");
307 
308  if (nj > 0)
309  {
310  T a = q*r;
311  for (octave_idx_type i = 0; i < js.numel (); i++)
312  a = octave::math::delete_col (a, js(i));
313  init (a, get_type ());
314  }
315  }
316 
317  template <typename T>
318  void
319  qr<T>::insert_row (const RV_T& u, octave_idx_type j)
320  {
322 
323  octave_idx_type m = r.rows ();
324  octave_idx_type n = r.columns ();
325 
326  if (! q.is_square () || u.numel () != n)
327  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
328  if (j < 0 || j > m)
329  (*current_liboctave_error_handler) ("qrinsert: index out of range");
330 
331  init (octave::math::insert_row (q*r, j, u), get_type ());
332  }
333 
334  template <typename T>
335  void
337  {
339 
340  octave_idx_type m = r.rows ();
341 
342  if (! q.is_square ())
343  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
344  if (j < 0 || j > m-1)
345  (*current_liboctave_error_handler) ("qrdelete: index out of range");
346 
347  init (octave::math::delete_row (q*r, j), get_type ());
348  }
349 
350  template <typename T>
351  void
353  {
355 
356  octave_idx_type n = r.columns ();
357 
358  if (i < 0 || i > n-1 || j < 0 || j > n-1)
359  (*current_liboctave_error_handler) ("qrshift: index out of range");
360 
361  init (octave::math::shift_cols (q*r, i, j), get_type ());
362  }
363 
364 #endif
365 
366  // Specializations.
367 
368  template <>
369  void
371  {
372  octave_idx_type m = afact.rows ();
373  octave_idx_type min_mn = std::min (m, n);
374  octave_idx_type info;
375 
376  if (qr_type == qr<Matrix>::raw)
377  {
378  for (octave_idx_type j = 0; j < min_mn; j++)
379  {
380  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
381  for (octave_idx_type i = limit + 1; i < m; i++)
382  afact.elem (i, j) *= tau[j];
383  }
384 
385  r = afact;
386  }
387  else
388  {
389  // Attempt to minimize copying.
390  if (m >= n)
391  {
392  // afact will become q.
393  q = afact;
394  octave_idx_type k = qr_type == qr<Matrix>::economy ? n : m;
395  r = Matrix (k, n);
396  for (octave_idx_type j = 0; j < n; j++)
397  {
398  octave_idx_type i = 0;
399  for (; i <= j; i++)
400  r.xelem (i, j) = afact.xelem (i, j);
401  for (; i < k; i++)
402  r.xelem (i, j) = 0;
403  }
404  afact = Matrix (); // optimize memory
405  }
406  else
407  {
408  // afact will become r.
409  q = Matrix (m, m);
410  for (octave_idx_type j = 0; j < m; j++)
411  for (octave_idx_type i = j + 1; i < m; i++)
412  {
413  q.xelem (i, j) = afact.xelem (i, j);
414  afact.xelem (i, j) = 0;
415  }
416  r = afact;
417  }
418 
419  if (m > 0)
420  {
421  octave_idx_type k = q.columns ();
422  // workspace query.
423  double rlwork;
424  F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
425  &rlwork, -1, info));
426 
427  // allocate buffer and do the job.
428  octave_idx_type lwork = rlwork;
429  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
430  OCTAVE_LOCAL_BUFFER (double, work, lwork);
431  F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
432  work, lwork, info));
433  }
434  }
435  }
436 
437  template <>
438  void
440  {
441  octave_idx_type m = a.rows ();
442  octave_idx_type n = a.cols ();
443 
444  octave_idx_type min_mn = m < n ? m : n;
445  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
446 
447  octave_idx_type info = 0;
448 
449  Matrix afact = a;
450  if (m > n && qr_type == qr<Matrix>::std)
451  afact.resize (m, m);
452 
453  if (m > 0)
454  {
455  // workspace query.
456  double rlwork;
457  F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
458  &rlwork, -1, info));
459 
460  // allocate buffer and do the job.
461  octave_idx_type lwork = rlwork;
462  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
463  OCTAVE_LOCAL_BUFFER (double, work, lwork);
464  F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
465  work, lwork, info));
466  }
467 
468  form (n, afact, tau, qr_type);
469  }
470 
471 #if defined (HAVE_QRUPDATE)
472 
473  template <>
474  void
476  {
477  octave_idx_type m = q.rows ();
478  octave_idx_type n = r.columns ();
479  octave_idx_type k = q.columns ();
480 
481  if (u.numel () != m || v.numel () != n)
482  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
483 
484  ColumnVector utmp = u;
485  ColumnVector vtmp = v;
486  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
487  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
488  m, r.fortran_vec (), k,
489  utmp.fortran_vec (), vtmp.fortran_vec (), w));
490  }
491 
492  template <>
493  void
494  qr<Matrix>::update (const Matrix& u, const Matrix& v)
495  {
496  octave_idx_type m = q.rows ();
497  octave_idx_type n = r.columns ();
498  octave_idx_type k = q.columns ();
499 
500  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
501  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
502 
503  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
504  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
505  {
506  ColumnVector utmp = u.column (i);
507  ColumnVector vtmp = v.column (i);
508  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (),
509  m, r.fortran_vec (), k,
510  utmp.fortran_vec (), vtmp.fortran_vec (),
511  w));
512  }
513  }
514 
515  template <>
516  void
518  {
519  octave_idx_type m = q.rows ();
520  octave_idx_type n = r.columns ();
521  octave_idx_type k = q.columns ();
522 
523  if (u.numel () != m)
524  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
525  if (j < 0 || j > n)
526  (*current_liboctave_error_handler) ("qrinsert: index out of range");
527 
528  if (k < m)
529  {
530  q.resize (m, k+1);
531  r.resize (k+1, n+1);
532  }
533  else
534  {
535  r.resize (k, n+1);
536  }
537 
538  ColumnVector utmp = u;
539  OCTAVE_LOCAL_BUFFER (double, w, k);
540  F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), q.rows (),
541  r.fortran_vec (), r.rows (), j + 1,
542  utmp.data (), w));
543  }
544 
545  template <>
546  void
548  {
549  octave_idx_type m = q.rows ();
550  octave_idx_type n = r.columns ();
551  octave_idx_type k = q.columns ();
552 
554  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
555  octave_idx_type nj = js.numel ();
556  bool dups = false;
557  for (octave_idx_type i = 0; i < nj - 1; i++)
558  dups = dups && js(i) == js(i+1);
559 
560  if (dups)
561  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
562  if (u.numel () != m || u.columns () != nj)
563  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
564  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
565  (*current_liboctave_error_handler) ("qrinsert: index out of range");
566 
567  if (nj > 0)
568  {
569  octave_idx_type kmax = std::min (k + nj, m);
570  if (k < m)
571  {
572  q.resize (m, kmax);
573  r.resize (kmax, n + nj);
574  }
575  else
576  {
577  r.resize (k, n + nj);
578  }
579 
580  OCTAVE_LOCAL_BUFFER (double, w, kmax);
581  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
582  {
583  octave_idx_type ii = i;
584  ColumnVector utmp = u.column (jsi(i));
585  F77_XFCN (dqrinc, DQRINC, (m, n + ii, std::min (kmax, k + ii),
586  q.fortran_vec (), q.rows (),
587  r.fortran_vec (), r.rows (), js(ii) + 1,
588  utmp.data (), w));
589  }
590  }
591  }
592 
593  template <>
594  void
596  {
597  octave_idx_type m = q.rows ();
598  octave_idx_type k = r.rows ();
599  octave_idx_type n = r.columns ();
600 
601  if (j < 0 || j > n-1)
602  (*current_liboctave_error_handler) ("qrdelete: index out of range");
603 
604  OCTAVE_LOCAL_BUFFER (double, w, k);
605  F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
606  r.fortran_vec (), r.rows (), j + 1, w));
607 
608  if (k < m)
609  {
610  q.resize (m, k-1);
611  r.resize (k-1, n-1);
612  }
613  else
614  {
615  r.resize (k, n-1);
616  }
617  }
618 
619  template <>
620  void
622  {
623  octave_idx_type m = q.rows ();
624  octave_idx_type n = r.columns ();
625  octave_idx_type k = q.columns ();
626 
628  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
629  octave_idx_type nj = js.numel ();
630  bool dups = false;
631  for (octave_idx_type i = 0; i < nj - 1; i++)
632  dups = dups && js(i) == js(i+1);
633 
634  if (dups)
635  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
636  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
637  (*current_liboctave_error_handler) ("qrinsert: index out of range");
638 
639  if (nj > 0)
640  {
641  OCTAVE_LOCAL_BUFFER (double, w, k);
642  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
643  {
644  octave_idx_type ii = i;
645  F77_XFCN (dqrdec, DQRDEC, (m, n - ii, k == m ? k : k - ii,
646  q.fortran_vec (), q.rows (),
647  r.fortran_vec (), r.rows (),
648  js(ii) + 1, w));
649  }
650  if (k < m)
651  {
652  q.resize (m, k - nj);
653  r.resize (k - nj, n - nj);
654  }
655  else
656  {
657  r.resize (k, n - nj);
658  }
659 
660  }
661  }
662 
663  template <>
664  void
666  {
667  octave_idx_type m = r.rows ();
668  octave_idx_type n = r.columns ();
669  octave_idx_type k = std::min (m, n);
670 
671  if (! q.is_square () || u.numel () != n)
672  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
673  if (j < 0 || j > m)
674  (*current_liboctave_error_handler) ("qrinsert: index out of range");
675 
676  q.resize (m + 1, m + 1);
677  r.resize (m + 1, n);
678  RowVector utmp = u;
679  OCTAVE_LOCAL_BUFFER (double, w, k);
680  F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), q.rows (),
681  r.fortran_vec (), r.rows (),
682  j + 1, utmp.fortran_vec (), w));
683 
684  }
685 
686  template <>
687  void
689  {
690  octave_idx_type m = r.rows ();
691  octave_idx_type n = r.columns ();
692 
693  if (! q.is_square ())
694  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
695  if (j < 0 || j > m-1)
696  (*current_liboctave_error_handler) ("qrdelete: index out of range");
697 
698  OCTAVE_LOCAL_BUFFER (double, w, 2*m);
699  F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), q.rows (),
700  r.fortran_vec (), r.rows (), j + 1,
701  w));
702 
703  q.resize (m - 1, m - 1);
704  r.resize (m - 1, n);
705  }
706 
707  template <>
708  void
710  {
711  octave_idx_type m = q.rows ();
712  octave_idx_type k = r.rows ();
713  octave_idx_type n = r.columns ();
714 
715  if (i < 0 || i > n-1 || j < 0 || j > n-1)
716  (*current_liboctave_error_handler) ("qrshift: index out of range");
717 
718  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
719  F77_XFCN (dqrshc, DQRSHC, (m, n, k,
720  q.fortran_vec (), q.rows (),
721  r.fortran_vec (), r.rows (),
722  i + 1, j + 1, w));
723  }
724 
725 #endif
726 
727  template <>
728  void
730  type qr_type)
731  {
732  octave_idx_type m = afact.rows ();
733  octave_idx_type min_mn = std::min (m, n);
734  octave_idx_type info;
735 
736  if (qr_type == qr<FloatMatrix>::raw)
737  {
738  for (octave_idx_type j = 0; j < min_mn; j++)
739  {
740  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
741  for (octave_idx_type i = limit + 1; i < m; i++)
742  afact.elem (i, j) *= tau[j];
743  }
744 
745  r = afact;
746  }
747  else
748  {
749  // Attempt to minimize copying.
750  if (m >= n)
751  {
752  // afact will become q.
753  q = afact;
754  octave_idx_type k = qr_type == qr<FloatMatrix>::economy ? n : m;
755  r = FloatMatrix (k, n);
756  for (octave_idx_type j = 0; j < n; j++)
757  {
758  octave_idx_type i = 0;
759  for (; i <= j; i++)
760  r.xelem (i, j) = afact.xelem (i, j);
761  for (; i < k; i++)
762  r.xelem (i, j) = 0;
763  }
764  afact = FloatMatrix (); // optimize memory
765  }
766  else
767  {
768  // afact will become r.
769  q = FloatMatrix (m, m);
770  for (octave_idx_type j = 0; j < m; j++)
771  for (octave_idx_type i = j + 1; i < m; i++)
772  {
773  q.xelem (i, j) = afact.xelem (i, j);
774  afact.xelem (i, j) = 0;
775  }
776  r = afact;
777  }
778 
779  if (m > 0)
780  {
781  octave_idx_type k = q.columns ();
782  // workspace query.
783  float rlwork;
784  F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
785  &rlwork, -1, info));
786 
787  // allocate buffer and do the job.
788  octave_idx_type lwork = rlwork;
789  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
790  OCTAVE_LOCAL_BUFFER (float, work, lwork);
791  F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m, tau,
792  work, lwork, info));
793  }
794  }
795  }
796 
797  template <>
798  void
800  {
801  octave_idx_type m = a.rows ();
802  octave_idx_type n = a.cols ();
803 
804  octave_idx_type min_mn = m < n ? m : n;
805  OCTAVE_LOCAL_BUFFER (float, tau, min_mn);
806 
807  octave_idx_type info = 0;
808 
809  FloatMatrix afact = a;
810  if (m > n && qr_type == qr<FloatMatrix>::std)
811  afact.resize (m, m);
812 
813  if (m > 0)
814  {
815  // workspace query.
816  float rlwork;
817  F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
818  &rlwork, -1, info));
819 
820  // allocate buffer and do the job.
821  octave_idx_type lwork = rlwork;
822  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
823  OCTAVE_LOCAL_BUFFER (float, work, lwork);
824  F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
825  work, lwork, info));
826  }
827 
828  form (n, afact, tau, qr_type);
829  }
830 
831 #if defined (HAVE_QRUPDATE)
832 
833  template <>
834  void
836  {
837  octave_idx_type m = q.rows ();
838  octave_idx_type n = r.columns ();
839  octave_idx_type k = q.columns ();
840 
841  if (u.numel () != m || v.numel () != n)
842  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
843 
844  FloatColumnVector utmp = u;
845  FloatColumnVector vtmp = v;
846  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
847  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
848  m, r.fortran_vec (), k,
849  utmp.fortran_vec (), vtmp.fortran_vec (), w));
850  }
851 
852  template <>
853  void
855  {
856  octave_idx_type m = q.rows ();
857  octave_idx_type n = r.columns ();
858  octave_idx_type k = q.columns ();
859 
860  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
861  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
862 
863  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
864  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
865  {
866  FloatColumnVector utmp = u.column (i);
867  FloatColumnVector vtmp = v.column (i);
868  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (),
869  m, r.fortran_vec (), k,
870  utmp.fortran_vec (), vtmp.fortran_vec (),
871  w));
872  }
873  }
874 
875  template <>
876  void
878  {
879  octave_idx_type m = q.rows ();
880  octave_idx_type n = r.columns ();
881  octave_idx_type k = q.columns ();
882 
883  if (u.numel () != m)
884  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
885  if (j < 0 || j > n)
886  (*current_liboctave_error_handler) ("qrinsert: index out of range");
887 
888  if (k < m)
889  {
890  q.resize (m, k+1);
891  r.resize (k+1, n+1);
892  }
893  else
894  {
895  r.resize (k, n+1);
896  }
897 
898  FloatColumnVector utmp = u;
899  OCTAVE_LOCAL_BUFFER (float, w, k);
900  F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), q.rows (),
901  r.fortran_vec (), r.rows (), j + 1,
902  utmp.data (), w));
903  }
904 
905  template <>
906  void
908  const Array<octave_idx_type>& j)
909  {
910  octave_idx_type m = q.rows ();
911  octave_idx_type n = r.columns ();
912  octave_idx_type k = q.columns ();
913 
915  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
916  octave_idx_type nj = js.numel ();
917  bool dups = false;
918  for (octave_idx_type i = 0; i < nj - 1; i++)
919  dups = dups && js(i) == js(i+1);
920 
921  if (dups)
922  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
923  if (u.numel () != m || u.columns () != nj)
924  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
925  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
926  (*current_liboctave_error_handler) ("qrinsert: index out of range");
927 
928  if (nj > 0)
929  {
930  octave_idx_type kmax = std::min (k + nj, m);
931  if (k < m)
932  {
933  q.resize (m, kmax);
934  r.resize (kmax, n + nj);
935  }
936  else
937  {
938  r.resize (k, n + nj);
939  }
940 
941  OCTAVE_LOCAL_BUFFER (float, w, kmax);
942  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
943  {
944  octave_idx_type ii = i;
945  FloatColumnVector utmp = u.column (jsi(i));
946  F77_XFCN (sqrinc, SQRINC, (m, n + ii, std::min (kmax, k + ii),
947  q.fortran_vec (), q.rows (),
948  r.fortran_vec (), r.rows (), js(ii) + 1,
949  utmp.data (), w));
950  }
951  }
952  }
953 
954  template <>
955  void
957  {
958  octave_idx_type m = q.rows ();
959  octave_idx_type k = r.rows ();
960  octave_idx_type n = r.columns ();
961 
962  if (j < 0 || j > n-1)
963  (*current_liboctave_error_handler) ("qrdelete: index out of range");
964 
965  OCTAVE_LOCAL_BUFFER (float, w, k);
966  F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), q.rows (),
967  r.fortran_vec (), r.rows (), j + 1, w));
968 
969  if (k < m)
970  {
971  q.resize (m, k-1);
972  r.resize (k-1, n-1);
973  }
974  else
975  {
976  r.resize (k, n-1);
977  }
978  }
979 
980  template <>
981  void
983  {
984  octave_idx_type m = q.rows ();
985  octave_idx_type n = r.columns ();
986  octave_idx_type k = q.columns ();
987 
989  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
990  octave_idx_type nj = js.numel ();
991  bool dups = false;
992  for (octave_idx_type i = 0; i < nj - 1; i++)
993  dups = dups && js(i) == js(i+1);
994 
995  if (dups)
996  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
997  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
998  (*current_liboctave_error_handler) ("qrinsert: index out of range");
999 
1000  if (nj > 0)
1001  {
1002  OCTAVE_LOCAL_BUFFER (float, w, k);
1003  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
1004  {
1005  octave_idx_type ii = i;
1006  F77_XFCN (sqrdec, SQRDEC, (m, n - ii, k == m ? k : k - ii,
1007  q.fortran_vec (), q.rows (),
1008  r.fortran_vec (), r.rows (),
1009  js(ii) + 1, w));
1010  }
1011  if (k < m)
1012  {
1013  q.resize (m, k - nj);
1014  r.resize (k - nj, n - nj);
1015  }
1016  else
1017  {
1018  r.resize (k, n - nj);
1019  }
1020 
1021  }
1022  }
1023 
1024  template <>
1025  void
1027  {
1028  octave_idx_type m = r.rows ();
1029  octave_idx_type n = r.columns ();
1030  octave_idx_type k = std::min (m, n);
1031 
1032  if (! q.is_square () || u.numel () != n)
1033  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1034  if (j < 0 || j > m)
1035  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1036 
1037  q.resize (m + 1, m + 1);
1038  r.resize (m + 1, n);
1039  FloatRowVector utmp = u;
1040  OCTAVE_LOCAL_BUFFER (float, w, k);
1041  F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), q.rows (),
1042  r.fortran_vec (), r.rows (),
1043  j + 1, utmp.fortran_vec (), w));
1044 
1045  }
1046 
1047  template <>
1048  void
1050  {
1051  octave_idx_type m = r.rows ();
1052  octave_idx_type n = r.columns ();
1053 
1054  if (! q.is_square ())
1055  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
1056  if (j < 0 || j > m-1)
1057  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1058 
1059  OCTAVE_LOCAL_BUFFER (float, w, 2*m);
1060  F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), q.rows (),
1061  r.fortran_vec (), r.rows (), j + 1,
1062  w));
1063 
1064  q.resize (m - 1, m - 1);
1065  r.resize (m - 1, n);
1066  }
1067 
1068  template <>
1069  void
1071  {
1072  octave_idx_type m = q.rows ();
1073  octave_idx_type k = r.rows ();
1074  octave_idx_type n = r.columns ();
1075 
1076  if (i < 0 || i > n-1 || j < 0 || j > n-1)
1077  (*current_liboctave_error_handler) ("qrshift: index out of range");
1078 
1079  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
1080  F77_XFCN (sqrshc, SQRSHC, (m, n, k,
1081  q.fortran_vec (), q.rows (),
1082  r.fortran_vec (), r.rows (),
1083  i + 1, j + 1, w));
1084  }
1085 
1086 #endif
1087 
1088  template <>
1089  void
1091  Complex *tau, type qr_type)
1092  {
1093  octave_idx_type m = afact.rows ();
1094  octave_idx_type min_mn = std::min (m, n);
1095  octave_idx_type info;
1096 
1097  if (qr_type == qr<ComplexMatrix>::raw)
1098  {
1099  for (octave_idx_type j = 0; j < min_mn; j++)
1100  {
1101  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
1102  for (octave_idx_type i = limit + 1; i < m; i++)
1103  afact.elem (i, j) *= tau[j];
1104  }
1105 
1106  r = afact;
1107  }
1108  else
1109  {
1110  // Attempt to minimize copying.
1111  if (m >= n)
1112  {
1113  // afact will become q.
1114  q = afact;
1115  octave_idx_type k = qr_type == qr<ComplexMatrix>::economy ? n : m;
1116  r = ComplexMatrix (k, n);
1117  for (octave_idx_type j = 0; j < n; j++)
1118  {
1119  octave_idx_type i = 0;
1120  for (; i <= j; i++)
1121  r.xelem (i, j) = afact.xelem (i, j);
1122  for (; i < k; i++)
1123  r.xelem (i, j) = 0;
1124  }
1125  afact = ComplexMatrix (); // optimize memory
1126  }
1127  else
1128  {
1129  // afact will become r.
1130  q = ComplexMatrix (m, m);
1131  for (octave_idx_type j = 0; j < m; j++)
1132  for (octave_idx_type i = j + 1; i < m; i++)
1133  {
1134  q.xelem (i, j) = afact.xelem (i, j);
1135  afact.xelem (i, j) = 0;
1136  }
1137  r = afact;
1138  }
1139 
1140  if (m > 0)
1141  {
1142  octave_idx_type k = q.columns ();
1143  // workspace query.
1144  Complex clwork;
1145  F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1146  m, F77_DBLE_CMPLX_ARG (tau),
1147  F77_DBLE_CMPLX_ARG (&clwork), -1, info));
1148 
1149  // allocate buffer and do the job.
1150  octave_idx_type lwork = clwork.real ();
1151  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
1152  OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
1153  F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1154  m, F77_DBLE_CMPLX_ARG (tau),
1155  F77_DBLE_CMPLX_ARG (work), lwork, info));
1156  }
1157  }
1158  }
1159 
1160  template <>
1161  void
1163  {
1164  octave_idx_type m = a.rows ();
1165  octave_idx_type n = a.cols ();
1166 
1167  octave_idx_type min_mn = m < n ? m : n;
1168  OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
1169 
1170  octave_idx_type info = 0;
1171 
1172  ComplexMatrix afact = a;
1173  if (m > n && qr_type == qr<ComplexMatrix>::std)
1174  afact.resize (m, m);
1175 
1176  if (m > 0)
1177  {
1178  // workspace query.
1179  Complex clwork;
1180  F77_XFCN (zgeqrf, ZGEQRF, (m, n, F77_DBLE_CMPLX_ARG (afact.fortran_vec ()), m,
1181  F77_DBLE_CMPLX_ARG (tau),
1182  F77_DBLE_CMPLX_ARG (&clwork), -1, info));
1183 
1184  // allocate buffer and do the job.
1185  octave_idx_type lwork = clwork.real ();
1186  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
1187  OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
1188  F77_XFCN (zgeqrf, ZGEQRF, (m, n, F77_DBLE_CMPLX_ARG (afact.fortran_vec ()), m,
1189  F77_DBLE_CMPLX_ARG (tau),
1190  F77_DBLE_CMPLX_ARG (work), lwork, info));
1191  }
1192 
1193  form (n, afact, tau, qr_type);
1194  }
1195 
1196 #if defined (HAVE_QRUPDATE)
1197 
1198  template <>
1199  void
1201  const ComplexColumnVector& v)
1202  {
1203  octave_idx_type m = q.rows ();
1204  octave_idx_type n = r.columns ();
1205  octave_idx_type k = q.columns ();
1206 
1207  if (u.numel () != m || v.numel () != n)
1208  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1209 
1210  ComplexColumnVector utmp = u;
1211  ComplexColumnVector vtmp = v;
1213  OCTAVE_LOCAL_BUFFER (double, rw, k);
1214  F77_XFCN (zqr1up, ZQR1UP, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1215  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
1216  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
1217  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()),
1218  F77_DBLE_CMPLX_ARG (w), rw));
1219  }
1220 
1221  template <>
1222  void
1224  {
1225  octave_idx_type m = q.rows ();
1226  octave_idx_type n = r.columns ();
1227  octave_idx_type k = q.columns ();
1228 
1229  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
1230  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1231 
1233  OCTAVE_LOCAL_BUFFER (double, rw, k);
1234  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
1235  {
1236  ComplexColumnVector utmp = u.column (i);
1237  ComplexColumnVector vtmp = v.column (i);
1238  F77_XFCN (zqr1up, ZQR1UP, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1239  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
1240  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
1241  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()),
1242  F77_DBLE_CMPLX_ARG (w), rw));
1243  }
1244  }
1245 
1246  template <>
1247  void
1249  {
1250  octave_idx_type m = q.rows ();
1251  octave_idx_type n = r.columns ();
1252  octave_idx_type k = q.columns ();
1253 
1254  if (u.numel () != m)
1255  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1256  if (j < 0 || j > n)
1257  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1258 
1259  if (k < m)
1260  {
1261  q.resize (m, k+1);
1262  r.resize (k+1, n+1);
1263  }
1264  else
1265  {
1266  r.resize (k, n+1);
1267  }
1268 
1269  ComplexColumnVector utmp = u;
1270  OCTAVE_LOCAL_BUFFER (double, rw, k);
1271  F77_XFCN (zqrinc, ZQRINC, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1272  q.rows (),
1273  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (), j + 1,
1274  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), rw));
1275  }
1276 
1277  template <>
1278  void
1280  const Array<octave_idx_type>& j)
1281  {
1282  octave_idx_type m = q.rows ();
1283  octave_idx_type n = r.columns ();
1284  octave_idx_type k = q.columns ();
1285 
1287  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
1288  octave_idx_type nj = js.numel ();
1289  bool dups = false;
1290  for (octave_idx_type i = 0; i < nj - 1; i++)
1291  dups = dups && js(i) == js(i+1);
1292 
1293  if (dups)
1294  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1295  if (u.numel () != m || u.columns () != nj)
1296  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1297  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
1298  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1299 
1300  if (nj > 0)
1301  {
1302  octave_idx_type kmax = std::min (k + nj, m);
1303  if (k < m)
1304  {
1305  q.resize (m, kmax);
1306  r.resize (kmax, n + nj);
1307  }
1308  else
1309  {
1310  r.resize (k, n + nj);
1311  }
1312 
1313  OCTAVE_LOCAL_BUFFER (double, rw, kmax);
1314  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
1315  {
1316  octave_idx_type ii = i;
1317  ComplexColumnVector utmp = u.column (jsi(i));
1318  F77_XFCN (zqrinc, ZQRINC, (m, n + ii, std::min (kmax, k + ii),
1319  F77_DBLE_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1320  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (), js(ii) + 1,
1321  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), rw));
1322  }
1323  }
1324  }
1325 
1326  template <>
1327  void
1329  {
1330  octave_idx_type m = q.rows ();
1331  octave_idx_type k = r.rows ();
1332  octave_idx_type n = r.columns ();
1333 
1334  if (j < 0 || j > n-1)
1335  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1336 
1337  OCTAVE_LOCAL_BUFFER (double, rw, k);
1338  F77_XFCN (zqrdec, ZQRDEC, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1339  q.rows (),
1340  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (), j + 1, rw));
1341 
1342  if (k < m)
1343  {
1344  q.resize (m, k-1);
1345  r.resize (k-1, n-1);
1346  }
1347  else
1348  {
1349  r.resize (k, n-1);
1350  }
1351  }
1352 
1353  template <>
1354  void
1356  {
1357  octave_idx_type m = q.rows ();
1358  octave_idx_type n = r.columns ();
1359  octave_idx_type k = q.columns ();
1360 
1362  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
1363  octave_idx_type nj = js.numel ();
1364  bool dups = false;
1365  for (octave_idx_type i = 0; i < nj - 1; i++)
1366  dups = dups && js(i) == js(i+1);
1367 
1368  if (dups)
1369  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1370  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
1371  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1372 
1373  if (nj > 0)
1374  {
1375  OCTAVE_LOCAL_BUFFER (double, rw, k);
1376  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
1377  {
1378  octave_idx_type ii = i;
1379  F77_XFCN (zqrdec, ZQRDEC, (m, n - ii, k == m ? k : k - ii,
1380  F77_DBLE_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1381  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (),
1382  js(ii) + 1, rw));
1383  }
1384  if (k < m)
1385  {
1386  q.resize (m, k - nj);
1387  r.resize (k - nj, n - nj);
1388  }
1389  else
1390  {
1391  r.resize (k, n - nj);
1392  }
1393 
1394  }
1395  }
1396 
1397  template <>
1398  void
1400  {
1401  octave_idx_type m = r.rows ();
1402  octave_idx_type n = r.columns ();
1403  octave_idx_type k = std::min (m, n);
1404 
1405  if (! q.is_square () || u.numel () != n)
1406  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1407  if (j < 0 || j > m)
1408  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1409 
1410  q.resize (m + 1, m + 1);
1411  r.resize (m + 1, n);
1412  ComplexRowVector utmp = u;
1413  OCTAVE_LOCAL_BUFFER (double, rw, k);
1414  F77_XFCN (zqrinr, ZQRINR, (m, n, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1415  q.rows (),
1416  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (),
1417  j + 1, F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw));
1418 
1419  }
1420 
1421  template <>
1422  void
1424  {
1425  octave_idx_type m = r.rows ();
1426  octave_idx_type n = r.columns ();
1427 
1428  if (! q.is_square ())
1429  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
1430  if (j < 0 || j > m-1)
1431  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1432 
1434  OCTAVE_LOCAL_BUFFER (double, rw, m);
1435  F77_XFCN (zqrder, ZQRDER, (m, n, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1436  q.rows (),
1437  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (), j + 1,
1438  F77_DBLE_CMPLX_ARG (w), rw));
1439 
1440  q.resize (m - 1, m - 1);
1441  r.resize (m - 1, n);
1442  }
1443 
1444  template <>
1445  void
1447  {
1448  octave_idx_type m = q.rows ();
1449  octave_idx_type k = r.rows ();
1450  octave_idx_type n = r.columns ();
1451 
1452  if (i < 0 || i > n-1 || j < 0 || j > n-1)
1453  (*current_liboctave_error_handler) ("qrshift: index out of range");
1454 
1456  OCTAVE_LOCAL_BUFFER (double, rw, k);
1457  F77_XFCN (zqrshc, ZQRSHC, (m, n, k,
1458  F77_DBLE_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1459  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), r.rows (),
1460  i + 1, j + 1, F77_DBLE_CMPLX_ARG (w), rw));
1461  }
1462 
1463 #endif
1464 
1465  template <>
1466  void
1468  FloatComplex *tau, type qr_type)
1469  {
1470  octave_idx_type m = afact.rows ();
1471  octave_idx_type min_mn = std::min (m, n);
1472  octave_idx_type info;
1473 
1474  if (qr_type == qr<FloatComplexMatrix>::raw)
1475  {
1476  for (octave_idx_type j = 0; j < min_mn; j++)
1477  {
1478  octave_idx_type limit = j < min_mn - 1 ? j : min_mn - 1;
1479  for (octave_idx_type i = limit + 1; i < m; i++)
1480  afact.elem (i, j) *= tau[j];
1481  }
1482 
1483  r = afact;
1484  }
1485  else
1486  {
1487  // Attempt to minimize copying.
1488  if (m >= n)
1489  {
1490  // afact will become q.
1491  q = afact;
1493  r = FloatComplexMatrix (k, n);
1494  for (octave_idx_type j = 0; j < n; j++)
1495  {
1496  octave_idx_type i = 0;
1497  for (; i <= j; i++)
1498  r.xelem (i, j) = afact.xelem (i, j);
1499  for (; i < k; i++)
1500  r.xelem (i, j) = 0;
1501  }
1502  afact = FloatComplexMatrix (); // optimize memory
1503  }
1504  else
1505  {
1506  // afact will become r.
1507  q = FloatComplexMatrix (m, m);
1508  for (octave_idx_type j = 0; j < m; j++)
1509  for (octave_idx_type i = j + 1; i < m; i++)
1510  {
1511  q.xelem (i, j) = afact.xelem (i, j);
1512  afact.xelem (i, j) = 0;
1513  }
1514  r = afact;
1515  }
1516 
1517  if (m > 0)
1518  {
1519  octave_idx_type k = q.columns ();
1520  // workspace query.
1521  FloatComplex clwork;
1522  F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, F77_CMPLX_ARG (q.fortran_vec ()), m,
1523  F77_CMPLX_ARG (tau),
1524  F77_CMPLX_ARG (&clwork), -1, info));
1525 
1526  // allocate buffer and do the job.
1527  octave_idx_type lwork = clwork.real ();
1528  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
1529  OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
1530  F77_XFCN (cungqr, CUNGQR, (m, k, min_mn, F77_CMPLX_ARG (q.fortran_vec ()), m,
1531  F77_CMPLX_ARG (tau),
1532  F77_CMPLX_ARG (work), lwork, info));
1533  }
1534  }
1535  }
1536 
1537  template <>
1538  void
1540  {
1541  octave_idx_type m = a.rows ();
1542  octave_idx_type n = a.cols ();
1543 
1544  octave_idx_type min_mn = m < n ? m : n;
1545  OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
1546 
1547  octave_idx_type info = 0;
1548 
1549  FloatComplexMatrix afact = a;
1550  if (m > n && qr_type == qr<FloatComplexMatrix>::std)
1551  afact.resize (m, m);
1552 
1553  if (m > 0)
1554  {
1555  // workspace query.
1556  FloatComplex clwork;
1557  F77_XFCN (cgeqrf, CGEQRF, (m, n, F77_CMPLX_ARG (afact.fortran_vec ()), m,
1558  F77_CMPLX_ARG (tau),
1559  F77_CMPLX_ARG (&clwork), -1, info));
1560 
1561  // allocate buffer and do the job.
1562  octave_idx_type lwork = clwork.real ();
1563  lwork = std::max (lwork, static_cast<octave_idx_type> (1));
1564  OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
1565  F77_XFCN (cgeqrf, CGEQRF, (m, n, F77_CMPLX_ARG (afact.fortran_vec ()), m,
1566  F77_CMPLX_ARG (tau),
1567  F77_CMPLX_ARG (work), lwork, info));
1568  }
1569 
1570  form (n, afact, tau, qr_type);
1571  }
1572 
1573 #if defined (HAVE_QRUPDATE)
1574 
1575  template <>
1576  void
1578  const FloatComplexColumnVector& v)
1579  {
1580  octave_idx_type m = q.rows ();
1581  octave_idx_type n = r.columns ();
1582  octave_idx_type k = q.columns ();
1583 
1584  if (u.numel () != m || v.numel () != n)
1585  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1586 
1587  FloatComplexColumnVector utmp = u;
1588  FloatComplexColumnVector vtmp = v;
1590  OCTAVE_LOCAL_BUFFER (float, rw, k);
1591  F77_XFCN (cqr1up, CQR1UP, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()),
1592  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
1593  F77_CMPLX_ARG (utmp.fortran_vec ()), F77_CMPLX_ARG (vtmp.fortran_vec ()),
1594  F77_CMPLX_ARG (w), rw));
1595  }
1596 
1597  template <>
1598  void
1600  const FloatComplexMatrix& v)
1601  {
1602  octave_idx_type m = q.rows ();
1603  octave_idx_type n = r.columns ();
1604  octave_idx_type k = q.columns ();
1605 
1606  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
1607  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1608 
1610  OCTAVE_LOCAL_BUFFER (float, rw, k);
1611  for (volatile octave_idx_type i = 0; i < u.cols (); i++)
1612  {
1613  FloatComplexColumnVector utmp = u.column (i);
1614  FloatComplexColumnVector vtmp = v.column (i);
1615  F77_XFCN (cqr1up, CQR1UP, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()),
1616  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
1617  F77_CMPLX_ARG (utmp.fortran_vec ()), F77_CMPLX_ARG (vtmp.fortran_vec ()),
1618  F77_CMPLX_ARG (w), rw));
1619  }
1620  }
1621 
1622  template <>
1623  void
1625  octave_idx_type j)
1626  {
1627  octave_idx_type m = q.rows ();
1628  octave_idx_type n = r.columns ();
1629  octave_idx_type k = q.columns ();
1630 
1631  if (u.numel () != m)
1632  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1633  if (j < 0 || j > n)
1634  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1635 
1636  if (k < m)
1637  {
1638  q.resize (m, k+1);
1639  r.resize (k+1, n+1);
1640  }
1641  else
1642  {
1643  r.resize (k, n+1);
1644  }
1645 
1646  FloatComplexColumnVector utmp = u;
1647  OCTAVE_LOCAL_BUFFER (float, rw, k);
1648  F77_XFCN (cqrinc, CQRINC, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1649  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (), j + 1,
1650  F77_CONST_CMPLX_ARG (utmp.data ()), rw));
1651  }
1652 
1653  template <>
1654  void
1656  const Array<octave_idx_type>& j)
1657  {
1658  octave_idx_type m = q.rows ();
1659  octave_idx_type n = r.columns ();
1660  octave_idx_type k = q.columns ();
1661 
1663  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
1664  octave_idx_type nj = js.numel ();
1665  bool dups = false;
1666  for (octave_idx_type i = 0; i < nj - 1; i++)
1667  dups = dups && js(i) == js(i+1);
1668 
1669  if (dups)
1670  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1671  if (u.numel () != m || u.columns () != nj)
1672  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1673  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
1674  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1675 
1676  if (nj > 0)
1677  {
1678  octave_idx_type kmax = std::min (k + nj, m);
1679  if (k < m)
1680  {
1681  q.resize (m, kmax);
1682  r.resize (kmax, n + nj);
1683  }
1684  else
1685  {
1686  r.resize (k, n + nj);
1687  }
1688 
1689  OCTAVE_LOCAL_BUFFER (float, rw, kmax);
1690  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
1691  {
1692  octave_idx_type ii = i;
1693  F77_XFCN (cqrinc, CQRINC, (m, n + ii, std::min (kmax, k + ii),
1694  F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1695  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (), js(ii) + 1,
1696  F77_CONST_CMPLX_ARG (u.column (jsi(i)).data ()), rw));
1697  }
1698  }
1699  }
1700 
1701  template <>
1702  void
1704  {
1705  octave_idx_type m = q.rows ();
1706  octave_idx_type k = r.rows ();
1707  octave_idx_type n = r.columns ();
1708 
1709  if (j < 0 || j > n-1)
1710  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1711 
1712  OCTAVE_LOCAL_BUFFER (float, rw, k);
1713  F77_XFCN (cqrdec, CQRDEC, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1714  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (), j + 1, rw));
1715 
1716  if (k < m)
1717  {
1718  q.resize (m, k-1);
1719  r.resize (k-1, n-1);
1720  }
1721  else
1722  {
1723  r.resize (k, n-1);
1724  }
1725  }
1726 
1727  template <>
1728  void
1730  {
1731  octave_idx_type m = q.rows ();
1732  octave_idx_type n = r.columns ();
1733  octave_idx_type k = q.columns ();
1734 
1736  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
1737  octave_idx_type nj = js.numel ();
1738  bool dups = false;
1739  for (octave_idx_type i = 0; i < nj - 1; i++)
1740  dups = dups && js(i) == js(i+1);
1741 
1742  if (dups)
1743  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1744  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
1745  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1746 
1747  if (nj > 0)
1748  {
1749  OCTAVE_LOCAL_BUFFER (float, rw, k);
1750  for (volatile octave_idx_type i = 0; i < js.numel (); i++)
1751  {
1752  octave_idx_type ii = i;
1753  F77_XFCN (cqrdec, CQRDEC, (m, n - ii, k == m ? k : k - ii,
1754  F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1755  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (),
1756  js(ii) + 1, rw));
1757  }
1758  if (k < m)
1759  {
1760  q.resize (m, k - nj);
1761  r.resize (k - nj, n - nj);
1762  }
1763  else
1764  {
1765  r.resize (k, n - nj);
1766  }
1767 
1768  }
1769  }
1770 
1771  template <>
1772  void
1774  octave_idx_type j)
1775  {
1776  octave_idx_type m = r.rows ();
1777  octave_idx_type n = r.columns ();
1778  octave_idx_type k = std::min (m, n);
1779 
1780  if (! q.is_square () || u.numel () != n)
1781  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1782  if (j < 0 || j > m)
1783  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1784 
1785  q.resize (m + 1, m + 1);
1786  r.resize (m + 1, n);
1787  FloatComplexRowVector utmp = u;
1788  OCTAVE_LOCAL_BUFFER (float, rw, k);
1789  F77_XFCN (cqrinr, CQRINR, (m, n, F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1790  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (),
1791  j + 1, F77_CMPLX_ARG (utmp.fortran_vec ()), rw));
1792 
1793  }
1794 
1795  template <>
1796  void
1798  {
1799  octave_idx_type m = r.rows ();
1800  octave_idx_type n = r.columns ();
1801 
1802  if (! q.is_square ())
1803  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
1804  if (j < 0 || j > m-1)
1805  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1806 
1808  OCTAVE_LOCAL_BUFFER (float, rw, m);
1809  F77_XFCN (cqrder, CQRDER, (m, n, F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1810  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (), j + 1,
1811  F77_CMPLX_ARG (w), rw));
1812 
1813  q.resize (m - 1, m - 1);
1814  r.resize (m - 1, n);
1815  }
1816 
1817  template <>
1818  void
1820  {
1821  octave_idx_type m = q.rows ();
1822  octave_idx_type k = r.rows ();
1823  octave_idx_type n = r.columns ();
1824 
1825  if (i < 0 || i > n-1 || j < 0 || j > n-1)
1826  (*current_liboctave_error_handler) ("qrshift: index out of range");
1827 
1829  OCTAVE_LOCAL_BUFFER (float, rw, k);
1830  F77_XFCN (cqrshc, CQRSHC, (m, n, k,
1831  F77_CMPLX_ARG (q.fortran_vec ()), q.rows (),
1832  F77_CMPLX_ARG (r.fortran_vec ()), r.rows (),
1833  i + 1, j + 1, F77_CMPLX_ARG (w), rw));
1834  }
1835 
1836 #endif
1837 
1838  // Instantiations we need.
1839 
1840  template class qr<Matrix>;
1841 
1842  template class qr<FloatMatrix>;
1843 
1844  template class qr<ComplexMatrix>;
1845 
1846  template class qr<FloatComplexMatrix>;
1847  }
1848 }
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition: fCMatrix.h:189
qr(void)
Definition: qr.h:53
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
Definition: fMatrix.h:149
Octave interface to the compression and uncompression libraries.
Definition: aepbalance.cc:47
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
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:145
static const idx_vector colon
Definition: idx-vector.h:482
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:363
OCTAVE_EXPORT octave_value_list while another program execution is suspended until the graphics object the function returns immediately In the second form
Definition: graphics.cc:11575
void delete_col(octave_idx_type j)
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition: lo-error.c:38
#define F77_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:345
for large enough k
Definition: lu.cc:606
bool regular(void) const
Definition: qr.cc:85
void shift_cols(octave_idx_type i, octave_idx_type j)
void insert_row(const RV_T &u, octave_idx_type j)
T & elem(octave_idx_type n)
Definition: Array.h:482
u
Definition: lu.cc:138
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:52
void warn_qrupdate_once(void)
octave_idx_type rows(void) const
Definition: Array.h:401
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
void insert_col(const CV_T &u, octave_idx_type j)
void update(const CV_T &u, const CV_T &v)
void init(const T &a, type qr_type)
octave_idx_type columns(void) const
Definition: ov.h:491
FloatComplexColumnVector column(octave_idx_type i) const
Definition: fCMatrix.cc:714
Array< T > sort(int dim=0, sortmode mode=ASCENDING) const
Definition: Array.cc:1775
octave_idx_type index(const T *src, octave_idx_type n, T *dest) const
Definition: idx-vector.h:611
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
std::complex< double > w(std::complex< double > z, double relerr=0)
const T * data(void) const
Definition: Array.h:582
octave_value & assign(assign_op op, const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
Definition: ov.cc:1556
octave_value retval
Definition: data.cc:6294
Array< T > column(octave_idx_type k) const
Extract column: A(:,k+1).
Definition: Array.cc:269
idx type
Definition: ov.cc:3129
Definition: dMatrix.h:37
void form(octave_idx_type n, T &afact, ELT_T *tau, type qr_type)
FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:429
void delete_row(octave_idx_type j)
T & xelem(octave_idx_type n)
Definition: Array.h:455
type get_type(void) const
Definition: qr.cc:69
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:339
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:348
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:184
Definition: mx-defs.h:77
octave_idx_type rows(void) const
Definition: oct-map.h:375
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
p
Definition: lu.cc:138
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:423
T::element_type ELT_T
Definition: qr.h:42
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:342
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::complex< double > Complex
Definition: oct-cmplx.h:31
const T * fortran_vec(void) const
Definition: Array.h:584
octave_idx_type cols(void) const
Definition: Array.h:409
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:711
octave_idx_type columns(void) const
Definition: Array.h:410
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
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205