GNU Octave 10.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-2025 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
53OCTAVE_BEGIN_NAMESPACE(octave)
54
55OCTAVE_BEGIN_NAMESPACE(math)
56
57template <typename T>
58qr<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
71template <typename T>
72typename qr<T>::type
73qr<T>::get_type () const
74{
75 type retval;
76
77 if (! m_q.isempty () && m_q.issquare ())
78 retval = qr<T>::std;
79 else if (m_q.rows () > m_q.cols () && m_r.issquare ())
80 retval = qr<T>::economy;
81 else
82 retval = qr<T>::raw;
83
84 return retval;
85}
86
87template <typename T>
88bool
89qr<T>::regular () const
90{
91 bool retval = true;
92
93 octave_idx_type k = std::min (m_r.rows (), m_r.cols ());
94
95 for (octave_idx_type i = 0; i < k; i++)
96 {
97 if (m_r(i, i) == ELT_T ())
98 {
99 retval = false;
100 break;
101 }
102 }
103
104 return retval;
105}
106
107#if ! defined (HAVE_QRUPDATE)
108
109// Replacement update methods.
110
111void
112warn_qrupdate_once ()
113{
114 static bool warned = false;
115
116 if (! warned)
117 {
118 (*current_liboctave_warning_with_id_handler)
119 ("Octave:missing-dependency",
120 "In this version of Octave, QR & Cholesky updating routines "
121 "simply update the matrix and recalculate factorizations. "
122 "To use fast algorithms, link Octave with the qrupdate library. "
123 "See <http://sourceforge.net/projects/qrupdate>.");
124
125 warned = true;
126 }
127}
128
129template <typename T>
130void
131qr<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
144template <typename T>
145void
146qr<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
159template <typename T, typename CV_T>
160static
161T
162insert_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
173template <typename T, typename RV_T>
174static
175T
176insert_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
187template <typename T>
188static
189T
190delete_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
197template <typename T>
198static
199T
200delete_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
207template <typename T>
208static
209T
210shift_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
229template <typename T>
230void
231qr<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
247template <typename T>
248void
249qr<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
282template <typename T>
283void
284qr<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
296template <typename T>
297void
298qr<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
327template <typename T>
328void
329qr<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
345template <typename T>
346void
347qr<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
362template <typename T>
363void
364qr<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
380template <>
381OCTAVE_API void
382qr<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.rwdata (), 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.rwdata (), m,
446 tau, work, lwork, info));
447 }
448 }
449}
450
451template <>
452OCTAVE_API void
453qr<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.rwdata (), 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.rwdata (), m, tau,
479 work, lwork, info));
480 }
481
482 form (n, afact, tau, qr_type);
483}
484
485#if defined (HAVE_QRUPDATE)
486
487template <>
488OCTAVE_API void
489qr<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.rwdata (), m,
505 m_r.rwdata (), k, utmp.rwdata (),
506 vtmp.rwdata (), w));
507}
508
509template <>
510OCTAVE_API void
511qr<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 (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.rwdata (), m,
532 m_r.rwdata (), k, utmp.rwdata (),
533 vtmp.rwdata (), w));
534 }
535}
536
537template <>
538OCTAVE_API void
539qr<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.rwdata (), ldq,
569 m_r.rwdata (), ldr, j + 1,
570 utmp.data (), w));
571}
572
573template <>
574OCTAVE_API void
575qr<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 (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.rwdata (), ldq,
625 m_r.rwdata (), ldr, js_elt + 1,
626 utmp.data (), w));
627 }
628 }
629}
630
631template <>
632OCTAVE_API void
633qr<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.rwdata (), ldq,
649 m_r.rwdata (), 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
660template <>
661OCTAVE_API void
662qr<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 (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.rwdata (), ldq,
696 m_r.rwdata (), 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
710template <>
711OCTAVE_API void
712qr<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.rwdata (), ldq,
737 m_r.rwdata (), ldr,
738 j + 1, utmp.rwdata (), w));
739
740}
741
742template <>
743OCTAVE_API void
744qr<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.rwdata (), ldq,
762 m_r.rwdata (), ldr, j + 1, w));
763
764 m_q.resize (m - 1, m - 1);
765 m_r.resize (m - 1, n);
766}
767
768template <>
769OCTAVE_API void
770qr<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.rwdata (), ldq,
788 m_r.rwdata (), ldr,
789 i + 1, j + 1, w));
790}
791
792#endif
793
794template <>
795OCTAVE_API void
796qr<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.rwdata (), 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.rwdata (), m,
860 tau, work, lwork, info));
861 }
862 }
863}
864
865template <>
866OCTAVE_API void
867qr<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.rwdata (), 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.rwdata (), m, tau,
893 work, lwork, info));
894 }
895
896 form (n, afact, tau, qr_type);
897}
898
899#if defined (HAVE_QRUPDATE)
900
901template <>
902OCTAVE_API void
903qr<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.rwdata (), m,
919 m_r.rwdata (), k, utmp.rwdata (),
920 vtmp.rwdata (), w));
921}
922
923template <>
924OCTAVE_API void
925qr<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 (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.rwdata (), m,
946 m_r.rwdata (), k, utmp.rwdata (),
947 vtmp.rwdata (), w));
948 }
949}
950
951template <>
952OCTAVE_API void
953qr<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.rwdata (), ldq,
984 m_r.rwdata (), ldr, j + 1,
985 utmp.data (), w));
986}
987
988template <>
989OCTAVE_API void
990qr<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 (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.rwdata (), ldq,
1041 m_r.rwdata (), ldr, js_elt + 1,
1042 utmp.data (), w));
1043 }
1044 }
1045}
1046
1047template <>
1048OCTAVE_API void
1049qr<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.rwdata (), ldq,
1065 m_r.rwdata (), 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
1076template <>
1077OCTAVE_API void
1078qr<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 (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.rwdata (), ldq,
1112 m_r.rwdata (), 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
1126template <>
1127OCTAVE_API void
1128qr<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.rwdata (), ldq,
1154 m_r.rwdata (), ldr,
1155 j + 1, utmp.rwdata (), w));
1156
1157}
1158
1159template <>
1160OCTAVE_API void
1161qr<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.rwdata (), ldq,
1179 m_r.rwdata (), ldr, j + 1,
1180 w));
1181
1182 m_q.resize (m - 1, m - 1);
1183 m_r.resize (m - 1, n);
1184}
1185
1186template <>
1187OCTAVE_API void
1188qr<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.rwdata (), ldq,
1206 m_r.rwdata (), ldr,
1207 i + 1, j + 1, w));
1208}
1209
1210#endif
1211
1212template <>
1213OCTAVE_API void
1214qr<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.rwdata ()),
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.rwdata ()),
1282 m, F77_DBLE_CMPLX_ARG (tau),
1283 F77_DBLE_CMPLX_ARG (work), lwork,
1284 info));
1285 }
1286 }
1287}
1288
1289template <>
1290OCTAVE_API void
1291qr<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.rwdata ()),
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.rwdata ()),
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
1329template <>
1330OCTAVE_API void
1331qr<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.rwdata ()),
1349 m, F77_DBLE_CMPLX_ARG (m_r.rwdata ()), k,
1350 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
1351 F77_DBLE_CMPLX_ARG (vtmp.rwdata ()),
1352 F77_DBLE_CMPLX_ARG (w), rw));
1353}
1354
1355template <>
1356OCTAVE_API void
1357qr<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 (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.rwdata ()),
1380 m, F77_DBLE_CMPLX_ARG (m_r.rwdata ()), k,
1381 F77_DBLE_CMPLX_ARG (utmp.rwdata ()),
1382 F77_DBLE_CMPLX_ARG (vtmp.rwdata ()),
1383 F77_DBLE_CMPLX_ARG (w), rw));
1384 }
1385}
1386
1387template <>
1388OCTAVE_API void
1389qr<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.rwdata ()),
1420 ldq, F77_DBLE_CMPLX_ARG (m_r.rwdata ()),
1421 ldr, j + 1,
1422 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()), rw));
1423}
1424
1425template <>
1426OCTAVE_API void
1427qr<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 (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.rwdata ()),
1478 ldq,
1479 F77_DBLE_CMPLX_ARG (m_r.rwdata ()),
1480 ldr, js_elt + 1,
1481 F77_CONST_DBLE_CMPLX_ARG (utmp.data ()),
1482 rw));
1483 }
1484 }
1485}
1486
1487template <>
1488OCTAVE_API void
1489qr<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.rwdata ()),
1505 ldq, F77_DBLE_CMPLX_ARG (m_r.rwdata ()),
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
1517template <>
1518OCTAVE_API void
1519qr<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 (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.rwdata ()),
1553 ldq,
1554 F77_DBLE_CMPLX_ARG (m_r.rwdata ()),
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
1568template <>
1569OCTAVE_API void
1570qr<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.rwdata ()),
1596 ldq, F77_DBLE_CMPLX_ARG (m_r.rwdata ()),
1597 ldr, j + 1,
1598 F77_DBLE_CMPLX_ARG (utmp.rwdata ()), rw));
1599
1600}
1601
1602template <>
1603OCTAVE_API void
1604qr<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.rwdata ()),
1623 ldq, F77_DBLE_CMPLX_ARG (m_r.rwdata ()),
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
1630template <>
1631OCTAVE_API void
1632qr<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.rwdata ()), ldq,
1652 F77_DBLE_CMPLX_ARG (m_r.rwdata ()), ldr,
1653 i + 1, j + 1, F77_DBLE_CMPLX_ARG (w), rw));
1654}
1655
1656#endif
1657
1658template <>
1659OCTAVE_API void
1660qr<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.rwdata ()), 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.rwdata ()), m,
1727 F77_CMPLX_ARG (tau),
1728 F77_CMPLX_ARG (work), lwork, info));
1729 }
1730 }
1731}
1732
1733template <>
1734OCTAVE_API void
1735qr<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.rwdata ()),
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.rwdata ()),
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
1771template <>
1772OCTAVE_API void
1773qr<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.rwdata ()),
1791 m, F77_CMPLX_ARG (m_r.rwdata ()), k,
1792 F77_CMPLX_ARG (utmp.rwdata ()),
1793 F77_CMPLX_ARG (vtmp.rwdata ()),
1794 F77_CMPLX_ARG (w), rw));
1795}
1796
1797template <>
1798OCTAVE_API void
1799qr<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 (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.rwdata ()),
1822 m, F77_CMPLX_ARG (m_r.rwdata ()), k,
1823 F77_CMPLX_ARG (utmp.rwdata ()),
1824 F77_CMPLX_ARG (vtmp.rwdata ()),
1825 F77_CMPLX_ARG (w), rw));
1826 }
1827}
1828
1829template <>
1830OCTAVE_API void
1831qr<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.rwdata ()), ldq,
1862 F77_CMPLX_ARG (m_r.rwdata ()), ldr, j + 1,
1863 F77_CONST_CMPLX_ARG (utmp.data ()), rw));
1864}
1865
1866template <>
1867OCTAVE_API void
1868qr<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 (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.rwdata ()), ldq,
1918 F77_CMPLX_ARG (m_r.rwdata ()), ldr,
1919 js_elt + 1,
1920 F77_CONST_CMPLX_ARG (u.column (jsi(i)).data ()),
1921 rw));
1922 }
1923 }
1924}
1925
1926template <>
1927OCTAVE_API void
1928qr<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.rwdata ()), ldq,
1944 F77_CMPLX_ARG (m_r.rwdata ()), 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
1956template <>
1957OCTAVE_API void
1958qr<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 (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.rwdata ()), ldq,
1992 F77_CMPLX_ARG (m_r.rwdata ()), 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
2006template <>
2007OCTAVE_API void
2008qr<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.rwdata ()), ldq,
2034 F77_CMPLX_ARG (m_r.rwdata ()), ldr,
2035 j + 1, F77_CMPLX_ARG (utmp.rwdata ()),
2036 rw));
2037
2038}
2039
2040template <>
2041OCTAVE_API void
2042qr<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.rwdata ()), ldq,
2061 F77_CMPLX_ARG (m_r.rwdata ()), 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
2068template <>
2069OCTAVE_API void
2070qr<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.rwdata ()), ldq,
2090 F77_CMPLX_ARG (m_r.rwdata ()), ldr,
2091 i + 1, j + 1, F77_CMPLX_ARG (w), rw));
2092}
2093
2094#endif
2095
2096// Instantiations we need.
2097
2098template class qr<Matrix>;
2099
2100template class qr<FloatMatrix>;
2101
2102template class qr<ComplexMatrix>;
2103
2104template class qr<FloatComplexMatrix>;
2105
2106OCTAVE_END_NAMESPACE(math)
2107OCTAVE_END_NAMESPACE(octave)
Array
N Dimensional Array with copy-on-write semantics.
Definition Array.h:130
Array::xelem
T & xelem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:525
Array::elem
T & elem(octave_idx_type n)
Size of the specified dimension.
Definition Array.h:563
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.h:463
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::sort
Array< T, Alloc > sort(int dim=0, sortmode mode=ASCENDING) const
Size of the specified dimension.
Definition Array-base.cc:1804
Array::data
const T * data() const
Size of the specified dimension.
Definition Array.h:665
Array::rwdata
T * rwdata()
Size of the specified dimension.
Definition Array-base.cc:1787
Array::numel
octave_idx_type numel() const
Number of elements in the array.
Definition Array.h:418
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:90
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:584
octave_idx_type
qr
Definition qr.h:39
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:42
qr::delete_row
void delete_row(octave_idx_type j)
qr::qr
qr()
Definition qr.h:53
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:47
qr::regular
bool regular() const
Definition qr.cc:89
qr::form
void form(octave_idx_type n, T &afact, ELT_T *tau, type qr_type)
qr::insert_col
void insert_col(const CV_T &u, octave_idx_type j)
qr::get_type
type get_type() const
Definition qr.cc:73
qr::shift_cols
void shift_cols(octave_idx_type i, octave_idx_type j)
qr::m_r
T m_r
Definition qr.h:109
qr::init
void init(const T &a, type qr_type)
qr::m_q
T m_q
Definition qr.h:108
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
F77_CONST_CMPLX_ARG
#define F77_CONST_CMPLX_ARG(x)
Definition f77-fcn.h:313
F77_DBLE_CMPLX_ARG
#define F77_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:316
F77_CMPLX_ARG
#define F77_CMPLX_ARG(x)
Definition f77-fcn.h:310
F77_XFCN
#define F77_XFCN(f, F, args)
Definition f77-fcn.h:45
F77_INT
octave_f77_int_type F77_INT
Definition f77-fcn.h:306
F77_CONST_DBLE_CMPLX_ARG
#define F77_CONST_DBLE_CMPLX_ARG(x)
Definition f77-fcn.h:319
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition lo-error.c:41
x
F77_RET_T const F77_DBLE * x
Definition lo-slatec-proto.h:38
OCTAVE_API
#define OCTAVE_API
Definition main.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()