GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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()