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