GNU Octave 7.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-2022 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
53namespace octave
54{
55 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 (void) 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 (void) 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 (void)
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}
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::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:504
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:411
Array::cols
octave_idx_type cols(void) const
Definition: Array.h:457
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition: Array.h:534
Array::rows
octave_idx_type rows(void) const
Definition: Array.h:449
Array::column
OCTARRAY_API Array< T, Alloc > column(octave_idx_type k) const
Extract column: A(:,k+1).
Definition: Array.cc:261
Array::data
const T * data(void) const
Size of the specified dimension.
Definition: Array.h:616
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
Array::sort
OCTARRAY_API Array< T, Alloc > sort(int dim=0, sortmode mode=ASCENDING) const
Size of the specified dimension.
Definition: Array.cc:1761
ComplexMatrix::resize
void resize(octave_idx_type nr, octave_idx_type nc, const Complex &rfv=Complex(0))
Definition: CMatrix.h:193
ComplexMatrix::column
OCTAVE_API ComplexColumnVector column(octave_idx_type i) const
Definition: CMatrix.cc:708
FloatComplexMatrix::column
OCTAVE_API 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
OCTAVE_API 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
OCTAVE_API 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
octave::idx_vector
Definition: idx-vector.h:60
octave::idx_vector::index
octave_idx_type index(const T *src, octave_idx_type n, T *dest) const
Definition: idx-vector.h:602
octave::idx_vector::colon
static const idx_vector colon
Definition: idx-vector.h:483
octave::math::qr::update
OCTAVE_API void update(const CV_T &u, const CV_T &v)
octave::math::qr::delete_col
OCTAVE_API void delete_col(octave_idx_type j)
octave::math::qr::regular
OCTAVE_API bool regular(void) const
Definition: qr.cc:89
octave::math::qr::ELT_T
T::element_type ELT_T
Definition: qr.h:43
octave::math::qr::qr
qr(void)
Definition: qr.h:54
octave::math::qr::get_type
OCTAVE_API type get_type(void) const
Definition: qr.cc:73
octave::math::qr::shift_cols
OCTAVE_API void shift_cols(octave_idx_type i, octave_idx_type j)
octave::math::qr::form
OCTAVE_API void form(octave_idx_type n, T &afact, ELT_T *tau, type qr_type)
octave::math::qr::insert_row
OCTAVE_API void insert_row(const RV_T &u, octave_idx_type j)
octave::math::qr::delete_row
OCTAVE_API void delete_row(octave_idx_type j)
octave::math::qr::insert_col
OCTAVE_API void insert_col(const CV_T &u, octave_idx_type j)
octave::math::qr::init
OCTAVE_API void init(const T &a, type qr_type)
octave_idx_type
qr
Definition: mx-defs.h:53
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
qr_type
static math::qr< T >::type qr_type(int nargout, bool economy)
Definition: qr.cc:74
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.in.cc:55
Matrix
class OCTAVE_API Matrix
Definition: mx-fwd.h:31
ComplexMatrix
class OCTAVE_API ComplexMatrix
Definition: mx-fwd.h:32
FloatComplexMatrix
class OCTAVE_API FloatComplexMatrix
Definition: mx-fwd.h:34
FloatMatrix
class OCTAVE_API FloatMatrix
Definition: mx-fwd.h:33
Faddeeva::w
std::complex< double > w(std::complex< double > z, double relerr=0)
octave::math::warn_qrupdate_once
OCTAVE_API void warn_qrupdate_once(void)
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