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