GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
qr.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1994-2018 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
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13 
14 Octave is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License 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 <https://www.gnu.org/licenses/>.
22 
23 */
24 
25 #if defined (HAVE_CONFIG_H)
26 # include "config.h"
27 #endif
28 
29 #include <algorithm>
30 
31 #include "Array.h"
32 #include "CColVector.h"
33 #include "CMatrix.h"
34 #include "CRowVector.h"
35 #include "dColVector.h"
36 #include "dMatrix.h"
37 #include "dRowVector.h"
38 #include "fCColVector.h"
39 #include "fCMatrix.h"
40 #include "fCRowVector.h"
41 #include "fColVector.h"
42 #include "fMatrix.h"
43 #include "fRowVector.h"
44 #include "lo-error.h"
45 #include "lo-lapack-proto.h"
46 #include "lo-qrupdate-proto.h"
47 #include "oct-cmplx.h"
48 #include "oct-locbuf.h"
49 #include "oct-sort.h"
50 #include "qr.h"
51 
52 namespace octave
53 {
54  namespace math
55  {
56  template <typename T>
57  qr<T>::qr (const T& q_arg, const T& r_arg)
58  : q (q_arg), r (r_arg)
59  {
60  octave_idx_type q_nr = q.rows ();
61  octave_idx_type q_nc = q.cols ();
62 
63  octave_idx_type r_nr = r.rows ();
64  octave_idx_type r_nc = r.cols ();
65 
66  if (! (q_nc == r_nr && (q_nr == q_nc || (q_nr > q_nc && r_nr == r_nc))))
67  (*current_liboctave_error_handler) ("QR dimensions mismatch");
68  }
69 
70  template <typename T>
71  typename qr<T>::type
72  qr<T>::get_type (void) const
73  {
74  type retval;
75 
76  if (! q.isempty () && q.issquare ())
78  else if (q.rows () > q.cols () && r.issquare ())
80  else
82 
83  return retval;
84  }
85 
86  template <typename T>
87  bool
88  qr<T>::regular (void) const
89  {
90  bool retval = true;
91 
92  octave_idx_type k = std::min (r.rows (), r.cols ());
93 
94  for (octave_idx_type i = 0; i < k; i++)
95  {
96  if (r(i, i) == ELT_T ())
97  {
98  retval = false;
99  break;
100  }
101  }
102 
103  return retval;
104  }
105 
106 #if ! defined (HAVE_QRUPDATE)
107 
108  // Replacement update methods.
109 
110  void
111  warn_qrupdate_once (void)
112  {
113  static bool warned = false;
114 
115  if (! warned)
116  {
117  (*current_liboctave_warning_with_id_handler)
118  ("Octave:missing-dependency",
119  "In this version of Octave, QR & Cholesky updating routines "
120  "simply update the matrix and recalculate factorizations. "
121  "To use fast algorithms, link Octave with the qrupdate library. "
122  "See <http://sourceforge.net/projects/qrupdate>.");
123 
124  warned = true;
125  }
126  }
127 
128  template <typename T>
129  void
130  qr<T>::update (const CV_T& u, const CV_T& v)
131  {
133 
134  octave_idx_type m = q.rows ();
135  octave_idx_type n = r.cols ();
136 
137  if (u.numel () != m || v.numel () != n)
138  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
139 
140  init (q*r + T (u) * T (v).hermitian (), get_type ());
141  }
142 
143  template <typename T>
144  void
145  qr<T>::update (const T& u, const T& v)
146  {
148 
149  octave_idx_type m = q.rows ();
150  octave_idx_type n = r.cols ();
151 
152  if (u.rows () != m || v.rows () != n || u.cols () != v.cols ())
153  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
154 
155  init (q*r + u * v.hermitian (), get_type ());
156  }
157 
158  template <typename T, typename CV_T>
159  static
160  T
161  insert_col (const T& a, octave_idx_type i, const CV_T& x)
162  {
163  T retval (a.rows (), a.cols () + 1);
165  a.index (idx_vector::colon, idx_vector (0, i)));
168  a.index (idx_vector::colon, idx_vector (i, a.cols ())));
169  return retval;
170  }
171 
172  template <typename T, typename RV_T>
173  static
174  T
175  insert_row (const T& a, octave_idx_type i, const RV_T& x)
176  {
177  T retval (a.rows () + 1, a.cols ());
179  a.index (idx_vector (0, i), idx_vector::colon));
182  a.index (idx_vector (i, a.rows ()), idx_vector::colon));
183  return retval;
184  }
185 
186  template <typename T>
187  static
188  T
189  delete_col (const T& a, octave_idx_type i)
190  {
191  T retval = a;
192  retval.delete_elements (1, idx_vector (i));
193  return retval;
194  }
195 
196  template <typename T>
197  static
198  T
199  delete_row (const T& a, octave_idx_type i)
200  {
201  T retval = a;
202  retval.delete_elements (0, idx_vector (i));
203  return retval;
204  }
205 
206  template <typename T>
207  static
208  T
209  shift_cols (const T& a, octave_idx_type i, octave_idx_type j)
210  {
211  octave_idx_type n = a.cols ();
213  for (octave_idx_type k = 0; k < n; k++) p(k) = k;
214  if (i < j)
215  {
216  for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
217  p(j) = i;
218  }
219  else if (j < i)
220  {
221  p(j) = i;
222  for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
223  }
224 
225  return a.index (idx_vector::colon, idx_vector (p));
226  }
227 
228  template <typename T>
229  void
230  qr<T>::insert_col (const CV_T& u, octave_idx_type j)
231  {
233 
234  octave_idx_type m = q.rows ();
235  octave_idx_type n = r.cols ();
236 
237  if (u.numel () != m)
238  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
239 
240  if (j < 0 || j > n)
241  (*current_liboctave_error_handler) ("qrinsert: index out of range");
242 
243  init (math::insert_col (q*r, j, u), get_type ());
244  }
245 
246  template <typename T>
247  void
248  qr<T>::insert_col (const T& u, const Array<octave_idx_type>& j)
249  {
251 
252  octave_idx_type m = q.rows ();
253  octave_idx_type n = r.cols ();
254 
256  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
257  octave_idx_type nj = js.numel ();
258  bool dups = false;
259  for (octave_idx_type i = 0; i < nj - 1; i++)
260  dups = dups && js(i) == js(i+1);
261 
262  if (dups)
263  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
264 
265  if (u.numel () != m || u.cols () != nj)
266  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
267 
268  if (nj > 0 && (js(0) < 0 || js(nj-1) > n))
269  (*current_liboctave_error_handler) ("qrinsert: index out of range");
270 
271  if (nj > 0)
272  {
273  T a = q*r;
274  for (octave_idx_type i = 0; i < nj; i++)
275  a = math::insert_col (a, js(i), u.column (i));
276 
277  init (a, get_type ());
278  }
279  }
280 
281  template <typename T>
282  void
284  {
286 
287  octave_idx_type n = r.cols ();
288 
289  if (j < 0 || j > n-1)
290  (*current_liboctave_error_handler) ("qrdelete: index out of range");
291 
292  init (math::delete_col (q*r, j), get_type ());
293  }
294 
295  template <typename T>
296  void
298  {
300 
301  octave_idx_type n = r.cols ();
302 
304  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
305  octave_idx_type nj = js.numel ();
306  bool dups = false;
307  for (octave_idx_type i = 0; i < nj - 1; i++)
308  dups = dups && js(i) == js(i+1);
309 
310  if (dups)
311  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
312 
313  if (nj > 0 && (js(0) > n-1 || js(nj-1) < 0))
314  (*current_liboctave_error_handler) ("qrinsert: index out of range");
315 
316  if (nj > 0)
317  {
318  T a = q*r;
319  for (octave_idx_type i = 0; i < nj; i++)
320  a = math::delete_col (a, js(i));
321 
322  init (a, get_type ());
323  }
324  }
325 
326  template <typename T>
327  void
328  qr<T>::insert_row (const RV_T& u, octave_idx_type j)
329  {
331 
332  octave_idx_type m = r.rows ();
333  octave_idx_type n = r.cols ();
334 
335  if (! q.issquare () || u.numel () != n)
336  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
337 
338  if (j < 0 || j > m)
339  (*current_liboctave_error_handler) ("qrinsert: index out of range");
340 
341  init (math::insert_row (q*r, j, u), get_type ());
342  }
343 
344  template <typename T>
345  void
347  {
349 
350  octave_idx_type m = r.rows ();
351 
352  if (! q.issquare ())
353  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
354 
355  if (j < 0 || j > m-1)
356  (*current_liboctave_error_handler) ("qrdelete: index out of range");
357 
358  init (math::delete_row (q*r, j), get_type ());
359  }
360 
361  template <typename T>
362  void
364  {
366 
367  octave_idx_type n = r.cols ();
368 
369  if (i < 0 || i > n-1 || j < 0 || j > n-1)
370  (*current_liboctave_error_handler) ("qrshift: index out of range");
371 
372  init (math::shift_cols (q*r, i, j), get_type ());
373  }
374 
375 #endif
376 
377  // Specializations.
378 
379  template <>
380  void
381  qr<Matrix>::form (octave_idx_type n_arg, Matrix& afact, double *tau,
382  type qr_type)
383  {
384  F77_INT n = to_f77_int (n_arg);
385  F77_INT m = to_f77_int (afact.rows ());
386  F77_INT min_mn = std::min (m, n);
387  F77_INT info;
388 
389  if (qr_type == qr<Matrix>::raw)
390  {
391  for (F77_INT j = 0; j < min_mn; j++)
392  {
393  F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
394  for (F77_INT i = limit + 1; i < m; i++)
395  afact.elem (i, j) *= tau[j];
396  }
397 
398  r = afact;
399  }
400  else
401  {
402  // Attempt to minimize copying.
403  if (m >= n)
404  {
405  // afact will become q.
406  q = afact;
407  F77_INT k = (qr_type == qr<Matrix>::economy ? n : m);
408  r = Matrix (k, n);
409  for (F77_INT j = 0; j < n; j++)
410  {
411  F77_INT i = 0;
412  for (; i <= j; i++)
413  r.xelem (i, j) = afact.xelem (i, j);
414  for (; i < k; i++)
415  r.xelem (i, j) = 0;
416  }
417  afact = Matrix (); // optimize memory
418  }
419  else
420  {
421  // afact will become r.
422  q = Matrix (m, m);
423  for (F77_INT j = 0; j < m; j++)
424  for (F77_INT i = j + 1; i < m; i++)
425  {
426  q.xelem (i, j) = afact.xelem (i, j);
427  afact.xelem (i, j) = 0;
428  }
429  r = afact;
430  }
431 
432  if (m > 0)
433  {
434  F77_INT k = to_f77_int (q.cols ());
435  // workspace query.
436  double rlwork;
437  F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m,
438  tau, &rlwork, -1, info));
439 
440  // allocate buffer and do the job.
441  F77_INT lwork = static_cast<F77_INT> (rlwork);
442  lwork = std::max (lwork, static_cast<F77_INT> (1));
443  OCTAVE_LOCAL_BUFFER (double, work, lwork);
444  F77_XFCN (dorgqr, DORGQR, (m, k, min_mn, q.fortran_vec (), m,
445  tau, work, lwork, info));
446  }
447  }
448  }
449 
450  template <>
451  void
453  {
454  F77_INT m = to_f77_int (a.rows ());
455  F77_INT n = to_f77_int (a.cols ());
456 
457  F77_INT min_mn = (m < n ? m : n);
458  OCTAVE_LOCAL_BUFFER (double, tau, min_mn);
459 
460  F77_INT info = 0;
461 
462  Matrix afact = a;
463  if (m > n && qr_type == qr<Matrix>::std)
464  afact.resize (m, m);
465 
466  if (m > 0)
467  {
468  // workspace query.
469  double rlwork;
470  F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
471  &rlwork, -1, info));
472 
473  // allocate buffer and do the job.
474  F77_INT lwork = static_cast<F77_INT> (rlwork);
475  lwork = std::max (lwork, static_cast<F77_INT> (1));
476  OCTAVE_LOCAL_BUFFER (double, work, lwork);
477  F77_XFCN (dgeqrf, DGEQRF, (m, n, afact.fortran_vec (), m, tau,
478  work, lwork, info));
479  }
480 
481  form (n, afact, tau, qr_type);
482  }
483 
484 #if defined (HAVE_QRUPDATE)
485 
486  template <>
487  void
489  {
490  F77_INT m = to_f77_int (q.rows ());
491  F77_INT n = to_f77_int (r.cols ());
492  F77_INT k = to_f77_int (q.cols ());
493 
494  F77_INT u_nel = to_f77_int (u.numel ());
495  F77_INT v_nel = to_f77_int (v.numel ());
496 
497  if (u_nel != m || v_nel != n)
498  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
499 
500  ColumnVector utmp = u;
501  ColumnVector vtmp = v;
502  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
503  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m,
504  r.fortran_vec (), k, utmp.fortran_vec (),
505  vtmp.fortran_vec (), w));
506  }
507 
508  template <>
509  void
510  qr<Matrix>::update (const Matrix& u, const Matrix& v)
511  {
512  F77_INT m = to_f77_int (q.rows ());
513  F77_INT n = to_f77_int (r.cols ());
514  F77_INT k = to_f77_int (q.cols ());
515 
516  F77_INT u_rows = to_f77_int (u.rows ());
517  F77_INT u_cols = to_f77_int (u.cols ());
518 
519  F77_INT v_rows = to_f77_int (v.rows ());
520  F77_INT v_cols = to_f77_int (v.cols ());
521 
522  if (u_rows != m || v_rows != n || u_cols != v_cols)
523  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
524 
525  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
526  for (volatile F77_INT i = 0; i < u_cols; i++)
527  {
528  ColumnVector utmp = u.column (i);
529  ColumnVector vtmp = v.column (i);
530  F77_XFCN (dqr1up, DQR1UP, (m, n, k, q.fortran_vec (), m,
531  r.fortran_vec (), k, utmp.fortran_vec (),
532  vtmp.fortran_vec (), w));
533  }
534  }
535 
536  template <>
537  void
539  {
540  F77_INT j = to_f77_int (j_arg);
541 
542  F77_INT m = to_f77_int (q.rows ());
543  F77_INT n = to_f77_int (r.cols ());
544  F77_INT k = to_f77_int (q.cols ());
545 
546  F77_INT u_nel = to_f77_int (u.numel ());
547 
548  if (u_nel != m)
549  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
550 
551  if (j < 0 || j > n)
552  (*current_liboctave_error_handler) ("qrinsert: index out of range");
553 
554  if (k < m)
555  {
556  q.resize (m, k+1);
557  r.resize (k+1, n+1);
558  }
559  else
560  r.resize (k, n+1);
561 
562  F77_INT ldq = to_f77_int (q.rows ());
563  F77_INT ldr = to_f77_int (r.rows ());
564 
565  ColumnVector utmp = u;
566  OCTAVE_LOCAL_BUFFER (double, w, k);
567  F77_XFCN (dqrinc, DQRINC, (m, n, k, q.fortran_vec (), ldq,
568  r.fortran_vec (), ldr, j + 1,
569  utmp.data (), w));
570  }
571 
572  template <>
573  void
575  {
576  F77_INT m = to_f77_int (q.rows ());
577  F77_INT n = to_f77_int (r.cols ());
578  F77_INT k = to_f77_int (q.cols ());
579 
581  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
582  F77_INT nj = to_f77_int (js.numel ());
583  bool dups = false;
584  for (F77_INT i = 0; i < nj - 1; i++)
585  dups = dups && js(i) == js(i+1);
586 
587  if (dups)
588  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
589 
590  F77_INT u_nel = to_f77_int (u.numel ());
591  F77_INT u_cols = to_f77_int (u.cols ());
592 
593  if (u_nel != m || u_cols != nj)
594  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
595 
596  F77_INT js_beg = to_f77_int (js(0));
597  F77_INT js_end = to_f77_int (js(nj-1));
598 
599  if (nj > 0 && (js_beg < 0 || js_end > n))
600  (*current_liboctave_error_handler) ("qrinsert: index out of range");
601 
602  if (nj > 0)
603  {
604  F77_INT kmax = std::min (k + nj, m);
605  if (k < m)
606  {
607  q.resize (m, kmax);
608  r.resize (kmax, n + nj);
609  }
610  else
611  r.resize (k, n + nj);
612 
613  F77_INT ldq = to_f77_int (q.rows ());
614  F77_INT ldr = to_f77_int (r.rows ());
615 
616  OCTAVE_LOCAL_BUFFER (double, w, kmax);
617  for (volatile F77_INT i = 0; i < nj; i++)
618  {
619  F77_INT ii = i;
620  ColumnVector utmp = u.column (jsi(i));
621  F77_INT js_elt = to_f77_int (js(ii));
622  F77_XFCN (dqrinc, DQRINC, (m, n + ii, std::min (kmax, k + ii),
623  q.fortran_vec (), ldq,
624  r.fortran_vec (), ldr, js_elt + 1,
625  utmp.data (), w));
626  }
627  }
628  }
629 
630  template <>
631  void
633  {
634  F77_INT j = to_f77_int (j_arg);
635 
636  F77_INT m = to_f77_int (q.rows ());
637  F77_INT k = to_f77_int (r.rows ());
638  F77_INT n = to_f77_int (r.cols ());
639 
640  if (j < 0 || j > n-1)
641  (*current_liboctave_error_handler) ("qrdelete: index out of range");
642 
643  F77_INT ldq = to_f77_int (q.rows ());
644  F77_INT ldr = to_f77_int (r.rows ());
645 
646  OCTAVE_LOCAL_BUFFER (double, w, k);
647  F77_XFCN (dqrdec, DQRDEC, (m, n, k, q.fortran_vec (), ldq,
648  r.fortran_vec (), ldr, j + 1, w));
649 
650  if (k < m)
651  {
652  q.resize (m, k-1);
653  r.resize (k-1, n-1);
654  }
655  else
656  r.resize (k, n-1);
657  }
658 
659  template <>
660  void
662  {
663  F77_INT m = to_f77_int (q.rows ());
664  F77_INT n = to_f77_int (r.cols ());
665  F77_INT k = to_f77_int (q.cols ());
666 
668  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
669  F77_INT nj = to_f77_int (js.numel ());
670  bool dups = false;
671  for (F77_INT i = 0; i < nj - 1; i++)
672  dups = dups && js(i) == js(i+1);
673 
674  if (dups)
675  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
676 
677  F77_INT js_beg = to_f77_int (js(0));
678  F77_INT js_end = to_f77_int (js(nj-1));
679 
680  if (nj > 0 && (js_beg > n-1 || js_end < 0))
681  (*current_liboctave_error_handler) ("qrinsert: index out of range");
682 
683  if (nj > 0)
684  {
685  F77_INT ldq = to_f77_int (q.rows ());
686  F77_INT ldr = to_f77_int (r.rows ());
687 
688  OCTAVE_LOCAL_BUFFER (double, w, k);
689  for (volatile F77_INT i = 0; i < nj; i++)
690  {
691  F77_INT ii = i;
692  F77_INT js_elt = to_f77_int (js(ii));
693  F77_XFCN (dqrdec, DQRDEC, (m, n - ii, (k == m ? k : k - ii),
694  q.fortran_vec (), ldq,
695  r.fortran_vec (), ldr,
696  js_elt + 1, w));
697  }
698 
699  if (k < m)
700  {
701  q.resize (m, k - nj);
702  r.resize (k - nj, n - nj);
703  }
704  else
705  r.resize (k, n - nj);
706  }
707  }
708 
709  template <>
710  void
712  {
713  F77_INT j = to_f77_int (j_arg);
714 
715  F77_INT m = to_f77_int (r.rows ());
716  F77_INT n = to_f77_int (r.cols ());
717  F77_INT k = std::min (m, n);
718 
719  F77_INT u_nel = to_f77_int (u.numel ());
720 
721  if (! q.issquare () || u_nel != n)
722  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
723 
724  if (j < 0 || j > m)
725  (*current_liboctave_error_handler) ("qrinsert: index out of range");
726 
727  q.resize (m + 1, m + 1);
728  r.resize (m + 1, n);
729 
730  F77_INT ldq = to_f77_int (q.rows ());
731  F77_INT ldr = to_f77_int (r.rows ());
732 
733  RowVector utmp = u;
734  OCTAVE_LOCAL_BUFFER (double, w, k);
735  F77_XFCN (dqrinr, DQRINR, (m, n, q.fortran_vec (), ldq,
736  r.fortran_vec (), ldr,
737  j + 1, utmp.fortran_vec (), w));
738 
739  }
740 
741  template <>
742  void
744  {
745  F77_INT j = to_f77_int (j_arg);
746 
747  F77_INT m = to_f77_int (r.rows ());
748  F77_INT n = to_f77_int (r.cols ());
749 
750  if (! q.issquare ())
751  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
752 
753  if (j < 0 || j > m-1)
754  (*current_liboctave_error_handler) ("qrdelete: index out of range");
755 
756  F77_INT ldq = to_f77_int (q.rows ());
757  F77_INT ldr = to_f77_int (r.rows ());
758 
759  OCTAVE_LOCAL_BUFFER (double, w, 2*m);
760  F77_XFCN (dqrder, DQRDER, (m, n, q.fortran_vec (), ldq,
761  r.fortran_vec (), ldr, j + 1, w));
762 
763  q.resize (m - 1, m - 1);
764  r.resize (m - 1, n);
765  }
766 
767  template <>
768  void
770  {
771  F77_INT i = to_f77_int (i_arg);
772  F77_INT j = to_f77_int (j_arg);
773 
774  F77_INT m = to_f77_int (q.rows ());
775  F77_INT k = to_f77_int (r.rows ());
776  F77_INT n = to_f77_int (r.cols ());
777 
778  if (i < 0 || i > n-1 || j < 0 || j > n-1)
779  (*current_liboctave_error_handler) ("qrshift: index out of range");
780 
781  F77_INT ldq = to_f77_int (q.rows ());
782  F77_INT ldr = to_f77_int (r.rows ());
783 
784  OCTAVE_LOCAL_BUFFER (double, w, 2*k);
785  F77_XFCN (dqrshc, DQRSHC, (m, n, k,
786  q.fortran_vec (), ldq,
787  r.fortran_vec (), ldr,
788  i + 1, j + 1, w));
789  }
790 
791 #endif
792 
793  template <>
794  void
796  float *tau, type qr_type)
797  {
798  F77_INT n = to_f77_int (n_arg);
799  F77_INT m = to_f77_int (afact.rows ());
800  F77_INT min_mn = std::min (m, n);
801  F77_INT info;
802 
804  {
805  for (F77_INT j = 0; j < min_mn; j++)
806  {
807  F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
808  for (F77_INT i = limit + 1; i < m; i++)
809  afact.elem (i, j) *= tau[j];
810  }
811 
812  r = afact;
813  }
814  else
815  {
816  // Attempt to minimize copying.
817  if (m >= n)
818  {
819  // afact will become q.
820  q = afact;
821  F77_INT k = (qr_type == qr<FloatMatrix>::economy ? n : m);
822  r = FloatMatrix (k, n);
823  for (F77_INT j = 0; j < n; j++)
824  {
825  F77_INT i = 0;
826  for (; i <= j; i++)
827  r.xelem (i, j) = afact.xelem (i, j);
828  for (; i < k; i++)
829  r.xelem (i, j) = 0;
830  }
831  afact = FloatMatrix (); // optimize memory
832  }
833  else
834  {
835  // afact will become r.
836  q = FloatMatrix (m, m);
837  for (F77_INT j = 0; j < m; j++)
838  for (F77_INT i = j + 1; i < m; i++)
839  {
840  q.xelem (i, j) = afact.xelem (i, j);
841  afact.xelem (i, j) = 0;
842  }
843  r = afact;
844  }
845 
846  if (m > 0)
847  {
848  F77_INT k = to_f77_int (q.cols ());
849  // workspace query.
850  float rlwork;
851  F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m,
852  tau, &rlwork, -1, info));
853 
854  // allocate buffer and do the job.
855  F77_INT lwork = static_cast<F77_INT> (rlwork);
856  lwork = std::max (lwork, static_cast<F77_INT> (1));
857  OCTAVE_LOCAL_BUFFER (float, work, lwork);
858  F77_XFCN (sorgqr, SORGQR, (m, k, min_mn, q.fortran_vec (), m,
859  tau, work, lwork, info));
860  }
861  }
862  }
863 
864  template <>
865  void
867  {
868  F77_INT m = to_f77_int (a.rows ());
869  F77_INT n = to_f77_int (a.cols ());
870 
871  F77_INT min_mn = (m < n ? m : n);
872  OCTAVE_LOCAL_BUFFER (float, tau, min_mn);
873 
874  F77_INT info = 0;
875 
876  FloatMatrix afact = a;
877  if (m > n && qr_type == qr<FloatMatrix>::std)
878  afact.resize (m, m);
879 
880  if (m > 0)
881  {
882  // workspace query.
883  float rlwork;
884  F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
885  &rlwork, -1, info));
886 
887  // allocate buffer and do the job.
888  F77_INT lwork = static_cast<F77_INT> (rlwork);
889  lwork = std::max (lwork, static_cast<F77_INT> (1));
890  OCTAVE_LOCAL_BUFFER (float, work, lwork);
891  F77_XFCN (sgeqrf, SGEQRF, (m, n, afact.fortran_vec (), m, tau,
892  work, lwork, info));
893  }
894 
895  form (n, afact, tau, qr_type);
896  }
897 
898 #if defined (HAVE_QRUPDATE)
899 
900  template <>
901  void
903  {
904  F77_INT m = to_f77_int (q.rows ());
905  F77_INT n = to_f77_int (r.cols ());
906  F77_INT k = to_f77_int (q.cols ());
907 
908  F77_INT u_nel = to_f77_int (u.numel ());
909  F77_INT v_nel = to_f77_int (v.numel ());
910 
911  if (u_nel != m || v_nel != n)
912  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
913 
914  FloatColumnVector utmp = u;
915  FloatColumnVector vtmp = v;
916  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
917  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), m,
918  r.fortran_vec (), k, utmp.fortran_vec (),
919  vtmp.fortran_vec (), w));
920  }
921 
922  template <>
923  void
925  {
926  F77_INT m = to_f77_int (q.rows ());
927  F77_INT n = to_f77_int (r.cols ());
928  F77_INT k = to_f77_int (q.cols ());
929 
930  F77_INT u_rows = to_f77_int (u.rows ());
931  F77_INT u_cols = to_f77_int (u.cols ());
932 
933  F77_INT v_rows = to_f77_int (v.rows ());
934  F77_INT v_cols = to_f77_int (v.cols ());
935 
936  if (u_rows != m || v_rows != n || u_cols != v_cols)
937  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
938 
939  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
940  for (volatile F77_INT i = 0; i < u_cols; i++)
941  {
942  FloatColumnVector utmp = u.column (i);
943  FloatColumnVector vtmp = v.column (i);
944  F77_XFCN (sqr1up, SQR1UP, (m, n, k, q.fortran_vec (), m,
945  r.fortran_vec (), k, utmp.fortran_vec (),
946  vtmp.fortran_vec (), w));
947  }
948  }
949 
950  template <>
951  void
953  octave_idx_type j_arg)
954  {
955  F77_INT j = to_f77_int (j_arg);
956 
957  F77_INT m = to_f77_int (q.rows ());
958  F77_INT n = to_f77_int (r.cols ());
959  F77_INT k = to_f77_int (q.cols ());
960 
961  F77_INT u_nel = to_f77_int (u.numel ());
962 
963  if (u_nel != m)
964  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
965 
966  if (j < 0 || j > n)
967  (*current_liboctave_error_handler) ("qrinsert: index out of range");
968 
969  if (k < m)
970  {
971  q.resize (m, k+1);
972  r.resize (k+1, n+1);
973  }
974  else
975  r.resize (k, n+1);
976 
977  F77_INT ldq = to_f77_int (q.rows ());
978  F77_INT ldr = to_f77_int (r.rows ());
979 
980  FloatColumnVector utmp = u;
981  OCTAVE_LOCAL_BUFFER (float, w, k);
982  F77_XFCN (sqrinc, SQRINC, (m, n, k, q.fortran_vec (), ldq,
983  r.fortran_vec (), ldr, j + 1,
984  utmp.data (), w));
985  }
986 
987  template <>
988  void
990  const Array<octave_idx_type>& j)
991  {
992  F77_INT m = to_f77_int (q.rows ());
993  F77_INT n = to_f77_int (r.cols ());
994  F77_INT k = to_f77_int (q.cols ());
995 
997  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
998  F77_INT nj = to_f77_int (js.numel ());
999  bool dups = false;
1000  for (F77_INT i = 0; i < nj - 1; i++)
1001  dups = dups && js(i) == js(i+1);
1002 
1003  if (dups)
1004  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1005 
1006  F77_INT u_nel = to_f77_int (u.numel ());
1007  F77_INT u_cols = to_f77_int (u.cols ());
1008 
1009  if (u_nel != m || u_cols != nj)
1010  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1011 
1012  F77_INT js_beg = to_f77_int (js(0));
1013  F77_INT js_end = to_f77_int (js(nj-1));
1014 
1015  if (nj > 0 && (js_beg < 0 || js_end > n))
1016  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1017 
1018  if (nj > 0)
1019  {
1020  F77_INT kmax = std::min (k + nj, m);
1021  if (k < m)
1022  {
1023  q.resize (m, kmax);
1024  r.resize (kmax, n + nj);
1025  }
1026  else
1027  r.resize (k, n + nj);
1028 
1029  F77_INT ldq = to_f77_int (q.rows ());
1030  F77_INT ldr = to_f77_int (r.rows ());
1031 
1032  OCTAVE_LOCAL_BUFFER (float, w, kmax);
1033  for (volatile F77_INT i = 0; i < nj; i++)
1034  {
1035  F77_INT ii = i;
1036  FloatColumnVector utmp = u.column (jsi(i));
1037  F77_INT js_elt = to_f77_int (js(ii));
1038  F77_XFCN (sqrinc, SQRINC, (m, n + ii, std::min (kmax, k + ii),
1039  q.fortran_vec (), ldq,
1040  r.fortran_vec (), ldr, js_elt + 1,
1041  utmp.data (), w));
1042  }
1043  }
1044  }
1045 
1046  template <>
1047  void
1049  {
1050  F77_INT j = to_f77_int (j_arg);
1051 
1052  F77_INT m = to_f77_int (q.rows ());
1053  F77_INT k = to_f77_int (r.rows ());
1054  F77_INT n = to_f77_int (r.cols ());
1055 
1056  if (j < 0 || j > n-1)
1057  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1058 
1059  F77_INT ldq = to_f77_int (q.rows ());
1060  F77_INT ldr = to_f77_int (r.rows ());
1061 
1062  OCTAVE_LOCAL_BUFFER (float, w, k);
1063  F77_XFCN (sqrdec, SQRDEC, (m, n, k, q.fortran_vec (), ldq,
1064  r.fortran_vec (), ldr, j + 1, w));
1065 
1066  if (k < m)
1067  {
1068  q.resize (m, k-1);
1069  r.resize (k-1, n-1);
1070  }
1071  else
1072  r.resize (k, n-1);
1073  }
1074 
1075  template <>
1076  void
1078  {
1079  F77_INT m = to_f77_int (q.rows ());
1080  F77_INT n = to_f77_int (r.cols ());
1081  F77_INT k = to_f77_int (q.cols ());
1082 
1084  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
1085  F77_INT nj = to_f77_int (js.numel ());
1086  bool dups = false;
1087  for (F77_INT i = 0; i < nj - 1; i++)
1088  dups = dups && js(i) == js(i+1);
1089 
1090  if (dups)
1091  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1092 
1093  F77_INT js_beg = to_f77_int (js(0));
1094  F77_INT js_end = to_f77_int (js(nj-1));
1095 
1096  if (nj > 0 && (js_beg > n-1 || js_end < 0))
1097  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1098 
1099  if (nj > 0)
1100  {
1101  F77_INT ldq = to_f77_int (q.rows ());
1102  F77_INT ldr = to_f77_int (r.rows ());
1103 
1104  OCTAVE_LOCAL_BUFFER (float, w, k);
1105  for (volatile F77_INT i = 0; i < nj; i++)
1106  {
1107  F77_INT ii = i;
1108  F77_INT js_elt = to_f77_int (js(ii));
1109  F77_XFCN (sqrdec, SQRDEC, (m, n - ii, (k == m ? k : k - ii),
1110  q.fortran_vec (), ldq,
1111  r.fortran_vec (), ldr,
1112  js_elt + 1, w));
1113  }
1114 
1115  if (k < m)
1116  {
1117  q.resize (m, k - nj);
1118  r.resize (k - nj, n - nj);
1119  }
1120  else
1121  r.resize (k, n - nj);
1122  }
1123  }
1124 
1125  template <>
1126  void
1128  octave_idx_type j_arg)
1129  {
1130  F77_INT j = to_f77_int (j_arg);
1131 
1132  F77_INT m = to_f77_int (r.rows ());
1133  F77_INT n = to_f77_int (r.cols ());
1134  F77_INT k = std::min (m, n);
1135 
1136  F77_INT u_nel = to_f77_int (u.numel ());
1137 
1138  if (! q.issquare () || u_nel != n)
1139  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1140 
1141  if (j < 0 || j > m)
1142  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1143 
1144  q.resize (m + 1, m + 1);
1145  r.resize (m + 1, n);
1146 
1147  F77_INT ldq = to_f77_int (q.rows ());
1148  F77_INT ldr = to_f77_int (r.rows ());
1149 
1150  FloatRowVector utmp = u;
1151  OCTAVE_LOCAL_BUFFER (float, w, k);
1152  F77_XFCN (sqrinr, SQRINR, (m, n, q.fortran_vec (), ldq,
1153  r.fortran_vec (), ldr,
1154  j + 1, utmp.fortran_vec (), w));
1155 
1156  }
1157 
1158  template <>
1159  void
1161  {
1162  F77_INT j = to_f77_int (j_arg);
1163 
1164  F77_INT m = to_f77_int (r.rows ());
1165  F77_INT n = to_f77_int (r.cols ());
1166 
1167  if (! q.issquare ())
1168  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
1169 
1170  if (j < 0 || j > m-1)
1171  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1172 
1173  F77_INT ldq = to_f77_int (q.rows ());
1174  F77_INT ldr = to_f77_int (r.rows ());
1175 
1176  OCTAVE_LOCAL_BUFFER (float, w, 2*m);
1177  F77_XFCN (sqrder, SQRDER, (m, n, q.fortran_vec (), ldq,
1178  r.fortran_vec (), ldr, j + 1,
1179  w));
1180 
1181  q.resize (m - 1, m - 1);
1182  r.resize (m - 1, n);
1183  }
1184 
1185  template <>
1186  void
1188  {
1189  F77_INT i = to_f77_int (i_arg);
1190  F77_INT j = to_f77_int (j_arg);
1191 
1192  F77_INT m = to_f77_int (q.rows ());
1193  F77_INT k = to_f77_int (r.rows ());
1194  F77_INT n = to_f77_int (r.cols ());
1195 
1196  if (i < 0 || i > n-1 || j < 0 || j > n-1)
1197  (*current_liboctave_error_handler) ("qrshift: index out of range");
1198 
1199  F77_INT ldq = to_f77_int (q.rows ());
1200  F77_INT ldr = to_f77_int (r.rows ());
1201 
1202  OCTAVE_LOCAL_BUFFER (float, w, 2*k);
1203  F77_XFCN (sqrshc, SQRSHC, (m, n, k,
1204  q.fortran_vec (), ldq,
1205  r.fortran_vec (), ldr,
1206  i + 1, j + 1, w));
1207  }
1208 
1209 #endif
1210 
1211  template <>
1212  void
1214  Complex *tau, type qr_type)
1215  {
1216  F77_INT n = to_f77_int (n_arg);
1217  F77_INT m = to_f77_int (afact.rows ());
1218  F77_INT min_mn = std::min (m, n);
1219  F77_INT info;
1220 
1222  {
1223  for (F77_INT j = 0; j < min_mn; j++)
1224  {
1225  F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
1226  for (F77_INT i = limit + 1; i < m; i++)
1227  afact.elem (i, j) *= tau[j];
1228  }
1229 
1230  r = afact;
1231  }
1232  else
1233  {
1234  // Attempt to minimize copying.
1235  if (m >= n)
1236  {
1237  // afact will become q.
1238  q = afact;
1240  r = ComplexMatrix (k, n);
1241  for (F77_INT j = 0; j < n; j++)
1242  {
1243  F77_INT i = 0;
1244  for (; i <= j; i++)
1245  r.xelem (i, j) = afact.xelem (i, j);
1246  for (; i < k; i++)
1247  r.xelem (i, j) = 0;
1248  }
1249  afact = ComplexMatrix (); // optimize memory
1250  }
1251  else
1252  {
1253  // afact will become r.
1254  q = ComplexMatrix (m, m);
1255  for (F77_INT j = 0; j < m; j++)
1256  for (F77_INT i = j + 1; i < m; i++)
1257  {
1258  q.xelem (i, j) = afact.xelem (i, j);
1259  afact.xelem (i, j) = 0;
1260  }
1261  r = afact;
1262  }
1263 
1264  if (m > 0)
1265  {
1266  F77_INT k = to_f77_int (q.cols ());
1267  // workspace query.
1268  Complex clwork;
1269  F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn,
1270  F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1271  m, F77_DBLE_CMPLX_ARG (tau),
1272  F77_DBLE_CMPLX_ARG (&clwork), -1,
1273  info));
1274 
1275  // allocate buffer and do the job.
1276  F77_INT lwork = static_cast<F77_INT> (clwork.real ());
1277  lwork = std::max (lwork, static_cast<F77_INT> (1));
1278  OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
1279  F77_XFCN (zungqr, ZUNGQR, (m, k, min_mn,
1280  F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1281  m, F77_DBLE_CMPLX_ARG (tau),
1282  F77_DBLE_CMPLX_ARG (work), lwork,
1283  info));
1284  }
1285  }
1286  }
1287 
1288  template <>
1289  void
1291  {
1292  F77_INT m = to_f77_int (a.rows ());
1293  F77_INT n = to_f77_int (a.cols ());
1294 
1295  F77_INT min_mn = (m < n ? m : n);
1296  OCTAVE_LOCAL_BUFFER (Complex, tau, min_mn);
1297 
1298  F77_INT info = 0;
1299 
1300  ComplexMatrix afact = a;
1301  if (m > n && qr_type == qr<ComplexMatrix>::std)
1302  afact.resize (m, m);
1303 
1304  if (m > 0)
1305  {
1306  // workspace query.
1307  Complex clwork;
1308  F77_XFCN (zgeqrf, ZGEQRF, (m, n,
1309  F77_DBLE_CMPLX_ARG (afact.fortran_vec ()),
1310  m, F77_DBLE_CMPLX_ARG (tau),
1311  F77_DBLE_CMPLX_ARG (&clwork), -1, info));
1312 
1313  // allocate buffer and do the job.
1314  F77_INT lwork = static_cast<F77_INT> (clwork.real ());
1315  lwork = std::max (lwork, static_cast<F77_INT> (1));
1316  OCTAVE_LOCAL_BUFFER (Complex, work, lwork);
1317  F77_XFCN (zgeqrf, ZGEQRF, (m, n,
1318  F77_DBLE_CMPLX_ARG (afact.fortran_vec ()),
1319  m, F77_DBLE_CMPLX_ARG (tau),
1320  F77_DBLE_CMPLX_ARG (work), lwork, info));
1321  }
1322 
1323  form (n, afact, tau, qr_type);
1324  }
1325 
1326 #if defined (HAVE_QRUPDATE)
1327 
1328  template <>
1329  void
1331  const ComplexColumnVector& v)
1332  {
1333  F77_INT m = to_f77_int (q.rows ());
1334  F77_INT n = to_f77_int (r.cols ());
1335  F77_INT k = to_f77_int (q.cols ());
1336 
1337  F77_INT u_nel = to_f77_int (u.numel ());
1338  F77_INT v_nel = to_f77_int (v.numel ());
1339 
1340  if (u_nel != m || v_nel != n)
1341  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1342 
1343  ComplexColumnVector utmp = u;
1344  ComplexColumnVector vtmp = v;
1346  OCTAVE_LOCAL_BUFFER (double, rw, k);
1347  F77_XFCN (zqr1up, ZQR1UP, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1348  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
1349  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
1350  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()),
1351  F77_DBLE_CMPLX_ARG (w), rw));
1352  }
1353 
1354  template <>
1355  void
1357  {
1358  F77_INT m = to_f77_int (q.rows ());
1359  F77_INT n = to_f77_int (r.cols ());
1360  F77_INT k = to_f77_int (q.cols ());
1361 
1362  F77_INT u_rows = to_f77_int (u.rows ());
1363  F77_INT u_cols = to_f77_int (u.cols ());
1364 
1365  F77_INT v_rows = to_f77_int (v.rows ());
1366  F77_INT v_cols = to_f77_int (v.cols ());
1367 
1368  if (u_rows != m || v_rows != n || u_cols != v_cols)
1369  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1370 
1372  OCTAVE_LOCAL_BUFFER (double, rw, k);
1373  for (volatile F77_INT i = 0; i < u_cols; i++)
1374  {
1375  ComplexColumnVector utmp = u.column (i);
1376  ComplexColumnVector vtmp = v.column (i);
1377  F77_XFCN (zqr1up, ZQR1UP, (m, n, k,
1378  F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1379  m, F77_DBLE_CMPLX_ARG (r.fortran_vec ()), k,
1380  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
1381  F77_DBLE_CMPLX_ARG (vtmp.fortran_vec ()),
1382  F77_DBLE_CMPLX_ARG (w), rw));
1383  }
1384  }
1385 
1386  template <>
1387  void
1389  octave_idx_type j_arg)
1390  {
1391  F77_INT j = to_f77_int (j_arg);
1392 
1393  F77_INT m = to_f77_int (q.rows ());
1394  F77_INT n = to_f77_int (r.cols ());
1395  F77_INT k = to_f77_int (q.cols ());
1396 
1397  F77_INT u_nel = to_f77_int (u.numel ());
1398 
1399  if (u_nel != m)
1400  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1401 
1402  if (j < 0 || j > n)
1403  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1404 
1405  if (k < m)
1406  {
1407  q.resize (m, k+1);
1408  r.resize (k+1, n+1);
1409  }
1410  else
1411  r.resize (k, n+1);
1412 
1413  F77_INT ldq = to_f77_int (q.rows ());
1414  F77_INT ldr = to_f77_int (r.rows ());
1415 
1416  ComplexColumnVector utmp = u;
1417  OCTAVE_LOCAL_BUFFER (double, rw, k);
1418  F77_XFCN (zqrinc, ZQRINC, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1419  ldq, F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
1420  ldr, j + 1,
1421  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), rw));
1422  }
1423 
1424  template <>
1425  void
1427  const Array<octave_idx_type>& j)
1428  {
1429  F77_INT m = to_f77_int (q.rows ());
1430  F77_INT n = to_f77_int (r.cols ());
1431  F77_INT k = to_f77_int (q.cols ());
1432 
1434  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
1435  F77_INT nj = to_f77_int (js.numel ());
1436  bool dups = false;
1437  for (F77_INT i = 0; i < nj - 1; i++)
1438  dups = dups && js(i) == js(i+1);
1439 
1440  if (dups)
1441  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1442 
1443  F77_INT u_nel = to_f77_int (u.numel ());
1444  F77_INT u_cols = to_f77_int (u.cols ());
1445 
1446  if (u_nel != m || u_cols != nj)
1447  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1448 
1449  F77_INT js_beg = to_f77_int (js(0));
1450  F77_INT js_end = to_f77_int (js(nj-1));
1451 
1452  if (nj > 0 && (js_beg < 0 || js_end > n))
1453  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1454 
1455  if (nj > 0)
1456  {
1457  F77_INT kmax = std::min (k + nj, m);
1458  if (k < m)
1459  {
1460  q.resize (m, kmax);
1461  r.resize (kmax, n + nj);
1462  }
1463  else
1464  r.resize (k, n + nj);
1465 
1466  F77_INT ldq = to_f77_int (q.rows ());
1467  F77_INT ldr = to_f77_int (r.rows ());
1468 
1469  OCTAVE_LOCAL_BUFFER (double, rw, kmax);
1470  for (volatile F77_INT i = 0; i < nj; i++)
1471  {
1472  F77_INT ii = i;
1473  ComplexColumnVector utmp = u.column (jsi(i));
1474  F77_INT js_elt = to_f77_int (js(ii));
1475  F77_XFCN (zqrinc, ZQRINC, (m, n + ii, std::min (kmax, k + ii),
1476  F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1477  ldq,
1478  F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
1479  ldr, js_elt + 1,
1480  F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
1481  rw));
1482  }
1483  }
1484  }
1485 
1486  template <>
1487  void
1489  {
1490  F77_INT j = to_f77_int (j_arg);
1491 
1492  F77_INT m = to_f77_int (q.rows ());
1493  F77_INT k = to_f77_int (r.rows ());
1494  F77_INT n = to_f77_int (r.cols ());
1495 
1496  if (j < 0 || j > n-1)
1497  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1498 
1499  F77_INT ldq = to_f77_int (q.rows ());
1500  F77_INT ldr = to_f77_int (r.rows ());
1501 
1502  OCTAVE_LOCAL_BUFFER (double, rw, k);
1503  F77_XFCN (zqrdec, ZQRDEC, (m, n, k, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1504  ldq, F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
1505  ldr, j + 1, rw));
1506 
1507  if (k < m)
1508  {
1509  q.resize (m, k-1);
1510  r.resize (k-1, n-1);
1511  }
1512  else
1513  r.resize (k, n-1);
1514  }
1515 
1516  template <>
1517  void
1519  {
1520  F77_INT m = to_f77_int (q.rows ());
1521  F77_INT n = to_f77_int (r.cols ());
1522  F77_INT k = to_f77_int (q.cols ());
1523 
1525  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
1526  F77_INT nj = to_f77_int (js.numel ());
1527  bool dups = false;
1528  for (F77_INT i = 0; i < nj - 1; i++)
1529  dups = dups && js(i) == js(i+1);
1530 
1531  if (dups)
1532  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1533 
1534  F77_INT js_beg = to_f77_int (js(0));
1535  F77_INT js_end = to_f77_int (js(nj-1));
1536 
1537  if (nj > 0 && (js_beg > n-1 || js_end < 0))
1538  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1539 
1540  if (nj > 0)
1541  {
1542  F77_INT ldq = to_f77_int (q.rows ());
1543  F77_INT ldr = to_f77_int (r.rows ());
1544 
1545  OCTAVE_LOCAL_BUFFER (double, rw, k);
1546  for (volatile F77_INT i = 0; i < nj; i++)
1547  {
1548  F77_INT ii = i;
1549  F77_INT js_elt = to_f77_int (js(ii));
1550  F77_XFCN (zqrdec, ZQRDEC, (m, n - ii, (k == m ? k : k - ii),
1551  F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1552  ldq,
1553  F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
1554  ldr, js_elt + 1, rw));
1555  }
1556 
1557  if (k < m)
1558  {
1559  q.resize (m, k - nj);
1560  r.resize (k - nj, n - nj);
1561  }
1562  else
1563  r.resize (k, n - nj);
1564  }
1565  }
1566 
1567  template <>
1568  void
1570  octave_idx_type j_arg)
1571  {
1572  F77_INT j = to_f77_int (j_arg);
1573 
1574  F77_INT m = to_f77_int (r.rows ());
1575  F77_INT n = to_f77_int (r.cols ());
1576  F77_INT k = std::min (m, n);
1577 
1578  F77_INT u_nel = to_f77_int (u.numel ());
1579 
1580  if (! q.issquare () || u_nel != n)
1581  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1582 
1583  if (j < 0 || j > m)
1584  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1585 
1586  q.resize (m + 1, m + 1);
1587  r.resize (m + 1, n);
1588 
1589  F77_INT ldq = to_f77_int (q.rows ());
1590  F77_INT ldr = to_f77_int (r.rows ());
1591 
1592  ComplexRowVector utmp = u;
1593  OCTAVE_LOCAL_BUFFER (double, rw, k);
1594  F77_XFCN (zqrinr, ZQRINR, (m, n, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1595  ldq, F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
1596  ldr, j + 1,
1597  F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()), rw));
1598 
1599  }
1600 
1601  template <>
1602  void
1604  {
1605  F77_INT j = to_f77_int (j_arg);
1606 
1607  F77_INT m = to_f77_int (r.rows ());
1608  F77_INT n = to_f77_int (r.cols ());
1609 
1610  if (! q.issquare ())
1611  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
1612 
1613  if (j < 0 || j > m-1)
1614  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1615 
1616  F77_INT ldq = to_f77_int (q.rows ());
1617  F77_INT ldr = to_f77_int (r.rows ());
1618 
1620  OCTAVE_LOCAL_BUFFER (double, rw, m);
1621  F77_XFCN (zqrder, ZQRDER, (m, n, F77_DBLE_CMPLX_ARG (q.fortran_vec ()),
1622  ldq, F77_DBLE_CMPLX_ARG (r.fortran_vec ()),
1623  ldr, j + 1, F77_DBLE_CMPLX_ARG (w), rw));
1624 
1625  q.resize (m - 1, m - 1);
1626  r.resize (m - 1, n);
1627  }
1628 
1629  template <>
1630  void
1632  octave_idx_type j_arg)
1633  {
1634  F77_INT i = to_f77_int (i_arg);
1635  F77_INT j = to_f77_int (j_arg);
1636 
1637  F77_INT m = to_f77_int (q.rows ());
1638  F77_INT k = to_f77_int (r.rows ());
1639  F77_INT n = to_f77_int (r.cols ());
1640 
1641  if (i < 0 || i > n-1 || j < 0 || j > n-1)
1642  (*current_liboctave_error_handler) ("qrshift: index out of range");
1643 
1644  F77_INT ldq = to_f77_int (q.rows ());
1645  F77_INT ldr = to_f77_int (r.rows ());
1646 
1648  OCTAVE_LOCAL_BUFFER (double, rw, k);
1649  F77_XFCN (zqrshc, ZQRSHC, (m, n, k,
1650  F77_DBLE_CMPLX_ARG (q.fortran_vec ()), ldq,
1651  F77_DBLE_CMPLX_ARG (r.fortran_vec ()), ldr,
1652  i + 1, j + 1, F77_DBLE_CMPLX_ARG (w), rw));
1653  }
1654 
1655 #endif
1656 
1657  template <>
1658  void
1660  FloatComplex *tau, type qr_type)
1661  {
1662  F77_INT n = to_f77_int (n_arg);
1663  F77_INT m = to_f77_int (afact.rows ());
1664  F77_INT min_mn = std::min (m, n);
1665  F77_INT info;
1666 
1668  {
1669  for (F77_INT j = 0; j < min_mn; j++)
1670  {
1671  F77_INT limit = (j < min_mn - 1 ? j : min_mn - 1);
1672  for (F77_INT i = limit + 1; i < m; i++)
1673  afact.elem (i, j) *= tau[j];
1674  }
1675 
1676  r = afact;
1677  }
1678  else
1679  {
1680  // Attempt to minimize copying.
1681  if (m >= n)
1682  {
1683  // afact will become q.
1684  q = afact;
1686  r = FloatComplexMatrix (k, n);
1687  for (F77_INT j = 0; j < n; j++)
1688  {
1689  F77_INT i = 0;
1690  for (; i <= j; i++)
1691  r.xelem (i, j) = afact.xelem (i, j);
1692  for (; i < k; i++)
1693  r.xelem (i, j) = 0;
1694  }
1695  afact = FloatComplexMatrix (); // optimize memory
1696  }
1697  else
1698  {
1699  // afact will become r.
1700  q = FloatComplexMatrix (m, m);
1701  for (F77_INT j = 0; j < m; j++)
1702  for (F77_INT i = j + 1; i < m; i++)
1703  {
1704  q.xelem (i, j) = afact.xelem (i, j);
1705  afact.xelem (i, j) = 0;
1706  }
1707  r = afact;
1708  }
1709 
1710  if (m > 0)
1711  {
1712  F77_INT k = to_f77_int (q.cols ());
1713  // workspace query.
1714  FloatComplex clwork;
1715  F77_XFCN (cungqr, CUNGQR, (m, k, min_mn,
1716  F77_CMPLX_ARG (q.fortran_vec ()), m,
1717  F77_CMPLX_ARG (tau),
1718  F77_CMPLX_ARG (&clwork), -1, info));
1719 
1720  // allocate buffer and do the job.
1721  F77_INT lwork = static_cast<F77_INT> (clwork.real ());
1722  lwork = std::max (lwork, static_cast<F77_INT> (1));
1723  OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
1724  F77_XFCN (cungqr, CUNGQR, (m, k, min_mn,
1725  F77_CMPLX_ARG (q.fortran_vec ()), m,
1726  F77_CMPLX_ARG (tau),
1727  F77_CMPLX_ARG (work), lwork, info));
1728  }
1729  }
1730  }
1731 
1732  template <>
1733  void
1735  {
1736  F77_INT m = to_f77_int (a.rows ());
1737  F77_INT n = to_f77_int (a.cols ());
1738 
1739  F77_INT min_mn = (m < n ? m : n);
1740  OCTAVE_LOCAL_BUFFER (FloatComplex, tau, min_mn);
1741 
1742  F77_INT info = 0;
1743 
1744  FloatComplexMatrix afact = a;
1745  if (m > n && qr_type == qr<FloatComplexMatrix>::std)
1746  afact.resize (m, m);
1747 
1748  if (m > 0)
1749  {
1750  // workspace query.
1751  FloatComplex clwork;
1752  F77_XFCN (cgeqrf, CGEQRF, (m, n, F77_CMPLX_ARG (afact.fortran_vec ()),
1753  m, F77_CMPLX_ARG (tau),
1754  F77_CMPLX_ARG (&clwork), -1, info));
1755 
1756  // allocate buffer and do the job.
1757  F77_INT lwork = static_cast<F77_INT> (clwork.real ());
1758  lwork = std::max (lwork, static_cast<F77_INT> (1));
1759  OCTAVE_LOCAL_BUFFER (FloatComplex, work, lwork);
1760  F77_XFCN (cgeqrf, CGEQRF, (m, n, F77_CMPLX_ARG (afact.fortran_vec ()),
1761  m, F77_CMPLX_ARG (tau),
1762  F77_CMPLX_ARG (work), lwork, info));
1763  }
1764 
1765  form (n, afact, tau, qr_type);
1766  }
1767 
1768 #if defined (HAVE_QRUPDATE)
1769 
1770  template <>
1771  void
1773  const FloatComplexColumnVector& v)
1774  {
1775  F77_INT m = to_f77_int (q.rows ());
1776  F77_INT n = to_f77_int (r.cols ());
1777  F77_INT k = to_f77_int (q.cols ());
1778 
1779  F77_INT u_nel = to_f77_int (u.numel ());
1780  F77_INT v_nel = to_f77_int (v.numel ());
1781 
1782  if (u_nel != m || v_nel != n)
1783  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1784 
1785  FloatComplexColumnVector utmp = u;
1786  FloatComplexColumnVector vtmp = v;
1788  OCTAVE_LOCAL_BUFFER (float, rw, k);
1789  F77_XFCN (cqr1up, CQR1UP, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()),
1790  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
1791  F77_CMPLX_ARG (utmp.fortran_vec ()),
1792  F77_CMPLX_ARG (vtmp.fortran_vec ()),
1793  F77_CMPLX_ARG (w), rw));
1794  }
1795 
1796  template <>
1797  void
1799  const FloatComplexMatrix& v)
1800  {
1801  F77_INT m = to_f77_int (q.rows ());
1802  F77_INT n = to_f77_int (r.cols ());
1803  F77_INT k = to_f77_int (q.cols ());
1804 
1805  F77_INT u_rows = to_f77_int (u.rows ());
1806  F77_INT u_cols = to_f77_int (u.cols ());
1807 
1808  F77_INT v_rows = to_f77_int (v.rows ());
1809  F77_INT v_cols = to_f77_int (v.cols ());
1810 
1811  if (u_rows != m || v_rows != n || u_cols != v_cols)
1812  (*current_liboctave_error_handler) ("qrupdate: dimensions mismatch");
1813 
1815  OCTAVE_LOCAL_BUFFER (float, rw, k);
1816  for (volatile F77_INT i = 0; i < u_cols; i++)
1817  {
1818  FloatComplexColumnVector utmp = u.column (i);
1819  FloatComplexColumnVector vtmp = v.column (i);
1820  F77_XFCN (cqr1up, CQR1UP, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()),
1821  m, F77_CMPLX_ARG (r.fortran_vec ()), k,
1822  F77_CMPLX_ARG (utmp.fortran_vec ()),
1823  F77_CMPLX_ARG (vtmp.fortran_vec ()),
1824  F77_CMPLX_ARG (w), rw));
1825  }
1826  }
1827 
1828  template <>
1829  void
1831  octave_idx_type j_arg)
1832  {
1833  F77_INT j = to_f77_int (j_arg);
1834 
1835  F77_INT m = to_f77_int (q.rows ());
1836  F77_INT n = to_f77_int (r.cols ());
1837  F77_INT k = to_f77_int (q.cols ());
1838 
1839  F77_INT u_nel = to_f77_int (u.numel ());
1840 
1841  if (u_nel != m)
1842  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1843 
1844  if (j < 0 || j > n)
1845  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1846 
1847  if (k < m)
1848  {
1849  q.resize (m, k+1);
1850  r.resize (k+1, n+1);
1851  }
1852  else
1853  r.resize (k, n+1);
1854 
1855  F77_INT ldq = to_f77_int (q.rows ());
1856  F77_INT ldr = to_f77_int (r.rows ());
1857 
1858  FloatComplexColumnVector utmp = u;
1859  OCTAVE_LOCAL_BUFFER (float, rw, k);
1860  F77_XFCN (cqrinc, CQRINC, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()), ldq,
1861  F77_CMPLX_ARG (r.fortran_vec ()), ldr, j + 1,
1862  F77_CONST_CMPLX_ARG (utmp.data ()), rw));
1863  }
1864 
1865  template <>
1866  void
1868  const Array<octave_idx_type>& j)
1869  {
1870  F77_INT m = to_f77_int (q.rows ());
1871  F77_INT n = to_f77_int (r.cols ());
1872  F77_INT k = to_f77_int (q.cols ());
1873 
1875  Array<octave_idx_type> js = j.sort (jsi, 0, ASCENDING);
1876  F77_INT nj = to_f77_int (js.numel ());
1877  bool dups = false;
1878  for (F77_INT i = 0; i < nj - 1; i++)
1879  dups = dups && js(i) == js(i+1);
1880 
1881  if (dups)
1882  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1883 
1884  F77_INT u_nel = to_f77_int (u.numel ());
1885  F77_INT u_cols = to_f77_int (u.cols ());
1886 
1887  if (u_nel != m || u_cols != nj)
1888  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
1889 
1890  F77_INT js_beg = to_f77_int (js(0));
1891  F77_INT js_end = to_f77_int (js(nj-1));
1892 
1893  if (nj > 0 && (js_beg < 0 || js_end > n))
1894  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1895 
1896  if (nj > 0)
1897  {
1898  F77_INT kmax = std::min (k + nj, m);
1899  if (k < m)
1900  {
1901  q.resize (m, kmax);
1902  r.resize (kmax, n + nj);
1903  }
1904  else
1905  r.resize (k, n + nj);
1906 
1907  F77_INT ldq = to_f77_int (q.rows ());
1908  F77_INT ldr = to_f77_int (r.rows ());
1909 
1910  OCTAVE_LOCAL_BUFFER (float, rw, kmax);
1911  for (volatile F77_INT i = 0; i < nj; i++)
1912  {
1913  F77_INT ii = i;
1914  F77_INT js_elt = to_f77_int (js(ii));
1915  F77_XFCN (cqrinc, CQRINC, (m, n + ii, std::min (kmax, k + ii),
1916  F77_CMPLX_ARG (q.fortran_vec ()), ldq,
1917  F77_CMPLX_ARG (r.fortran_vec ()), ldr,
1918  js_elt + 1,
1919  F77_CONST_CMPLX_ARG (u.column (jsi(i)).data ()),
1920  rw));
1921  }
1922  }
1923  }
1924 
1925  template <>
1926  void
1928  {
1929  F77_INT j = to_f77_int (j_arg);
1930 
1931  F77_INT m = to_f77_int (q.rows ());
1932  F77_INT k = to_f77_int (r.rows ());
1933  F77_INT n = to_f77_int (r.cols ());
1934 
1935  if (j < 0 || j > n-1)
1936  (*current_liboctave_error_handler) ("qrdelete: index out of range");
1937 
1938  F77_INT ldq = to_f77_int (q.rows ());
1939  F77_INT ldr = to_f77_int (r.rows ());
1940 
1941  OCTAVE_LOCAL_BUFFER (float, rw, k);
1942  F77_XFCN (cqrdec, CQRDEC, (m, n, k, F77_CMPLX_ARG (q.fortran_vec ()), ldq,
1943  F77_CMPLX_ARG (r.fortran_vec ()), ldr, j + 1,
1944  rw));
1945 
1946  if (k < m)
1947  {
1948  q.resize (m, k-1);
1949  r.resize (k-1, n-1);
1950  }
1951  else
1952  r.resize (k, n-1);
1953  }
1954 
1955  template <>
1956  void
1958  {
1959  F77_INT m = to_f77_int (q.rows ());
1960  F77_INT n = to_f77_int (r.cols ());
1961  F77_INT k = to_f77_int (q.cols ());
1962 
1964  Array<octave_idx_type> js = j.sort (jsi, 0, DESCENDING);
1965  F77_INT nj = to_f77_int (js.numel ());
1966  bool dups = false;
1967  for (F77_INT i = 0; i < nj - 1; i++)
1968  dups = dups && js(i) == js(i+1);
1969 
1970  if (dups)
1971  (*current_liboctave_error_handler) ("qrinsert: duplicate index detected");
1972 
1973  F77_INT js_beg = to_f77_int (js(0));
1974  F77_INT js_end = to_f77_int (js(nj-1));
1975 
1976  if (nj > 0 && (js_beg > n-1 || js_end < 0))
1977  (*current_liboctave_error_handler) ("qrinsert: index out of range");
1978 
1979  if (nj > 0)
1980  {
1981  F77_INT ldq = to_f77_int (q.rows ());
1982  F77_INT ldr = to_f77_int (r.rows ());
1983 
1984  OCTAVE_LOCAL_BUFFER (float, rw, k);
1985  for (volatile F77_INT i = 0; i < nj; i++)
1986  {
1987  F77_INT ii = i;
1988  F77_INT js_elt = to_f77_int (js(ii));
1989  F77_XFCN (cqrdec, CQRDEC, (m, n - ii, (k == m ? k : k - ii),
1990  F77_CMPLX_ARG (q.fortran_vec ()), ldq,
1991  F77_CMPLX_ARG (r.fortran_vec ()), ldr,
1992  js_elt + 1, rw));
1993  }
1994 
1995  if (k < m)
1996  {
1997  q.resize (m, k - nj);
1998  r.resize (k - nj, n - nj);
1999  }
2000  else
2001  r.resize (k, n - nj);
2002  }
2003  }
2004 
2005  template <>
2006  void
2008  octave_idx_type j_arg)
2009  {
2010  F77_INT j = to_f77_int (j_arg);
2011 
2012  F77_INT m = to_f77_int (r.rows ());
2013  F77_INT n = to_f77_int (r.cols ());
2014  F77_INT k = std::min (m, n);
2015 
2016  F77_INT u_nel = to_f77_int (u.numel ());
2017 
2018  if (! q.issquare () || u_nel != n)
2019  (*current_liboctave_error_handler) ("qrinsert: dimensions mismatch");
2020 
2021  if (j < 0 || j > m)
2022  (*current_liboctave_error_handler) ("qrinsert: index out of range");
2023 
2024  q.resize (m + 1, m + 1);
2025  r.resize (m + 1, n);
2026 
2027  F77_INT ldq = to_f77_int (q.rows ());
2028  F77_INT ldr = to_f77_int (r.rows ());
2029 
2030  FloatComplexRowVector utmp = u;
2031  OCTAVE_LOCAL_BUFFER (float, rw, k);
2032  F77_XFCN (cqrinr, CQRINR, (m, n, F77_CMPLX_ARG (q.fortran_vec ()), ldq,
2033  F77_CMPLX_ARG (r.fortran_vec ()), ldr,
2034  j + 1, F77_CMPLX_ARG (utmp.fortran_vec ()),
2035  rw));
2036 
2037  }
2038 
2039  template <>
2040  void
2042  {
2043  F77_INT j = to_f77_int (j_arg);
2044 
2045  F77_INT m = to_f77_int (r.rows ());
2046  F77_INT n = to_f77_int (r.cols ());
2047 
2048  if (! q.issquare ())
2049  (*current_liboctave_error_handler) ("qrdelete: dimensions mismatch");
2050 
2051  if (j < 0 || j > m-1)
2052  (*current_liboctave_error_handler) ("qrdelete: index out of range");
2053 
2054  F77_INT ldq = to_f77_int (q.rows ());
2055  F77_INT ldr = to_f77_int (r.rows ());
2056 
2058  OCTAVE_LOCAL_BUFFER (float, rw, m);
2059  F77_XFCN (cqrder, CQRDER, (m, n, F77_CMPLX_ARG (q.fortran_vec ()), ldq,
2060  F77_CMPLX_ARG (r.fortran_vec ()), ldr, j + 1,
2061  F77_CMPLX_ARG (w), rw));
2062 
2063  q.resize (m - 1, m - 1);
2064  r.resize (m - 1, n);
2065  }
2066 
2067  template <>
2068  void
2070  octave_idx_type j_arg)
2071  {
2072  F77_INT i = to_f77_int (i_arg);
2073  F77_INT j = to_f77_int (j_arg);
2074 
2075  F77_INT m = to_f77_int (q.rows ());
2076  F77_INT k = to_f77_int (r.rows ());
2077  F77_INT n = to_f77_int (r.cols ());
2078 
2079  if (i < 0 || i > n-1 || j < 0 || j > n-1)
2080  (*current_liboctave_error_handler) ("qrshift: index out of range");
2081 
2082  F77_INT ldq = to_f77_int (q.rows ());
2083  F77_INT ldr = to_f77_int (r.rows ());
2084 
2086  OCTAVE_LOCAL_BUFFER (float, rw, k);
2087  F77_XFCN (cqrshc, CQRSHC, (m, n, k,
2088  F77_CMPLX_ARG (q.fortran_vec ()), ldq,
2089  F77_CMPLX_ARG (r.fortran_vec ()), ldr,
2090  i + 1, j + 1, F77_CMPLX_ARG (w), rw));
2091  }
2092 
2093 #endif
2094 
2095  // Instantiations we need.
2096 
2097  template class qr<Matrix>;
2098 
2099  template class qr<FloatMatrix>;
2100 
2101  template class qr<ComplexMatrix>;
2102 
2103  template class qr<FloatComplexMatrix>;
2104  }
2105 }
octave_idx_type rows(void) const
Definition: Array.h:404
void resize(octave_idx_type nr, octave_idx_type nc, const FloatComplex &rfv=FloatComplex(0))
Definition: fCMatrix.h:192
qr(void)
Definition: qr.h:53
void resize(octave_idx_type nr, octave_idx_type nc, float rfv=0)
Definition: fMatrix.h:152
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
Definition: dMatrix.h:148
const T * data(void) const
Definition: Array.h:582
static const idx_vector colon
Definition: idx-vector.h:498
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:12177
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:315
for large enough k
Definition: lu.cc:617
const T * fortran_vec(void) const
Definition: Array.h:584
bool regular(void) const
Definition: qr.cc:88
type get_type(void) const
Definition: qr.cc:72
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:488
u
Definition: lu.cc:138
#define F77_XFCN(f, F, args)
Definition: f77-fcn.h:41
void warn_qrupdate_once(void)
octave_idx_type cols(void) const
Definition: Array.h:412
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:400
Definition: mx-defs.h:79
void insert_col(const CV_T &u, octave_idx_type j)
Array< T > sort(int dim=0, sortmode mode=ASCENDING) const
Definition: Array.cc:1756
void update(const CV_T &u, const CV_T &v)
void init(const T &a, type qr_type)
octave_value & assign(assign_op op, const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
static octave::math::qr< T >::type qr_type(int nargout, bool economy)
Definition: qr.cc:56
std::complex< double > w(std::complex< double > z, double relerr=0)
octave_idx_type rows(void) const
Definition: ov.h:472
octave_value retval
Definition: data.cc:6246
ColumnVector column(octave_idx_type i) const
Definition: dMatrix.cc:421
idx type
Definition: ov.cc:3114
Definition: dMatrix.h:36
void form(octave_idx_type n, T &afact, ELT_T *tau, type qr_type)
void delete_row(octave_idx_type j)
T & xelem(octave_idx_type n)
Definition: Array.h:458
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
#define F77_CMPLX_ARG(x)
Definition: f77-fcn.h:309
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition: f77-fcn.h:318
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:187
octave_f77_int_type F77_INT
Definition: f77-fcn.h:305
p
Definition: lu.cc:138
FloatComplexColumnVector column(octave_idx_type i) const
Definition: fCMatrix.cc:710
T::element_type ELT_T
Definition: qr.h:42
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:41
#define F77_CONST_CMPLX_ARG(x)
Definition: f77-fcn.h:312
ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:707
for i
Definition: data.cc:5264
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
std::complex< double > Complex
Definition: oct-cmplx.h:31
FloatColumnVector column(octave_idx_type i) const
Definition: fMatrix.cc:427
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:366
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204