GNU Octave 10.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 
Loading...
Searching...
No Matches
sparse-qr.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 2005-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 "CMatrix.h"
31#include "CSparse.h"
32#include "MArray.h"
33#include "dColVector.h"
34#include "dMatrix.h"
35#include "dSparse.h"
36#include "lo-error.h"
37#include "oct-locbuf.h"
38#include "oct-sparse.h"
39#include "quit.h"
40#include "sparse-qr.h"
41
42OCTAVE_BEGIN_NAMESPACE(octave)
43
44OCTAVE_BEGIN_NAMESPACE(math)
45
46#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
47// Decide once at runtime whether Octave must workaround SuiteSparse library.
48static constexpr bool octave_suitesparse_ptr_size_mismatch
49 = (sizeof (octave_idx_type) != sizeof (SuiteSparse_long));
50
51static constexpr bool suitesparse_integer_long_mismatch
52 = (sizeof (suitesparse_integer) != sizeof (SuiteSparse_long));
53#endif
54
55#if defined (HAVE_CXSPARSE)
56template <typename SPARSE_T>
57class cxsparse_types
58{ };
59
60template <>
61class cxsparse_types<SparseMatrix>
62{
63public:
64 typedef CXSPARSE_DNAME (s) symbolic_type;
65 typedef CXSPARSE_DNAME (n) numeric_type;
66};
67
68template <>
69class cxsparse_types<SparseComplexMatrix>
70{
71public:
72 typedef CXSPARSE_ZNAME (s) symbolic_type;
73 typedef CXSPARSE_ZNAME (n) numeric_type;
74};
75#endif
76
77template <typename SPARSE_T>
78class sparse_qr<SPARSE_T>::sparse_qr_rep
79{
80public:
81
82 sparse_qr_rep (const SPARSE_T& a, int order);
83
84 OCTAVE_DISABLE_CONSTRUCT_COPY_MOVE (sparse_qr_rep)
85
86 ~sparse_qr_rep ();
87
88 bool ok () const
89 {
90#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
91 return (m_H && m_Htau && m_HPinv && m_R && m_E);
92#elif defined (HAVE_CXSPARSE)
93 return (N && S);
94#else
95 return false;
96#endif
97 }
98
99 SPARSE_T V () const;
100
101 ColumnVector Pinv () const;
102
103 ColumnVector P () const;
104
105 ColumnVector E () const;
106
107 SPARSE_T R (bool econ) const;
108
109 typename SPARSE_T::dense_matrix_type
110 C (const typename SPARSE_T::dense_matrix_type& b, bool econ = false);
111
112 typename SPARSE_T::dense_matrix_type Q (bool econ = false);
113
114 octave_idx_type nrows;
115 octave_idx_type ncols;
116
117#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
118
119 template <typename RHS_T, typename RET_T>
120 RET_T solve (const RHS_T& b, octave_idx_type& info) const;
121
122#endif
123
124#if defined (HAVE_CXSPARSE)
125
126 typename cxsparse_types<SPARSE_T>::symbolic_type *S;
127 typename cxsparse_types<SPARSE_T>::numeric_type *N;
128
129#endif
130
131 template <typename RHS_T, typename RET_T>
132 RET_T tall_solve (const RHS_T& b, octave_idx_type& info);
133
134 template <typename RHS_T, typename RET_T>
135 RET_T wide_solve (const RHS_T& b, octave_idx_type& info) const;
136
137#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
138
139private:
140
141 cholmod_common m_cc;
142 cholmod_sparse *m_R; // R factor
143 // Column permutation for A. Fill-reducing ordering.
144 SuiteSparse_long *m_E;
145 cholmod_sparse *m_H; // Householder vectors
146 cholmod_dense *m_Htau; // beta scalars
147 SuiteSparse_long *m_HPinv;
148
149#endif
150};
151
152template <typename SPARSE_T>
153ColumnVector
154sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv () const
155{
156#if defined (HAVE_CXSPARSE)
157
158 ColumnVector ret (N->L->m);
159
160 for (octave_idx_type i = 0; i < N->L->m; i++)
161 ret.xelem (i) = S->pinv[i];
162
163 return ret;
164
165#else
166
167 return ColumnVector ();
168
169#endif
170}
171
172template <typename SPARSE_T>
173ColumnVector
174sparse_qr<SPARSE_T>::sparse_qr_rep::P () const
175{
176#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
177
178 ColumnVector ret (nrows);
179
180 // FIXME: Is ret.xelem (m_HPinv[i]) = i + 1 correct?
181 for (octave_idx_type i = 0; i < nrows; i++)
182 ret.xelem (from_suitesparse_long (m_HPinv[i])) = i + 1;
183
184 return ret;
185
186#elif defined (HAVE_CXSPARSE)
187
188 ColumnVector ret (N->L->m);
189
190 for (octave_idx_type i = 0; i < N->L->m; i++)
191 ret.xelem (S->pinv[i]) = i;
192
193 return ret;
194
195#else
196
197 return ColumnVector ();
198
199#endif
200}
201
202template <typename SPARSE_T>
203ColumnVector
204sparse_qr<SPARSE_T>::sparse_qr_rep::E () const
205{
206#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
207
208 ColumnVector ret (ncols);
209
210 if (m_E)
211 for (octave_idx_type i = 0; i < ncols; i++)
212 ret(i) = from_suitesparse_long (m_E[i]) + 1;
213 else
214 for (octave_idx_type i = 0; i < ncols; i++)
215 ret(i) = i + 1;
216
217 return ret;
218
219#elif defined (HAVE_CXSPARSE)
220
221 ColumnVector ret (ncols);
222
223 if (S->q)
224 for (octave_idx_type i = 0; i < ncols; i++)
225 ret(i) = from_suitesparse_long (S->q[i]) + 1;
226 else
227 for (octave_idx_type i = 0; i < ncols; i++)
228 ret(i) = i + 1;
229
230 return ret;
231
232#else
233
234 return ColumnVector ();
235
236#endif
237}
238
239#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
240
241// Convert octave real sparse matrix to cholmod real sparse matrix.
242// Returns a "shallow" copy of a.
243static cholmod_sparse
244ors2crs (const SparseMatrix& a)
245{
246 cholmod_sparse A;
247
248 octave_idx_type ncols = a.cols ();
249 octave_idx_type nnz = a.nnz ();
250
251 A.ncol = ncols;
252 A.nrow = a.rows ();
253 A.itype = CHOLMOD_LONG;
254 A.nzmax = nnz;
255 A.sorted = 0;
256 A.packed = 1;
257 A.stype = 0;
258 A.xtype = CHOLMOD_REAL;
259 A.dtype = CHOLMOD_DOUBLE;
260 A.nz = nullptr;
261 A.z = nullptr;
262 if (! octave_suitesparse_ptr_size_mismatch)
263 {
264 A.p = reinterpret_cast<SuiteSparse_long *> (a.cidx ());
265 A.i = reinterpret_cast<SuiteSparse_long *> (a.ridx ());
266 }
267 else
268 {
269 SuiteSparse_long *A_p;
270 A_p = new SuiteSparse_long[ncols+1];
271 for (octave_idx_type i = 0; i < ncols+1; i++)
272 A_p[i] = a.cidx (i);
273 A.p = A_p;
274 SuiteSparse_long *A_i;
275 A_i = new SuiteSparse_long[nnz];
276 for (octave_idx_type i = 0; i < nnz; i++)
277 A_i[i] = a.ridx (i);
278 A.i = A_i;
279 }
280 A.x = const_cast<double *> (a.data ());
281
282 return A;
283}
284
285// Convert octave complex sparse matrix to cholmod complex sparse matrix.
286// Returns a "shallow" copy of a.
287static cholmod_sparse
288ocs2ccs (const SparseComplexMatrix& a)
289{
290 cholmod_sparse A;
291
292 octave_idx_type ncols = a.cols ();
293 octave_idx_type nnz = a.nnz ();
294
295 A.ncol = ncols;
296 A.nrow = a.rows ();
297 A.itype = CHOLMOD_LONG;
298 A.nzmax = nnz;
299 A.sorted = 0;
300 A.packed = 1;
301 A.stype = 0;
302 A.xtype = CHOLMOD_COMPLEX;
303 A.dtype = CHOLMOD_DOUBLE;
304 A.nz = nullptr;
305 A.z = nullptr;
306 if (! octave_suitesparse_ptr_size_mismatch)
307 {
308 A.p = reinterpret_cast<SuiteSparse_long *> (a.cidx ());
309 A.i = reinterpret_cast<SuiteSparse_long *> (a.ridx ());
310 }
311 else
312 {
313 SuiteSparse_long *A_p;
314 A_p = new SuiteSparse_long[ncols+1];
315 for (octave_idx_type i = 0; i < ncols+1; i++)
316 A_p[i] = a.cidx (i);
317 A.p = A_p;
318 SuiteSparse_long *A_i;
319 A_i = new SuiteSparse_long[nnz];
320 for (octave_idx_type i = 0; i < nnz; i++)
321 A_i[i] = a.ridx (i);
322 A.i = A_i;
323 }
324 A.x = const_cast<Complex *>
325 (reinterpret_cast<const Complex *> (a.data ()));
326
327 return A;
328}
329
330// Convert octave real dense matrix to cholmod complex dense matrix.
331// Returns a "deep" copy of a.
332static cholmod_dense *
333ord2ccd (const MArray<double>& a, cholmod_common *cc1)
334{
335 cholmod_dense *A
336 = cholmod_l_allocate_dense (a.rows (), a.cols (), a.rows(),
337 CHOLMOD_COMPLEX, cc1);
338
339 const double *a_x = a.data ();
340
341 Complex *A_x = reinterpret_cast<Complex *> (A->x);
342 for (octave_idx_type j = 0; j < a.cols() * a.rows() ; j++)
343 A_x[j] = Complex (a_x[j], 0.0);
344
345 return A;
346}
347
348// Convert octave real dense matrix to cholmod real dense matrix.
349// Returns a "shallow" copy of a.
350static cholmod_dense
351ord2crd (const MArray<double>& a)
352{
353 cholmod_dense A;
354
355 A.ncol = a.cols ();
356 A.nrow = a.rows ();
357 A.nzmax = a.cols () * a.rows ();
358 A.xtype = CHOLMOD_REAL;
359 A.dtype = CHOLMOD_DOUBLE;
360 A.z = nullptr;
361 A.d = a.rows ();
362 A.x = const_cast<double *> (a.data ());
363
364 return A;
365}
366
367// Convert octave complex dense matrix to cholmod complex dense matrix.
368// Returns a "shallow" copy of a.
369static cholmod_dense
370ocd2ccd (const ComplexMatrix& a)
371{
372 cholmod_dense A;
373
374 A.ncol = a.cols ();
375 A.nrow = a.rows ();
376 A.nzmax = a.cols () * a.rows ();
377 A.xtype = CHOLMOD_COMPLEX;
378 A.dtype = CHOLMOD_DOUBLE;
379 A.z = nullptr;
380 A.d = a.rows ();
381 A.x = const_cast<Complex *> (reinterpret_cast<const Complex *> (a.data ()));
382
383 return A;
384}
385
386// Convert cholmod real sparse matrix to octave real sparse matrix.
387// Returns a "shallow" copy of y.
388static SparseMatrix
389crs2ors (cholmod_sparse *y, const cholmod_common *cc1)
390{
391 octave_idx_type nrow = from_size_t (y->nrow);
392 octave_idx_type ncol = from_size_t (y->ncol);
393 octave_idx_type nz
394 = from_suitesparse_long
395 (cholmod_l_nnz (y, const_cast<cholmod_common *> (cc1)));
396 SparseMatrix ret (nrow, ncol, nz);
397
398 SuiteSparse_long *y_p = reinterpret_cast<SuiteSparse_long *> (y->p);
399 for (octave_idx_type j = 0; j < ncol + 1; j++)
400 ret.xcidx (j) = from_suitesparse_long (y_p[j]);
401
402 SuiteSparse_long *y_i = reinterpret_cast<SuiteSparse_long *> (y->i);
403 double *y_x = reinterpret_cast<double *> (y->x);
404 for (octave_idx_type j = 0; j < nz; j++)
405 {
406 ret.xridx (j) = from_suitesparse_long (y_i[j]);
407 ret.xdata (j) = y_x[j];
408 }
409
410 return ret;
411}
412
413// Convert cholmod complex sparse matrix to octave complex sparse matrix.
414// Returns a "deep" copy of a.
415static SparseComplexMatrix
416ccs2ocs (cholmod_sparse *a, cholmod_common *cc1)
417{
418 octave_idx_type nrow = from_size_t (a->nrow);
419 octave_idx_type ncol = from_size_t (a->ncol);
420 octave_idx_type nz = from_suitesparse_long (cholmod_l_nnz (a, cc1));
421 SparseComplexMatrix ret (nrow, ncol, nz);
422
423 SuiteSparse_long *a_p = reinterpret_cast<SuiteSparse_long *> (a->p);
424 for (octave_idx_type j = 0; j < ncol + 1; j++)
425 ret.xcidx(j) = from_suitesparse_long (a_p[j]);
426
427 SuiteSparse_long *a_i = reinterpret_cast<SuiteSparse_long *> (a->i);
428 Complex *a_x = reinterpret_cast<Complex *> (a->x);
429 for (octave_idx_type j = 0; j < nz; j++)
430 {
431 ret.xridx(j) = from_suitesparse_long (a_i[j]);
432 ret.xdata(j) = a_x[j];
433 }
434
435 return ret;
436}
437
438// Convert octave real sparse matrix to cholmod complex sparse matrix.
439// Returns a "deep" copy of a.
440static cholmod_sparse *
441ors2ccs (const SparseMatrix& a, cholmod_common *cc)
442{
443 cholmod_sparse *A
444 = cholmod_l_allocate_sparse (a.rows (), a.cols (), a.nnz (), 0, 1, 0,
445 CHOLMOD_COMPLEX, cc);
446
447 octave_idx_type ncol = a.cols ();
448 SuiteSparse_long *A_p = reinterpret_cast<SuiteSparse_long *> (A->p);
449 for (octave_idx_type j = 0; j < ncol + 1; j++)
450 A_p[j] = a.cidx(j);
451
452 const double *a_x = a.data ();
453 Complex *A_x = reinterpret_cast<Complex *> (A->x);
454 SuiteSparse_long *A_i = reinterpret_cast<SuiteSparse_long *> (A->i);
455 for (octave_idx_type j = 0; j < a.nnz (); j++)
456 {
457 A_x[j] = Complex (a_x[j], 0.0);
458 A_i[j] = a.ridx(j);
459 }
460 return A;
461}
462
463#ifdef HAVE_CXSPARSE
464static suitesparse_integer
465suitesparse_long_to_suitesparse_integer (SuiteSparse_long x)
466{
467 if (x < std::numeric_limits<suitesparse_integer>::min ()
468 || x > std::numeric_limits<suitesparse_integer>::max ())
469 (*current_liboctave_error_handler)
470 ("integer dimension or index out of range for SuiteSparse's indexing type");
471
472 return static_cast<suitesparse_integer> (x);
473}
474#endif
475
476static void
477spqr_error_handler (const cholmod_common *cc)
478{
479 if (cc->status >= 0)
480 return;
481
482 switch (cc->status)
483 {
484 case CHOLMOD_OUT_OF_MEMORY:
485 (*current_liboctave_error_handler)
486 ("sparse_qr: sparse matrix QR factorization failed"
487 " - out of memory");
488 case CHOLMOD_TOO_LARGE:
489 (*current_liboctave_error_handler)
490 ("sparse_qr: sparse matrix QR factorization failed"
491 " - integer overflow occurred");
492 default:
493 (*current_liboctave_error_handler)
494 ("sparse_qr: sparse matrix QR factorization failed"
495 " - error %d", cc->status);
496 }
497
498 // FIXME: Free memory?
499 // FIXME: Can cc-status > 0 (CHOLMOD_NOT_POSDEF, CHOLMOD_DSMALL) occur?
500}
501#endif
502
503// Specializations.
504
505// Real-valued matrices.
506
507// Arguments for parameter order (taken from SuiteSparseQR documentation).
508// 0: fixed ordering 0 (no permutation of columns)
509// 1: natural ordering 1 (only singleton columns are permuted to the left of
510// the matrix)
511// 2: colamd
512// 3:
513// 4: CHOLMOD best-effort (COLAMD, METIS,...)
514// 5: AMD(a'*a)
515// 6: metis(a'*a)
516// 7: SuiteSparseQR default ordering
517// 8: try COLAMD, AMD, and METIS; pick best
518// 9: try COLAMD and AMD; pick best
519//FIXME: What is order = 3?
520template <>
521sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep
522(const SparseMatrix& a, int order)
523 : nrows (a.rows ()), ncols (a.columns ())
524#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
525 , m_cc (), m_R (nullptr), m_E (nullptr), m_H (nullptr), m_Htau (nullptr),
526 m_HPinv (nullptr)
527{
528 octave_idx_type nr = a.rows ();
529 octave_idx_type nc = a.cols ();
530
531 if (nr < 0 || nc < 0)
532 (*current_liboctave_error_handler)
533 ("matrix dimension with negative size");
534
535 if (order < 0 || order > 9)
536 (*current_liboctave_error_handler)
537 ("ordering %d is not supported by SPQR", order);
538
539 cholmod_l_start (&m_cc);
540 cholmod_sparse A = ors2crs (a);
541
542 SuiteSparseQR<double> (order, static_cast<double> (SPQR_DEFAULT_TOL),
543 static_cast<SuiteSparse_long> (A.nrow),
544 &A, &m_R, &m_E, &m_H, &m_HPinv, &m_Htau, &m_cc);
545 spqr_error_handler (&m_cc);
546
547 if (octave_suitesparse_ptr_size_mismatch)
548 {
549 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
550 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
551 }
552}
553
554#elif defined (HAVE_CXSPARSE)
555 , S (nullptr), N (nullptr)
556{
557 CXSPARSE_DNAME () A;
558
559 A.nzmax = a.nnz ();
560 A.m = nrows;
561 A.n = ncols;
562 // Cast away const on A, with full knowledge that CSparse won't touch it
563 // Prevents the methods below making a copy of the data.
564 A.p = const_cast<suitesparse_integer *>
565 (to_suitesparse_intptr (a.cidx ()));
566 A.i = const_cast<suitesparse_integer *>
567 (to_suitesparse_intptr (a.ridx ()));
568 A.x = const_cast<double *> (a.data ());
569 A.nz = -1;
570
571 S = CXSPARSE_DNAME (_sqr) (order, &A, 1);
572 N = CXSPARSE_DNAME (_qr) (&A, S);
573
574 if (! N)
575 (*current_liboctave_error_handler)
576 ("sparse_qr: sparse matrix QR factorization filled");
577
578}
579
580#else
581
582{
583 octave_unused_parameter (order);
584
585 (*current_liboctave_error_handler)
586 ("sparse_qr: support for SPQR or CXSparse was unavailable or disabled when liboctave was built");
587}
588
589#endif
590
591template <>
592sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep ()
593{
594#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
595
596 cholmod_l_free_sparse (&m_R, &m_cc);
597 cholmod_l_free_sparse (&m_H, &m_cc);
598 cholmod_l_free_dense (&m_Htau, &m_cc);
599 free (m_E); // FIXME: use cholmod_l_free
600 free (m_HPinv);
601 cholmod_l_finish (&m_cc);
602
603#elif defined (HAVE_CXSPARSE)
604
605 CXSPARSE_DNAME (_sfree) (S);
606 CXSPARSE_DNAME (_nfree) (N);
607
608#endif
609}
610
611template <>
612SparseMatrix
613sparse_qr<SparseMatrix>::sparse_qr_rep::V () const
614{
615#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
616
617 return crs2ors (m_H, &m_cc);
618
619#elif defined (HAVE_CXSPARSE)
620
621 // Drop zeros from V and sort
622 // FIXME: Is the double transpose to sort necessary?
623
624 CXSPARSE_DNAME (_dropzeros) (N->L);
625 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1);
626 CXSPARSE_DNAME (_spfree) (N->L);
627 N->L = CXSPARSE_DNAME (_transpose) (D, 1);
628 CXSPARSE_DNAME (_spfree) (D);
629
630 octave_idx_type nc = N->L->n;
631 octave_idx_type nz = N->L->p[nc];
632 SparseMatrix ret (N->L->m, nc, nz);
633
634 for (octave_idx_type j = 0; j < nc+1; j++)
635 ret.xcidx (j) = N->L->p[j];
636
637 for (octave_idx_type j = 0; j < nz; j++)
638 {
639 ret.xridx (j) = N->L->i[j];
640 ret.xdata (j) = N->L->x[j];
641 }
642
643 return ret;
644
645#else
646
647 return SparseMatrix ();
648
649#endif
650}
651
652template <>
653SparseMatrix
654sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const
655{
656#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
657
658 octave_idx_type nr = from_size_t (m_R->nrow);
659 octave_idx_type nc = from_size_t (m_R->ncol);
660 octave_idx_type nz
661 = from_suitesparse_long
662 (cholmod_l_nnz (m_R, const_cast<cholmod_common *> (&m_cc)));
663
664 // FIXME: Does this work if econ = true?
665 SparseMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz);
666 SuiteSparse_long *Rp = reinterpret_cast<SuiteSparse_long *> (m_R->p);
667 SuiteSparse_long *Ri = reinterpret_cast<SuiteSparse_long *> (m_R->i);
668
669 for (octave_idx_type j = 0; j < nc + 1; j++)
670 ret.xcidx (j) = from_suitesparse_long (Rp[j]);
671
672 for (octave_idx_type j = 0; j < nz; j++)
673 {
674 ret.xridx (j) = from_suitesparse_long (Ri[j]);
675 ret.xdata (j) = (reinterpret_cast<double *> (m_R->x))[j];
676 }
677
678 return ret;
679
680#elif defined (HAVE_CXSPARSE)
681
682 // Drop zeros from R and sort
683 // FIXME: Is the double transpose to sort necessary?
684
685 CXSPARSE_DNAME (_dropzeros) (N->U);
686 CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1);
687 CXSPARSE_DNAME (_spfree) (N->U);
688 N->U = CXSPARSE_DNAME (_transpose) (D, 1);
689 CXSPARSE_DNAME (_spfree) (D);
690
691 octave_idx_type nc = N->U->n;
692 octave_idx_type nz = N->U->p[nc];
693
694 SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz);
695
696 for (octave_idx_type j = 0; j < nc+1; j++)
697 ret.xcidx (j) = N->U->p[j];
698
699 for (octave_idx_type j = 0; j < nz; j++)
700 {
701 ret.xridx (j) = N->U->i[j];
702 ret.xdata (j) = N->U->x[j];
703 }
704
705 return ret;
706
707#else
708
709 octave_unused_parameter (econ);
710
711 return SparseMatrix ();
712
713#endif
714}
715
716template <>
717Matrix
718sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b, bool econ)
719{
720#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
721 octave_idx_type nr = (econ
722 ? (ncols > nrows ? nrows : ncols)
723 : nrows);
724 octave_idx_type b_nr = b.rows ();
725 octave_idx_type b_nc = b.cols ();
726 Matrix ret (nr, b_nc);
727
728 if (nrows != b_nr)
729 (*current_liboctave_error_handler)
730 ("sparse_qr: matrix dimension mismatch");
731 else if (b_nc < 0 || b_nr < 0)
732 (*current_liboctave_error_handler)
733 ("sparse_qr: matrix dimension with negative size");
734
735 cholmod_dense *QTB; // Q' * B
736 cholmod_dense B = ord2crd (b);
737
738 QTB = SuiteSparseQR_qmult<double> (SPQR_QTX, m_H, m_Htau, m_HPinv, &B,
739 &m_cc);
740 spqr_error_handler (&m_cc);
741
742 // copy QTB into ret
743 double *QTB_x = reinterpret_cast<double *> (QTB->x);
744 double *ret_vec = reinterpret_cast<double *> (ret.rwdata ());
745 for (octave_idx_type j = 0; j < b_nc; j++)
746 for (octave_idx_type i = 0; i < nr; i++)
747 ret_vec[j * nr + i] = QTB_x[j * b_nr + i];
748
749 cholmod_l_free_dense (&QTB, &m_cc);
750
751 return ret;
752
753#elif defined (HAVE_CXSPARSE)
754
755 if (econ)
756 (*current_liboctave_error_handler)
757 ("sparse-qr: economy mode with CXSparse not supported");
758
759 octave_idx_type b_nr = b.rows ();
760 octave_idx_type b_nc = b.cols ();
761
762 octave_idx_type nc = N->L->n;
763 octave_idx_type nr = nrows;
764
765 const double *bvec = b.data ();
766
767 Matrix ret (b_nr, b_nc);
768 double *vec = ret.rwdata ();
769
770 if (nr < 0 || nc < 0 || nr != b_nr)
771 (*current_liboctave_error_handler) ("matrix dimension mismatch");
772
773 if (nr == 0 || nc == 0 || b_nc == 0)
774 ret = Matrix (nc, b_nc, 0.0);
775 else
776 {
777 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
778
779 for (octave_idx_type j = 0, idx = 0;
780 j < b_nc;
781 j++, idx += b_nr)
782 {
783 octave_quit ();
784
785 for (octave_idx_type i = nr; i < S->m2; i++)
786 buf[i] = 0.;
787
788 octave_idx_type nm = (nr < nc ? nr : nc);
789
790 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr);
791
792 for (octave_idx_type i = 0; i < nm; i++)
793 {
794 octave_quit ();
795
796 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
797 }
798
799 for (octave_idx_type i = 0; i < b_nr; i++)
800 vec[i+idx] = buf[i];
801 }
802 }
803
804 return ret;
805
806#else
807
808 octave_unused_parameter (b);
809 octave_unused_parameter (econ);
810
811 return Matrix ();
812
813#endif
814}
815
816template <>
817Matrix
818sparse_qr<SparseMatrix>::sparse_qr_rep::Q (bool econ)
819{
820#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
821
822 octave_idx_type nc = (econ
823 ? (ncols > nrows ? nrows : ncols)
824 : nrows);
825 Matrix ret (nrows, nc);
826 cholmod_dense *q;
827
828 // I is nrows x nrows identity matrix
829 cholmod_dense *I_mat
830 = cholmod_l_allocate_dense (nrows, nrows, nrows, CHOLMOD_REAL, &m_cc);
831
832 for (octave_idx_type i = 0; i < nrows * nrows; i++)
833 (reinterpret_cast<double *> (I_mat->x))[i] = 0.0;
834
835 for (octave_idx_type i = 0; i < nrows; i++)
836 (reinterpret_cast<double *> (I_mat->x))[i * nrows + i] = 1.0;
837
838 q = SuiteSparseQR_qmult<double> (SPQR_QX, m_H, m_Htau, m_HPinv, I_mat,
839 &m_cc);
840 spqr_error_handler (&m_cc);
841
842 double *q_x = reinterpret_cast<double *> (q->x);
843 double *ret_vec = const_cast<double *> (ret.rwdata ());
844 for (octave_idx_type j = 0; j < nc; j++)
845 for (octave_idx_type i = 0; i < nrows; i++)
846 ret_vec[j * nrows + i] = q_x[j * nrows + i];
847
848 cholmod_l_free_dense (&q, &m_cc);
849 cholmod_l_free_dense (&I_mat, &m_cc);
850
851 return ret;
852
853#elif defined (HAVE_CXSPARSE)
854
855 if (econ)
856 (*current_liboctave_error_handler)
857 ("sparse-qr: economy mode with CXSparse not supported");
858
859 octave_idx_type nc = N->L->n;
860 octave_idx_type nr = nrows;
861 Matrix ret (nr, nr);
862 double *ret_vec = ret.rwdata ();
863
864 if (nr < 0 || nc < 0)
865 (*current_liboctave_error_handler) ("matrix dimension mismatch");
866
867 if (nr == 0 || nc == 0)
868 ret = Matrix (nc, nr, 0.0);
869 else
870 {
871 OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1);
872
873 for (octave_idx_type i = 0; i < nr; i++)
874 bvec[i] = 0.;
875
876 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
877
878 for (octave_idx_type j = 0, idx = 0; j < nr; j++, idx += nr)
879 {
880 octave_quit ();
881
882 bvec[j] = 1.0;
883 for (octave_idx_type i = nr; i < S->m2; i++)
884 buf[i] = 0.;
885
886 octave_idx_type nm = (nr < nc ? nr : nc);
887
888 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr);
889
890 for (octave_idx_type i = 0; i < nm; i++)
891 {
892 octave_quit ();
893
894 CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf);
895 }
896
897 for (octave_idx_type i = 0; i < nr; i++)
898 ret_vec[i+idx] = buf[i];
899
900 bvec[j] = 0.0;
901 }
902 }
903
904 return ret.transpose ();
905
906#else
907
908 octave_unused_parameter (econ);
909
910 return Matrix ();
911
912#endif
913}
914
915template <>
916template <>
917Matrix
918sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix>
919(const MArray<double>& b, octave_idx_type& info)
920{
921 info = -1;
922
923#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD) && defined (HAVE_CXSPARSE))
924
925 octave_idx_type b_nr = b.rows ();
926 octave_idx_type b_nc = b.cols ();
927 Matrix x (ncols, b_nc); // X = m_E'*(m_R\‍(Q'*B))
928
929 if (nrows < 0 || ncols < 0 || b_nc < 0 || b_nr < 0)
930 (*current_liboctave_error_handler)
931 ("matrix dimension with negative size");
932
933 if (nrows < 0 || ncols < 0 || nrows != b_nr)
934 (*current_liboctave_error_handler) ("matrix dimension mismatch");
935
936 cholmod_dense *QTB; // Q' * B
937 cholmod_dense B = ord2crd (b);
938
939 // FIXME: Process b column by column as in the CXSPARSE version below.
940 // This avoids a large dense matrix Q' * B in memory.
941 QTB = SuiteSparseQR_qmult<double> (SPQR_QTX, m_H, m_Htau, m_HPinv, &B, &m_cc);
942
943 spqr_error_handler (&m_cc);
944
945 // convert m_R into CXSPARSE matrix R2
946 CXSPARSE_DNAME (_sparse) R2;
947 R2.n = ncols;
948 R2.m = ncols;
949 R2.nzmax = m_R->nzmax;
950 R2.x = reinterpret_cast<double *> (m_R->x);
951 suitesparse_integer *R2_p;
952 suitesparse_integer *R2_i;
953 if (! suitesparse_integer_long_mismatch)
954 {
955 R2.p = reinterpret_cast<suitesparse_integer *> (m_R->p);
956 R2.i = reinterpret_cast<suitesparse_integer *> (m_R->i);
957 }
958 else
959 {
960 R2_p = new suitesparse_integer[ncols+1];
961 SuiteSparse_long *R_p = reinterpret_cast<SuiteSparse_long *> (m_R->p);
962 for (octave_idx_type i = 0; i < ncols+1; i++)
963 R2_p[i] = suitesparse_long_to_suitesparse_integer (R_p[i]);
964 R2.p = R2_p;
965 octave_idx_type nnz = from_suitesparse_long (cholmod_l_nnz (m_R, &m_cc));
966 R2_i = new suitesparse_integer[nnz];
967 SuiteSparse_long *R_i = reinterpret_cast<SuiteSparse_long *> (m_R->i);
968 for (octave_idx_type i = 0; i < nnz; i++)
969 R2_i[i] = suitesparse_long_to_suitesparse_integer (R_i[i]);
970 R2.i = R2_i;
971 }
972 R2.nz = -1;
973 double *x_vec = const_cast<double *> (x.rwdata ());
974 suitesparse_integer *E;
975 if (suitesparse_integer_long_mismatch)
976 {
977 E = new suitesparse_integer [ncols];
978 for (octave_idx_type i = 0; i < ncols; i++)
979 E[i] = suitesparse_long_to_suitesparse_integer (m_E[i]);
980 }
981 else
982 E = reinterpret_cast<suitesparse_integer *> (m_E);
983 for (octave_idx_type j = 0; j < b_nc; j++)
984 {
985 // fill x(:,j)
986 // solve (m_R\‍(Q'*B(:,j)) and store result in QTB(:,j)
987 CXSPARSE_DNAME (_usolve)
988 (&R2, &(reinterpret_cast<double *> (QTB->x)[j * b_nr]));
989 // x(:,j) = m_E' * (m_R\‍(Q'*B(:,j))
990 CXSPARSE_DNAME (_ipvec)
991 (E, &(reinterpret_cast<double *> (QTB->x)[j * b_nr]),
992 &x_vec[j * ncols], ncols);
993 }
994
995 if (suitesparse_integer_long_mismatch)
996 {
997 delete [] R2_p;
998 delete [] R2_i;
999 delete [] E;
1000 }
1001 cholmod_l_free_dense (&QTB, &m_cc);
1002
1003 info = 0;
1004
1005 return x;
1006
1007#elif defined (HAVE_CXSPARSE)
1008
1009 octave_idx_type nr = nrows;
1010 octave_idx_type nc = ncols;
1011
1012 octave_idx_type b_nc = b.cols ();
1013 octave_idx_type b_nr = b.rows ();
1014
1015 const double *bvec = b.data ();
1016
1017 Matrix x (nc, b_nc);
1018 double *vec = x.rwdata ();
1019
1020 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
1021
1022 for (octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
1023 i++, idx+=nc, bidx+=b_nr)
1024 {
1025 octave_quit ();
1026
1027 for (octave_idx_type j = nr; j < S->m2; j++)
1028 buf[j] = 0.;
1029
1030 CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
1031
1032 for (octave_idx_type j = 0; j < nc; j++)
1033 {
1034 octave_quit ();
1035
1036 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1037 }
1038
1039 CXSPARSE_DNAME (_usolve) (N->U, buf);
1040 CXSPARSE_DNAME (_ipvec) (S->q, buf, vec + idx, nc);
1041 }
1042
1043 info = 0;
1044
1045 return x;
1046
1047#else
1048
1049 octave_unused_parameter (b);
1050
1051 return Matrix ();
1052
1053#endif
1054}
1055
1056template <>
1057template <>
1058Matrix
1059sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>, Matrix>
1060(const MArray<double>& b, octave_idx_type& info) const
1061{
1062 info = -1;
1063#if defined (HAVE_CXSPARSE)
1064
1065 // These are swapped because the original matrix was transposed in
1066 // sparse_qr<SparseMatrix>::solve.
1067
1068 octave_idx_type nr = ncols;
1069 octave_idx_type nc = nrows;
1070
1071 octave_idx_type b_nc = b.cols ();
1072 octave_idx_type b_nr = b.rows ();
1073
1074 const double *bvec = b.data ();
1075
1076 Matrix x (nc, b_nc);
1077 double *vec = x.rwdata ();
1078
1079 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1080
1081 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
1082
1083 for (octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
1084 i++, idx+=nc, bidx+=b_nr)
1085 {
1086 octave_quit ();
1087
1088 for (octave_idx_type j = nr; j < nbuf; j++)
1089 buf[j] = 0.;
1090
1091 CXSPARSE_DNAME (_pvec) (S->q, bvec + bidx, buf, nr);
1092 CXSPARSE_DNAME (_utsolve) (N->U, buf);
1093
1094 for (octave_idx_type j = nr-1; j >= 0; j--)
1095 {
1096 octave_quit ();
1097
1098 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1099 }
1100
1101 CXSPARSE_DNAME (_pvec) (S->pinv, buf, vec + idx, nc);
1102 }
1103
1104 info = 0;
1105
1106 return x;
1107
1108#else
1109 octave_unused_parameter (b);
1110
1111 return Matrix ();
1112
1113#endif
1114}
1115
1116template <>
1117template <>
1118SparseMatrix
1119sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix>
1120(const SparseMatrix& b, octave_idx_type& info)
1121{
1122 info = -1;
1123
1124#if defined (HAVE_CXSPARSE)
1125
1126 octave_idx_type nr = nrows;
1127 octave_idx_type nc = ncols;
1128
1129 octave_idx_type b_nr = b.rows ();
1130 octave_idx_type b_nc = b.cols ();
1131
1132 SparseMatrix x (nc, b_nc, b.nnz ());
1133 x.xcidx (0) = 0;
1134
1135 octave_idx_type x_nz = b.nnz ();
1136 octave_idx_type ii = 0;
1137
1138 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1139 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
1140
1141 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1142 {
1143 octave_quit ();
1144
1145 for (octave_idx_type j = 0; j < b_nr; j++)
1146 Xx[j] = b.xelem (j, i);
1147
1148 for (octave_idx_type j = nr; j < S->m2; j++)
1149 buf[j] = 0.;
1150
1151 CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
1152
1153 for (octave_idx_type j = 0; j < nc; j++)
1154 {
1155 octave_quit ();
1156
1157 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1158 }
1159
1160 CXSPARSE_DNAME (_usolve) (N->U, buf);
1161 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
1162
1163 for (octave_idx_type j = 0; j < nc; j++)
1164 {
1165 double tmp = Xx[j];
1166
1167 if (tmp != 0.0)
1168 {
1169 if (ii == x_nz)
1170 {
1171 // Resize the sparse matrix
1172 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1173 sz = (sz > 10 ? sz : 10) + x_nz;
1174 x.change_capacity (sz);
1175 x_nz = sz;
1176 }
1177
1178 x.xdata (ii) = tmp;
1179 x.xridx (ii++) = j;
1180 }
1181 }
1182
1183 x.xcidx (i+1) = ii;
1184 }
1185
1186 info = 0;
1187
1188 return x;
1189
1190#else
1191
1192 octave_unused_parameter (b);
1193
1194 return SparseMatrix ();
1195
1196#endif
1197}
1198
1199template <>
1200template <>
1201SparseMatrix
1202sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseMatrix>
1203(const SparseMatrix& b, octave_idx_type& info) const
1204{
1205 info = -1;
1206
1207#if defined (HAVE_CXSPARSE)
1208
1209 // These are swapped because the original matrix was transposed in
1210 // sparse_qr<SparseMatrix>::solve.
1211
1212 octave_idx_type nr = ncols;
1213 octave_idx_type nc = nrows;
1214
1215 octave_idx_type b_nr = b.rows ();
1216 octave_idx_type b_nc = b.cols ();
1217
1218 SparseMatrix x (nc, b_nc, b.nnz ());
1219 x.xcidx (0) = 0;
1220
1221 octave_idx_type x_nz = b.nnz ();
1222 octave_idx_type ii = 0;
1223 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1224
1225 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1226 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
1227
1228 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1229 {
1230 octave_quit ();
1231
1232 for (octave_idx_type j = 0; j < b_nr; j++)
1233 Xx[j] = b.xelem (j, i);
1234
1235 for (octave_idx_type j = nr; j < nbuf; j++)
1236 buf[j] = 0.;
1237
1238 CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
1239 CXSPARSE_DNAME (_utsolve) (N->U, buf);
1240
1241 for (octave_idx_type j = nr-1; j >= 0; j--)
1242 {
1243 octave_quit ();
1244
1245 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1246 }
1247
1248 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
1249
1250 for (octave_idx_type j = 0; j < nc; j++)
1251 {
1252 double tmp = Xx[j];
1253
1254 if (tmp != 0.0)
1255 {
1256 if (ii == x_nz)
1257 {
1258 // Resize the sparse matrix
1259 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1260 sz = (sz > 10 ? sz : 10) + x_nz;
1261 x.change_capacity (sz);
1262 x_nz = sz;
1263 }
1264
1265 x.xdata (ii) = tmp;
1266 x.xridx (ii++) = j;
1267 }
1268 }
1269
1270 x.xcidx (i+1) = ii;
1271 }
1272
1273 info = 0;
1274
1275 x.maybe_compress ();
1276
1277 return x;
1278
1279#else
1280
1281 octave_unused_parameter (b);
1282
1283 return SparseMatrix ();
1284
1285#endif
1286}
1287
1288template <>
1289template <>
1290ComplexMatrix
1291sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix>
1292(const MArray<Complex>& b, octave_idx_type& info)
1293{
1294 info = -1;
1295
1296#if defined (HAVE_CXSPARSE)
1297
1298 octave_idx_type nr = nrows;
1299 octave_idx_type nc = ncols;
1300
1301 octave_idx_type b_nc = b.cols ();
1302 octave_idx_type b_nr = b.rows ();
1303
1304 ComplexMatrix x (nc, b_nc);
1305 Complex *vec = x.rwdata ();
1306
1307 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1308 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1309 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
1310
1311 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1312 {
1313 octave_quit ();
1314
1315 for (octave_idx_type j = 0; j < b_nr; j++)
1316 {
1317 Complex c = b.xelem (j, i);
1318 Xx[j] = c.real ();
1319 Xz[j] = c.imag ();
1320 }
1321
1322 for (octave_idx_type j = nr; j < S->m2; j++)
1323 buf[j] = 0.;
1324
1325 CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
1326
1327 for (octave_idx_type j = 0; j < nc; j++)
1328 {
1329 octave_quit ();
1330
1331 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1332 }
1333
1334 CXSPARSE_DNAME (_usolve) (N->U, buf);
1335 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
1336
1337 for (octave_idx_type j = nr; j < S->m2; j++)
1338 buf[j] = 0.;
1339
1340 CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
1341
1342 for (octave_idx_type j = 0; j < nc; j++)
1343 {
1344 octave_quit ();
1345
1346 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1347 }
1348
1349 CXSPARSE_DNAME (_usolve) (N->U, buf);
1350 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
1351
1352 for (octave_idx_type j = 0; j < nc; j++)
1353 vec[j+idx] = Complex (Xx[j], Xz[j]);
1354 }
1355
1356 info = 0;
1357
1358 return x;
1359
1360#else
1361
1362 octave_unused_parameter (b);
1363
1364 return ComplexMatrix ();
1365
1366#endif
1367}
1368
1369template <>
1370template <>
1371ComplexMatrix
1372sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix>
1373(const MArray<Complex>& b, octave_idx_type& info) const
1374{
1375 info = -1;
1376
1377#if defined (HAVE_CXSPARSE)
1378
1379 // These are swapped because the original matrix was transposed in
1380 // sparse_qr<SparseMatrix>::solve.
1381
1382 octave_idx_type nr = ncols;
1383 octave_idx_type nc = nrows;
1384
1385 octave_idx_type b_nc = b.cols ();
1386 octave_idx_type b_nr = b.rows ();
1387
1388 ComplexMatrix x (nc, b_nc);
1389 Complex *vec = x.rwdata ();
1390
1391 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1392
1393 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1394 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1395 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
1396
1397 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1398 {
1399 octave_quit ();
1400
1401 for (octave_idx_type j = 0; j < b_nr; j++)
1402 {
1403 Complex c = b.xelem (j, i);
1404 Xx[j] = c.real ();
1405 Xz[j] = c.imag ();
1406 }
1407
1408 for (octave_idx_type j = nr; j < nbuf; j++)
1409 buf[j] = 0.;
1410
1411 CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
1412 CXSPARSE_DNAME (_utsolve) (N->U, buf);
1413
1414 for (octave_idx_type j = nr-1; j >= 0; j--)
1415 {
1416 octave_quit ();
1417
1418 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1419 }
1420
1421 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
1422
1423 for (octave_idx_type j = nr; j < nbuf; j++)
1424 buf[j] = 0.;
1425
1426 CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
1427 CXSPARSE_DNAME (_utsolve) (N->U, buf);
1428
1429 for (octave_idx_type j = nr-1; j >= 0; j--)
1430 {
1431 octave_quit ();
1432
1433 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1434 }
1435
1436 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
1437
1438 for (octave_idx_type j = 0; j < nc; j++)
1439 vec[j+idx] = Complex (Xx[j], Xz[j]);
1440 }
1441
1442 info = 0;
1443
1444 return x;
1445
1446#else
1447
1448 octave_unused_parameter (b);
1449
1450 return ComplexMatrix ();
1451
1452#endif
1453}
1454
1455// Complex-valued matrices.
1456
1457template <>
1458sparse_qr<SparseComplexMatrix>::sparse_qr_rep::sparse_qr_rep
1459(const SparseComplexMatrix& a, int order)
1460 : nrows (a.rows ()), ncols (a.columns ())
1461#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
1462 , m_cc (), m_R (nullptr), m_E (nullptr), m_H (nullptr),
1463 m_Htau (nullptr), m_HPinv (nullptr)
1464{
1465 octave_idx_type nr = a.rows ();
1466 octave_idx_type nc = a.cols ();
1467
1468 if (nr < 0 || nc < 0)
1469 (*current_liboctave_error_handler)
1470 ("matrix dimension with negative size");
1471
1472 if (order < 0 || order > 9)
1473 (*current_liboctave_error_handler)
1474 ("ordering %d is not supported by SPQR", order);
1475
1476 cholmod_l_start (&m_cc);
1477 cholmod_sparse A = ocs2ccs (a);
1478
1479 SuiteSparseQR<Complex> (order, static_cast<double> (SPQR_DEFAULT_TOL),
1480 static_cast<SuiteSparse_long> (A.nrow),
1481 &A, &m_R, &m_E, &m_H,
1482 &m_HPinv, &m_Htau, &m_cc);
1483 spqr_error_handler (&m_cc);
1484
1485 if (octave_suitesparse_ptr_size_mismatch)
1486 {
1487 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
1488 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
1489 }
1490}
1491
1492#elif defined (HAVE_CXSPARSE)
1493 , S (nullptr), N (nullptr)
1494{
1495 CXSPARSE_ZNAME () A;
1496
1497 A.nzmax = a.nnz ();
1498 A.m = nrows;
1499 A.n = ncols;
1500 // Cast away const on A, with full knowledge that CSparse won't touch it
1501 // Prevents the methods below making a copy of the data.
1502 A.p = const_cast<suitesparse_integer *>
1503 (to_suitesparse_intptr (a.cidx ()));
1504 A.i = const_cast<suitesparse_integer *>
1505 (to_suitesparse_intptr (a.ridx ()));
1506 A.x = const_cast<cs_complex_t *>
1507 (reinterpret_cast<const cs_complex_t *> (a.data ()));
1508 A.nz = -1;
1509
1510 S = CXSPARSE_ZNAME (_sqr) (order, &A, 1);
1511 N = CXSPARSE_ZNAME (_qr) (&A, S);
1512
1513 if (! N)
1514 (*current_liboctave_error_handler)
1515 ("sparse_qr: sparse matrix QR factorization filled");
1516
1517}
1518
1519#else
1520
1521{
1522 octave_unused_parameter (order);
1523
1524 (*current_liboctave_error_handler)
1525 ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
1526}
1527
1528#endif
1529
1530template <>
1531sparse_qr<SparseComplexMatrix>::sparse_qr_rep::~sparse_qr_rep ()
1532{
1533#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
1534
1535 cholmod_l_free_sparse (&m_R, &m_cc);
1536 cholmod_l_free_sparse (&m_H, &m_cc);
1537 cholmod_l_free_dense (&m_Htau, &m_cc);
1538 free (m_E); // FIXME: use cholmod_l_free
1539 free (m_HPinv);
1540 cholmod_l_finish (&m_cc);
1541
1542#elif defined (HAVE_CXSPARSE)
1543
1544 CXSPARSE_ZNAME (_sfree) (S);
1545 CXSPARSE_ZNAME (_nfree) (N);
1546
1547#endif
1548}
1549
1550template <>
1551SparseComplexMatrix
1552sparse_qr<SparseComplexMatrix>::sparse_qr_rep::V () const
1553{
1554#if defined (HAVE_CXSPARSE)
1555 // Drop zeros from V and sort
1556 // FIXME: Is the double transpose to sort necessary?
1557
1558 CXSPARSE_ZNAME (_dropzeros) (N->L);
1559 CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1);
1560 CXSPARSE_ZNAME (_spfree) (N->L);
1561 N->L = CXSPARSE_ZNAME (_transpose) (D, 1);
1562 CXSPARSE_ZNAME (_spfree) (D);
1563
1564 octave_idx_type nc = N->L->n;
1565 octave_idx_type nz = N->L->p[nc];
1566 SparseComplexMatrix ret (N->L->m, nc, nz);
1567
1568 for (octave_idx_type j = 0; j < nc+1; j++)
1569 ret.xcidx (j) = N->L->p[j];
1570
1571 for (octave_idx_type j = 0; j < nz; j++)
1572 {
1573 ret.xridx (j) = N->L->i[j];
1574 ret.xdata (j) = reinterpret_cast<Complex *> (N->L->x)[j];
1575 }
1576
1577 return ret;
1578
1579#else
1580
1581 return SparseComplexMatrix ();
1582
1583#endif
1584}
1585
1586template <>
1587SparseComplexMatrix
1588sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const
1589{
1590#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
1591
1592 octave_idx_type nr = from_size_t (m_R->nrow);
1593 octave_idx_type nc = from_size_t (m_R->ncol);
1594 octave_idx_type nz
1595 = from_suitesparse_long
1596 (cholmod_l_nnz (m_R, const_cast<cholmod_common *> (&m_cc)));
1597
1598 // FIXME: Does this work if econ = true?
1599 SparseComplexMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz);
1600 SuiteSparse_long *Rp = reinterpret_cast<SuiteSparse_long *> (m_R->p);
1601 SuiteSparse_long *Ri = reinterpret_cast<SuiteSparse_long *> (m_R->i);
1602
1603 for (octave_idx_type j = 0; j < nc + 1; j++)
1604 ret.xcidx (j) = from_suitesparse_long (Rp[j]);
1605
1606 for (octave_idx_type j = 0; j < nz; j++)
1607 {
1608 ret.xridx (j) = from_suitesparse_long (Ri[j]);
1609 ret.xdata (j) = (reinterpret_cast<Complex *> (m_R->x))[j];
1610 }
1611
1612 return ret;
1613
1614#elif defined (HAVE_CXSPARSE)
1615
1616 // Drop zeros from R and sort
1617 // FIXME: Is the double transpose to sort necessary?
1618
1619 CXSPARSE_ZNAME (_dropzeros) (N->U);
1620 CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1);
1621 CXSPARSE_ZNAME (_spfree) (N->U);
1622 N->U = CXSPARSE_ZNAME (_transpose) (D, 1);
1623 CXSPARSE_ZNAME (_spfree) (D);
1624
1625 octave_idx_type nc = N->U->n;
1626 octave_idx_type nz = N->U->p[nc];
1627
1628 SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows),
1629 nc, nz);
1630
1631 for (octave_idx_type j = 0; j < nc+1; j++)
1632 ret.xcidx (j) = N->U->p[j];
1633
1634 for (octave_idx_type j = 0; j < nz; j++)
1635 {
1636 ret.xridx (j) = N->U->i[j];
1637 ret.xdata (j) = reinterpret_cast<Complex *> (N->U->x)[j];
1638 }
1639
1640 return ret;
1641
1642#else
1643
1644 octave_unused_parameter (econ);
1645
1646 return SparseComplexMatrix ();
1647
1648#endif
1649}
1650
1651template <>
1652ComplexMatrix
1653sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C
1654(const ComplexMatrix& b, bool econ)
1655{
1656#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
1657
1658 // FIXME: not tested
1659 octave_idx_type nr = (econ
1660 ? (ncols > nrows ? nrows : ncols)
1661 : nrows);
1662 octave_idx_type b_nr = b.rows ();
1663 octave_idx_type b_nc = b.cols ();
1664 ComplexMatrix ret (nr, b_nc);
1665
1666 if (nrows != b_nr)
1667 (*current_liboctave_error_handler) ("matrix dimension mismatch");
1668
1669 if (b_nc < 0 || b_nr < 0)
1670 (*current_liboctave_error_handler)
1671 ("matrix dimension with negative size");
1672
1673 cholmod_dense *QTB; // Q' * B
1674 cholmod_dense B = ocd2ccd (b);
1675
1676 QTB = SuiteSparseQR_qmult<Complex> (SPQR_QTX, m_H, m_Htau, m_HPinv, &B,
1677 &m_cc);
1678 spqr_error_handler (&m_cc);
1679
1680 // copy QTB into ret
1681 Complex *QTB_x = reinterpret_cast<Complex *> (QTB->x);
1682 Complex *ret_vec = reinterpret_cast<Complex *> (ret.rwdata ());
1683 for (octave_idx_type j = 0; j < b_nc; j++)
1684 for (octave_idx_type i = 0; i < nr; i++)
1685 ret_vec[j * nr + i] = QTB_x[j * b_nr + i];
1686
1687 cholmod_l_free_dense (&QTB, &m_cc);
1688
1689 return ret;
1690
1691#elif defined (HAVE_CXSPARSE)
1692
1693 if (econ)
1694 (*current_liboctave_error_handler)
1695 ("sparse-qr: economy mode with CXSparse not supported");
1696
1697 octave_idx_type b_nr = b.rows ();
1698 octave_idx_type b_nc = b.cols ();
1699 octave_idx_type nc = N->L->n;
1700 octave_idx_type nr = nrows;
1701 const cs_complex_t *bvec
1702 = reinterpret_cast<const cs_complex_t *> (b.data ());
1703 ComplexMatrix ret (b_nr, b_nc);
1704 Complex *vec = ret.rwdata ();
1705
1706 if (nr < 0 || nc < 0 || nr != b_nr)
1707 (*current_liboctave_error_handler) ("matrix dimension mismatch");
1708
1709 if (nr == 0 || nc == 0 || b_nc == 0)
1710 ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0));
1711 else
1712 {
1713 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
1714
1715 for (octave_idx_type j = 0, idx = 0;
1716 j < b_nc;
1717 j++, idx += b_nr)
1718 {
1719 octave_quit ();
1720
1721 octave_idx_type nm = (nr < nc ? nr : nc);
1722
1723 CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + idx,
1724 reinterpret_cast<cs_complex_t *> (buf),
1725 b_nr);
1726
1727 for (octave_idx_type i = 0; i < nm; i++)
1728 {
1729 octave_quit ();
1730
1731 CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i],
1732 reinterpret_cast<cs_complex_t *> (buf));
1733 }
1734
1735 for (octave_idx_type i = 0; i < b_nr; i++)
1736 vec[i+idx] = buf[i];
1737 }
1738 }
1739
1740 return ret;
1741
1742#else
1743
1744 octave_unused_parameter (b);
1745 octave_unused_parameter (econ);
1746
1747 return ComplexMatrix ();
1748
1749#endif
1750}
1751
1752template <>
1753ComplexMatrix
1754sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (bool econ)
1755{
1756#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
1757
1758 octave_idx_type nc = (econ
1759 ? (ncols > nrows ? nrows : ncols)
1760 : nrows);
1761 ComplexMatrix ret (nrows, nc);
1762 cholmod_dense *q;
1763
1764 // I is nrows x nrows identity matrix
1765 cholmod_dense *I_mat
1766 = reinterpret_cast<cholmod_dense *>
1767 (cholmod_l_allocate_dense (nrows, nrows, nrows, CHOLMOD_COMPLEX, &m_cc));
1768
1769 for (octave_idx_type i = 0; i < nrows * nrows; i++)
1770 (reinterpret_cast<Complex *> (I_mat->x))[i] = 0.0;
1771
1772 for (octave_idx_type i = 0; i < nrows; i++)
1773 (reinterpret_cast<Complex *> (I_mat->x))[i * nrows + i] = 1.0;
1774
1775 q = SuiteSparseQR_qmult<Complex> (SPQR_QX, m_H, m_Htau, m_HPinv, I_mat,
1776 &m_cc);
1777 spqr_error_handler (&m_cc);
1778
1779 Complex *q_x = reinterpret_cast<Complex *> (q->x);
1780 Complex *ret_vec = const_cast<Complex *> (ret.rwdata ());
1781
1782 for (octave_idx_type j = 0; j < nc; j++)
1783 for (octave_idx_type i = 0; i < nrows; i++)
1784 ret_vec[j * nrows + i] = q_x[j * nrows + i];
1785
1786 cholmod_l_free_dense (&q, &m_cc);
1787 cholmod_l_free_dense (&I_mat, &m_cc);
1788
1789 return ret;
1790
1791#elif defined (HAVE_CXSPARSE)
1792
1793 if (econ)
1794 (*current_liboctave_error_handler)
1795 ("sparse-qr: economy mode with CXSparse not supported");
1796
1797 octave_idx_type nc = N->L->n;
1798 octave_idx_type nr = nrows;
1799 ComplexMatrix ret (nr, nr);
1800 Complex *vec = ret.rwdata ();
1801
1802 if (nr < 0 || nc < 0)
1803 (*current_liboctave_error_handler) ("matrix dimension mismatch");
1804
1805 if (nr == 0 || nc == 0)
1806 ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0));
1807 else
1808 {
1809 OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr);
1810
1811 for (octave_idx_type i = 0; i < nr; i++)
1812 bvec[i] = 0.0;
1813
1814 OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2);
1815
1816 for (octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr)
1817 {
1818 octave_quit ();
1819
1820 bvec[j] = 1.0;
1821
1822 octave_idx_type nm = (nr < nc ? nr : nc);
1823
1824 CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec,
1825 reinterpret_cast<cs_complex_t *> (buf),
1826 nr);
1827
1828 for (octave_idx_type i = 0; i < nm; i++)
1829 {
1830 octave_quit ();
1831
1832 CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i],
1833 reinterpret_cast<cs_complex_t *> (buf));
1834 }
1835
1836 for (octave_idx_type i = 0; i < nr; i++)
1837 vec[i+idx] = buf[i];
1838
1839 bvec[j] = 0.0;
1840 }
1841 }
1842
1843 return ret.hermitian ();
1844
1845#else
1846
1847 octave_unused_parameter (econ);
1848
1849 return ComplexMatrix ();
1850
1851#endif
1852}
1853
1854template <>
1855template <>
1856SparseComplexMatrix
1857sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
1858 SparseComplexMatrix>
1859 (const SparseComplexMatrix& b, octave_idx_type& info)
1860{
1861 info = -1;
1862
1863#if defined (HAVE_CXSPARSE)
1864
1865 octave_idx_type nr = nrows;
1866 octave_idx_type nc = ncols;
1867
1868 octave_idx_type b_nr = b.rows ();
1869 octave_idx_type b_nc = b.cols ();
1870
1871 SparseComplexMatrix x (nc, b_nc, b.nnz ());
1872 x.xcidx (0) = 0;
1873
1874 octave_idx_type x_nz = b.nnz ();
1875 octave_idx_type ii = 0;
1876
1877 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1878 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1879 OCTAVE_LOCAL_BUFFER (double, buf, S->m2);
1880
1881 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1882 {
1883 octave_quit ();
1884
1885 for (octave_idx_type j = 0; j < b_nr; j++)
1886 {
1887 Complex c = b.xelem (j, i);
1888 Xx[j] = c.real ();
1889 Xz[j] = c.imag ();
1890 }
1891
1892 for (octave_idx_type j = nr; j < S->m2; j++)
1893 buf[j] = 0.;
1894
1895 CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr);
1896
1897 for (octave_idx_type j = 0; j < nc; j++)
1898 {
1899 octave_quit ();
1900
1901 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1902 }
1903
1904 CXSPARSE_DNAME (_usolve) (N->U, buf);
1905 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc);
1906
1907 for (octave_idx_type j = nr; j < S->m2; j++)
1908 buf[j] = 0.;
1909
1910 CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr);
1911
1912 for (octave_idx_type j = 0; j < nc; j++)
1913 {
1914 octave_quit ();
1915
1916 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
1917 }
1918
1919 CXSPARSE_DNAME (_usolve) (N->U, buf);
1920 CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc);
1921
1922 for (octave_idx_type j = 0; j < nc; j++)
1923 {
1924 Complex tmp = Complex (Xx[j], Xz[j]);
1925
1926 if (tmp != 0.0)
1927 {
1928 if (ii == x_nz)
1929 {
1930 // Resize the sparse matrix
1931 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
1932 sz = (sz > 10 ? sz : 10) + x_nz;
1933 x.change_capacity (sz);
1934 x_nz = sz;
1935 }
1936
1937 x.xdata (ii) = tmp;
1938 x.xridx (ii++) = j;
1939 }
1940 }
1941
1942 x.xcidx (i+1) = ii;
1943 }
1944
1945 info = 0;
1946
1947 return x;
1948
1949#else
1950
1951 octave_unused_parameter (b);
1952
1953 return SparseComplexMatrix ();
1954
1955#endif
1956}
1957
1958template <>
1959template <>
1960SparseComplexMatrix
1961sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix,
1962 SparseComplexMatrix>
1963 (const SparseComplexMatrix& b, octave_idx_type& info) const
1964{
1965 info = -1;
1966
1967#if defined (HAVE_CXSPARSE)
1968
1969 // These are swapped because the original matrix was transposed in
1970 // sparse_qr<SparseMatrix>::solve.
1971
1972 octave_idx_type nr = ncols;
1973 octave_idx_type nc = nrows;
1974
1975 octave_idx_type b_nr = b.rows ();
1976 octave_idx_type b_nc = b.cols ();
1977
1978 SparseComplexMatrix x (nc, b_nc, b.nnz ());
1979 x.xcidx (0) = 0;
1980
1981 octave_idx_type x_nz = b.nnz ();
1982 octave_idx_type ii = 0;
1983 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
1984
1985 OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc));
1986 OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc));
1987 OCTAVE_LOCAL_BUFFER (double, buf, nbuf);
1988
1989 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
1990 {
1991 octave_quit ();
1992
1993 for (octave_idx_type j = 0; j < b_nr; j++)
1994 {
1995 Complex c = b.xelem (j, i);
1996 Xx[j] = c.real ();
1997 Xz[j] = c.imag ();
1998 }
1999
2000 for (octave_idx_type j = nr; j < nbuf; j++)
2001 buf[j] = 0.;
2002
2003 CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr);
2004 CXSPARSE_DNAME (_utsolve) (N->U, buf);
2005
2006 for (octave_idx_type j = nr-1; j >= 0; j--)
2007 {
2008 octave_quit ();
2009
2010 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
2011 }
2012
2013 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc);
2014
2015 for (octave_idx_type j = nr; j < nbuf; j++)
2016 buf[j] = 0.;
2017
2018 CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr);
2019 CXSPARSE_DNAME (_utsolve) (N->U, buf);
2020
2021 for (octave_idx_type j = nr-1; j >= 0; j--)
2022 {
2023 octave_quit ();
2024
2025 CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf);
2026 }
2027
2028 CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc);
2029
2030 for (octave_idx_type j = 0; j < nc; j++)
2031 {
2032 Complex tmp = Complex (Xx[j], Xz[j]);
2033
2034 if (tmp != 0.0)
2035 {
2036 if (ii == x_nz)
2037 {
2038 // Resize the sparse matrix
2039 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2040 sz = (sz > 10 ? sz : 10) + x_nz;
2041 x.change_capacity (sz);
2042 x_nz = sz;
2043 }
2044
2045 x.xdata (ii) = tmp;
2046 x.xridx (ii++) = j;
2047 }
2048 }
2049
2050 x.xcidx (i+1) = ii;
2051 }
2052
2053 info = 0;
2054
2055 x.maybe_compress ();
2056
2057 return x;
2058
2059#else
2060
2061 octave_unused_parameter (b);
2062
2063 return SparseComplexMatrix ();
2064
2065#endif
2066}
2067
2068template <>
2069template <>
2070ComplexMatrix
2071sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>,
2072 ComplexMatrix>
2073 (const MArray<double>& b, octave_idx_type& info)
2074{
2075 info = -1;
2076
2077#if defined (HAVE_CXSPARSE)
2078
2079 octave_idx_type nr = nrows;
2080 octave_idx_type nc = ncols;
2081
2082 octave_idx_type b_nc = b.cols ();
2083 octave_idx_type b_nr = b.rows ();
2084
2085 ComplexMatrix x (nc, b_nc);
2086 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.rwdata ());
2087
2088 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
2089 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
2090
2091 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2092 {
2093 octave_quit ();
2094
2095 for (octave_idx_type j = 0; j < b_nr; j++)
2096 Xx[j] = b.xelem (j, i);
2097
2098 for (octave_idx_type j = nr; j < S->m2; j++)
2099 buf[j] = 0.0;
2100
2101 CXSPARSE_ZNAME (_ipvec) (S->pinv,
2102 reinterpret_cast<cs_complex_t *> (Xx),
2103 buf, nr);
2104
2105 for (octave_idx_type j = 0; j < nc; j++)
2106 {
2107 octave_quit ();
2108
2109 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
2110 }
2111
2112 CXSPARSE_ZNAME (_usolve) (N->U, buf);
2113 CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
2114 }
2115
2116 info = 0;
2117
2118 return x;
2119
2120#else
2121
2122 octave_unused_parameter (b);
2123
2124 return ComplexMatrix ();
2125
2126#endif
2127}
2128
2129template <>
2130template <>
2131ComplexMatrix
2132sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<double>,
2133 ComplexMatrix>
2134 (const MArray<double>& b, octave_idx_type& info) const
2135{
2136 info = -1;
2137
2138#if defined (HAVE_CXSPARSE)
2139
2140 // These are swapped because the original matrix was transposed in
2141 // sparse_qr<SparseComplexMatrix>::solve.
2142
2143 octave_idx_type nr = ncols;
2144 octave_idx_type nc = nrows;
2145
2146 octave_idx_type b_nc = b.cols ();
2147 octave_idx_type b_nr = b.rows ();
2148
2149 ComplexMatrix x (nc, b_nc);
2150 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.rwdata ());
2151
2152 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
2153
2154 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
2155 OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr);
2156 OCTAVE_LOCAL_BUFFER (double, B, nr);
2157
2158 for (octave_idx_type i = 0; i < nr; i++)
2159 B[i] = N->B[i];
2160
2161 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2162 {
2163 octave_quit ();
2164
2165 for (octave_idx_type j = 0; j < b_nr; j++)
2166 Xx[j] = b.xelem (j, i);
2167
2168 for (octave_idx_type j = nr; j < nbuf; j++)
2169 buf[j] = 0.0;
2170
2171 CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *> (Xx),
2172 buf, nr);
2173 CXSPARSE_ZNAME (_utsolve) (N->U, buf);
2174
2175 for (octave_idx_type j = nr-1; j >= 0; j--)
2176 {
2177 octave_quit ();
2178
2179 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
2180 }
2181
2182 CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
2183 }
2184
2185 info = 0;
2186
2187 return x;
2188
2189#else
2190
2191 octave_unused_parameter (b);
2192
2193 return ComplexMatrix ();
2194
2195#endif
2196}
2197
2198template <>
2199template <>
2200SparseComplexMatrix
2201sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix,
2202 SparseComplexMatrix>
2203 (const SparseMatrix& b, octave_idx_type& info)
2204{
2205 info = -1;
2206
2207#if defined (HAVE_CXSPARSE)
2208
2209 octave_idx_type nr = nrows;
2210 octave_idx_type nc = ncols;
2211
2212 octave_idx_type b_nc = b.cols ();
2213 octave_idx_type b_nr = b.rows ();
2214
2215 SparseComplexMatrix x (nc, b_nc, b.nnz ());
2216 x.xcidx (0) = 0;
2217
2218 octave_idx_type x_nz = b.nnz ();
2219 octave_idx_type ii = 0;
2220
2221 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
2222 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
2223
2224 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2225 {
2226 octave_quit ();
2227
2228 for (octave_idx_type j = 0; j < b_nr; j++)
2229 Xx[j] = b.xelem (j, i);
2230
2231 for (octave_idx_type j = nr; j < S->m2; j++)
2232 buf[j] = 0.0;
2233
2234 CXSPARSE_ZNAME (_ipvec) (S->pinv,
2235 reinterpret_cast<cs_complex_t *> (Xx),
2236 buf, nr);
2237
2238 for (octave_idx_type j = 0; j < nc; j++)
2239 {
2240 octave_quit ();
2241
2242 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
2243 }
2244
2245 CXSPARSE_ZNAME (_usolve) (N->U, buf);
2246 CXSPARSE_ZNAME (_ipvec) (S->q, buf,
2247 reinterpret_cast<cs_complex_t *> (Xx),
2248 nc);
2249
2250 for (octave_idx_type j = 0; j < nc; j++)
2251 {
2252 Complex tmp = Xx[j];
2253
2254 if (tmp != 0.0)
2255 {
2256 if (ii == x_nz)
2257 {
2258 // Resize the sparse matrix
2259 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2260 sz = (sz > 10 ? sz : 10) + x_nz;
2261 x.change_capacity (sz);
2262 x_nz = sz;
2263 }
2264
2265 x.xdata (ii) = tmp;
2266 x.xridx (ii++) = j;
2267 }
2268 }
2269
2270 x.xcidx (i+1) = ii;
2271 }
2272
2273 info = 0;
2274
2275 x.maybe_compress ();
2276
2277 return x;
2278
2279#else
2280
2281 octave_unused_parameter (b);
2282
2283 return SparseComplexMatrix ();
2284
2285#endif
2286}
2287
2288template <>
2289template <>
2290SparseComplexMatrix
2291sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseMatrix,
2292 SparseComplexMatrix>
2293 (const SparseMatrix& b, octave_idx_type& info) const
2294{
2295 info = -1;
2296
2297#if defined (HAVE_CXSPARSE)
2298
2299 // These are swapped because the original matrix was transposed in
2300 // sparse_qr<SparseComplexMatrix>::solve.
2301
2302 octave_idx_type nr = ncols;
2303 octave_idx_type nc = nrows;
2304
2305 octave_idx_type b_nc = b.cols ();
2306 octave_idx_type b_nr = b.rows ();
2307
2308 SparseComplexMatrix x (nc, b_nc, b.nnz ());
2309 x.xcidx (0) = 0;
2310
2311 octave_idx_type x_nz = b.nnz ();
2312 octave_idx_type ii = 0;
2313 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
2314
2315 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
2316 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
2317 OCTAVE_LOCAL_BUFFER (double, B, nr);
2318
2319 for (octave_idx_type i = 0; i < nr; i++)
2320 B[i] = N->B[i];
2321
2322 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2323 {
2324 octave_quit ();
2325
2326 for (octave_idx_type j = 0; j < b_nr; j++)
2327 Xx[j] = b.xelem (j, i);
2328
2329 for (octave_idx_type j = nr; j < nbuf; j++)
2330 buf[j] = 0.0;
2331
2332 CXSPARSE_ZNAME (_pvec) (S->q,
2333 reinterpret_cast<cs_complex_t *> (Xx),
2334 buf, nr);
2335 CXSPARSE_ZNAME (_utsolve) (N->U, buf);
2336
2337 for (octave_idx_type j = nr-1; j >= 0; j--)
2338 {
2339 octave_quit ();
2340
2341 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
2342 }
2343
2344 CXSPARSE_ZNAME (_pvec) (S->pinv, buf,
2345 reinterpret_cast<cs_complex_t *> (Xx),
2346 nc);
2347
2348 for (octave_idx_type j = 0; j < nc; j++)
2349 {
2350 Complex tmp = Xx[j];
2351
2352 if (tmp != 0.0)
2353 {
2354 if (ii == x_nz)
2355 {
2356 // Resize the sparse matrix
2357 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2358 sz = (sz > 10 ? sz : 10) + x_nz;
2359 x.change_capacity (sz);
2360 x_nz = sz;
2361 }
2362
2363 x.xdata (ii) = tmp;
2364 x.xridx (ii++) = j;
2365 }
2366 }
2367
2368 x.xcidx (i+1) = ii;
2369 }
2370
2371 info = 0;
2372
2373 x.maybe_compress ();
2374
2375 return x;
2376
2377#else
2378
2379 octave_unused_parameter (b);
2380
2381 return SparseComplexMatrix ();
2382
2383#endif
2384}
2385
2386template <>
2387template <>
2388ComplexMatrix
2389sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>,
2390 ComplexMatrix>
2391 (const MArray<Complex>& b, octave_idx_type& info)
2392{
2393 info = -1;
2394
2395#if defined (HAVE_CXSPARSE)
2396
2397 octave_idx_type nr = nrows;
2398 octave_idx_type nc = ncols;
2399
2400 octave_idx_type b_nc = b.cols ();
2401 octave_idx_type b_nr = b.rows ();
2402
2403 const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>
2404 (b.data ());
2405
2406 ComplexMatrix x (nc, b_nc);
2407 cs_complex_t *vec = reinterpret_cast<cs_complex_t *>
2408 (x.rwdata ());
2409
2410 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
2411
2412 for (octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
2413 i++, idx+=nc, bidx+=b_nr)
2414 {
2415 octave_quit ();
2416
2417 for (octave_idx_type j = nr; j < S->m2; j++)
2418 buf[j] = 0.0;
2419
2420 CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr);
2421
2422 for (octave_idx_type j = 0; j < nc; j++)
2423 {
2424 octave_quit ();
2425
2426 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
2427 }
2428
2429 CXSPARSE_ZNAME (_usolve) (N->U, buf);
2430 CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc);
2431 }
2432
2433 info = 0;
2434
2435 return x;
2436
2437#else
2438
2439 octave_unused_parameter (b);
2440
2441 return ComplexMatrix ();
2442
2443#endif
2444}
2445
2446template <>
2447template <>
2448ComplexMatrix
2449sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>,
2450 ComplexMatrix>
2451 (const MArray<Complex>& b, octave_idx_type& info) const
2452{
2453 info = -1;
2454
2455#if defined (HAVE_CXSPARSE)
2456
2457 // These are swapped because the original matrix was transposed in
2458 // sparse_qr<SparseComplexMatrix>::solve.
2459
2460 octave_idx_type nr = ncols;
2461 octave_idx_type nc = nrows;
2462
2463 octave_idx_type b_nc = b.cols ();
2464 octave_idx_type b_nr = b.rows ();
2465
2466 const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>
2467 (b.data ());
2468
2469 ComplexMatrix x (nc, b_nc);
2470 cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.rwdata ());
2471
2472 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
2473
2474 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
2475 OCTAVE_LOCAL_BUFFER (double, B, nr);
2476
2477 for (octave_idx_type i = 0; i < nr; i++)
2478 B[i] = N->B[i];
2479
2480 for (octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
2481 i++, idx+=nc, bidx+=b_nr)
2482 {
2483 octave_quit ();
2484
2485 for (octave_idx_type j = nr; j < nbuf; j++)
2486 buf[j] = 0.0;
2487
2488 CXSPARSE_ZNAME (_pvec) (S->q, bvec + bidx, buf, nr);
2489 CXSPARSE_ZNAME (_utsolve) (N->U, buf);
2490
2491 for (octave_idx_type j = nr-1; j >= 0; j--)
2492 {
2493 octave_quit ();
2494
2495 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
2496 }
2497
2498 CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc);
2499 }
2500
2501 info = 0;
2502
2503 return x;
2504
2505#else
2506
2507 octave_unused_parameter (b);
2508
2509 return ComplexMatrix ();
2510
2511#endif
2512}
2513
2514template <>
2515template <>
2516SparseComplexMatrix
2517sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
2518 SparseComplexMatrix>
2519 (const SparseComplexMatrix& b, octave_idx_type& info)
2520{
2521 info = -1;
2522
2523#if defined (HAVE_CXSPARSE)
2524
2525 octave_idx_type nr = nrows;
2526 octave_idx_type nc = ncols;
2527
2528 octave_idx_type b_nc = b.cols ();
2529 octave_idx_type b_nr = b.rows ();
2530
2531 SparseComplexMatrix x (nc, b_nc, b.nnz ());
2532 x.xcidx (0) = 0;
2533
2534 octave_idx_type x_nz = b.nnz ();
2535 octave_idx_type ii = 0;
2536
2537 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
2538 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2);
2539
2540 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2541 {
2542 octave_quit ();
2543
2544 for (octave_idx_type j = 0; j < b_nr; j++)
2545 Xx[j] = b.xelem (j, i);
2546
2547 for (octave_idx_type j = nr; j < S->m2; j++)
2548 buf[j] = 0.0;
2549
2550 CXSPARSE_ZNAME (_ipvec) (S->pinv,
2551 reinterpret_cast<cs_complex_t *> (Xx),
2552 buf, nr);
2553
2554 for (octave_idx_type j = 0; j < nc; j++)
2555 {
2556 octave_quit ();
2557
2558 CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf);
2559 }
2560
2561 CXSPARSE_ZNAME (_usolve) (N->U, buf);
2562 CXSPARSE_ZNAME (_ipvec) (S->q, buf,
2563 reinterpret_cast<cs_complex_t *> (Xx),
2564 nc);
2565
2566 for (octave_idx_type j = 0; j < nc; j++)
2567 {
2568 Complex tmp = Xx[j];
2569
2570 if (tmp != 0.0)
2571 {
2572 if (ii == x_nz)
2573 {
2574 // Resize the sparse matrix
2575 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2576 sz = (sz > 10 ? sz : 10) + x_nz;
2577 x.change_capacity (sz);
2578 x_nz = sz;
2579 }
2580
2581 x.xdata (ii) = tmp;
2582 x.xridx (ii++) = j;
2583 }
2584 }
2585
2586 x.xcidx (i+1) = ii;
2587 }
2588
2589 info = 0;
2590
2591 x.maybe_compress ();
2592
2593 return x;
2594
2595#else
2596
2597 octave_unused_parameter (b);
2598
2599 return SparseComplexMatrix ();
2600
2601#endif
2602}
2603
2604template <>
2605template <>
2606SparseComplexMatrix
2607sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix,
2608 SparseComplexMatrix>
2609 (const SparseComplexMatrix& b, octave_idx_type& info) const
2610{
2611 info = -1;
2612
2613#if defined (HAVE_CXSPARSE)
2614
2615 // These are swapped because the original matrix was transposed in
2616 // sparse_qr<SparseComplexMatrix>::solve.
2617
2618 octave_idx_type nr = ncols;
2619 octave_idx_type nc = nrows;
2620
2621 octave_idx_type b_nc = b.cols ();
2622 octave_idx_type b_nr = b.rows ();
2623
2624 SparseComplexMatrix x (nc, b_nc, b.nnz ());
2625 x.xcidx (0) = 0;
2626
2627 octave_idx_type x_nz = b.nnz ();
2628 octave_idx_type ii = 0;
2629 octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2);
2630
2631 OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc));
2632 OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf);
2633 OCTAVE_LOCAL_BUFFER (double, B, nr);
2634
2635 for (octave_idx_type i = 0; i < nr; i++)
2636 B[i] = N->B[i];
2637
2638 for (octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc)
2639 {
2640 octave_quit ();
2641
2642 for (octave_idx_type j = 0; j < b_nr; j++)
2643 Xx[j] = b.xelem (j, i);
2644
2645 for (octave_idx_type j = nr; j < nbuf; j++)
2646 buf[j] = 0.0;
2647
2648 CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *> (Xx),
2649 buf, nr);
2650 CXSPARSE_ZNAME (_utsolve) (N->U, buf);
2651
2652 for (octave_idx_type j = nr-1; j >= 0; j--)
2653 {
2654 octave_quit ();
2655
2656 CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf);
2657 }
2658
2659 CXSPARSE_ZNAME (_pvec) (S->pinv, buf,
2660 reinterpret_cast<cs_complex_t *> (Xx), nc);
2661
2662 for (octave_idx_type j = 0; j < nc; j++)
2663 {
2664 Complex tmp = Xx[j];
2665
2666 if (tmp != 0.0)
2667 {
2668 if (ii == x_nz)
2669 {
2670 // Resize the sparse matrix
2671 octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
2672 sz = (sz > 10 ? sz : 10) + x_nz;
2673 x.change_capacity (sz);
2674 x_nz = sz;
2675 }
2676
2677 x.xdata (ii) = tmp;
2678 x.xridx (ii++) = j;
2679 }
2680 }
2681
2682 x.xcidx (i+1) = ii;
2683 }
2684
2685 info = 0;
2686
2687 x.maybe_compress ();
2688
2689 return x;
2690
2691#else
2692
2693 octave_unused_parameter (b);
2694
2695 return SparseComplexMatrix ();
2696
2697#endif
2698}
2699
2700template <typename SPARSE_T>
2701sparse_qr<SPARSE_T>::sparse_qr ()
2702 : m_rep (new sparse_qr_rep (SPARSE_T (), 0))
2703{ }
2704
2705template <typename SPARSE_T>
2706sparse_qr<SPARSE_T>::sparse_qr (const SPARSE_T& a, int order)
2707 : m_rep (new sparse_qr_rep (a, order))
2708{ }
2709
2710template <typename SPARSE_T>
2711bool
2712sparse_qr<SPARSE_T>::ok () const
2713{
2714 return m_rep->ok ();
2715}
2716
2717template <typename SPARSE_T>
2718SPARSE_T
2719sparse_qr<SPARSE_T>::V () const
2720{
2721 return m_rep->V ();
2722}
2723
2724template <typename SPARSE_T>
2725ColumnVector
2726sparse_qr<SPARSE_T>::Pinv () const
2727{
2728 return m_rep->P ();
2729}
2730
2731template <typename SPARSE_T>
2732ColumnVector
2733sparse_qr<SPARSE_T>::P () const
2734{
2735 return m_rep->P ();
2736}
2737
2738template <typename SPARSE_T>
2739ColumnVector
2740sparse_qr<SPARSE_T>::E () const
2741{
2742 return m_rep->E();
2743}
2744
2745
2746template <typename SPARSE_T>
2747SparseMatrix
2748sparse_qr<SPARSE_T>::E_MAT () const
2749{
2750 ColumnVector perm = m_rep->E ();
2751 octave_idx_type nrows = perm.rows ();
2752 SparseMatrix ret (nrows, nrows, nrows);
2753 for (octave_idx_type i = 0; i < nrows; i++)
2754 ret(perm(i) - 1, i) = 1.0;
2755 return ret;
2756}
2757
2758
2759template <typename SPARSE_T>
2760SPARSE_T
2761sparse_qr<SPARSE_T>::R (bool econ) const
2762{
2763 return m_rep->R (econ);
2764}
2765
2766template <typename SPARSE_T>
2767typename SPARSE_T::dense_matrix_type
2768sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b,
2769 bool econ) const
2770{
2771 return m_rep->C (b, econ);
2772}
2773
2774template <typename SPARSE_T>
2775typename SPARSE_T::dense_matrix_type
2776sparse_qr<SPARSE_T>::Q (bool econ) const
2777{
2778 return m_rep->Q (econ);
2779}
2780
2781#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
2782//specializations of function min2norm_solve
2783template <>
2784template <>
2785OCTAVE_API Matrix
2786sparse_qr<SparseMatrix>::min2norm_solve<MArray<double>, Matrix>
2787 (const SparseMatrix& a, const MArray<double>& b, octave_idx_type& info, int order)
2788{
2789 info = -1;
2790
2791 octave_idx_type b_nc = b.cols ();
2792 octave_idx_type nc = a.cols ();
2793 Matrix x (nc, b_nc);
2794 cholmod_common cc;
2795
2796 cholmod_l_start (&cc);
2797 cholmod_sparse A = ors2crs (a);
2798 cholmod_dense B = ord2crd (b);
2799 cholmod_dense *X;
2800
2801 X = SuiteSparseQR_min2norm<double> (order, SPQR_DEFAULT_TOL, &A, &B, &cc);
2802 spqr_error_handler (&cc);
2803
2804 double *xdata = x.rwdata ();
2805 for (octave_idx_type i = 0; i < nc * b_nc; i++)
2806 xdata[i] = reinterpret_cast<double *> (X->x)[i];
2807 info = 0;
2808
2809 if (octave_suitesparse_ptr_size_mismatch)
2810 {
2811 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
2812 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
2813 }
2814 cholmod_l_free_dense (&X, &cc);
2815 cholmod_l_finish (&cc);
2816
2817 return x;
2818}
2819
2820template <>
2821template <>
2822OCTAVE_API SparseMatrix
2823sparse_qr<SparseMatrix>::min2norm_solve<SparseMatrix, SparseMatrix>
2824 (const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info, int order)
2825{
2826 info = -1;
2827
2828 SparseMatrix x;
2829 cholmod_common cc;
2830
2831 cholmod_l_start (&cc);
2832 cholmod_sparse A = ors2crs (a);
2833 cholmod_sparse B = ors2crs (b);
2834 cholmod_sparse *X;
2835
2836 X = SuiteSparseQR_min2norm<double> (order, SPQR_DEFAULT_TOL, &A, &B, &cc);
2837 spqr_error_handler (&cc);
2838
2839 x = crs2ors (X, &cc);
2840 info = 0;
2841
2842 if (octave_suitesparse_ptr_size_mismatch)
2843 {
2844 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
2845 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
2846 delete [] reinterpret_cast<SuiteSparse_long *> (B.p);
2847 delete [] reinterpret_cast<SuiteSparse_long *> (B.i);
2848 }
2849 cholmod_l_free_sparse (&X, &cc);
2850 cholmod_l_finish (&cc);
2851
2852 return x;
2853}
2854
2855template <>
2856template <>
2857OCTAVE_API ComplexMatrix
2858sparse_qr<SparseMatrix>::min2norm_solve<MArray<Complex>, ComplexMatrix>
2859 (const SparseMatrix& a, const MArray<Complex>& b, octave_idx_type& info, int order)
2860{
2861 info = -1;
2862
2863 octave_idx_type b_nc = b.cols ();
2864 octave_idx_type nc = a.cols ();
2865 ComplexMatrix x (nc, b_nc);
2866 cholmod_common cc;
2867
2868 cholmod_l_start (&cc);
2869 cholmod_sparse *A = ors2ccs (a, &cc);
2870 cholmod_dense B = ocd2ccd (b);
2871 cholmod_dense *X;
2872
2873 X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, A, &B, &cc);
2874 spqr_error_handler (&cc);
2875
2876 Complex *xdata = x.rwdata ();
2877 for (octave_idx_type i = 0; i < nc * b_nc; i++)
2878 xdata[i] = reinterpret_cast<Complex *> (X->x)[i];
2879 info = 0;
2880
2881 cholmod_l_free_sparse (&A, &cc);
2882 cholmod_l_free_dense (&X, &cc);
2883 cholmod_l_finish (&cc);
2884
2885 return x;
2886}
2887
2888template <>
2889template <>
2890OCTAVE_API SparseComplexMatrix
2891sparse_qr<SparseMatrix>::min2norm_solve<SparseComplexMatrix, SparseComplexMatrix>
2892 (const SparseMatrix& a, const SparseComplexMatrix& b, octave_idx_type& info, int order)
2893{
2894 info = -1;
2895
2896 cholmod_common cc;
2897
2898 cholmod_l_start (&cc);
2899 cholmod_sparse *A = ors2ccs (a, &cc);
2900 cholmod_sparse B = ocs2ccs (b);
2901 cholmod_sparse *X;
2902
2903 X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, A, &B, &cc);
2904 spqr_error_handler (&cc);
2905
2906 SparseComplexMatrix ret = ccs2ocs (X, &cc);
2907 info = 0;
2908
2909 if (octave_suitesparse_ptr_size_mismatch)
2910 {
2911 delete [] reinterpret_cast<SuiteSparse_long *> (B.p);
2912 delete [] reinterpret_cast<SuiteSparse_long *> (B.i);
2913 }
2914 cholmod_l_free_sparse (&A, &cc);
2915 cholmod_l_free_sparse (&X, &cc);
2916 cholmod_l_finish (&cc);
2917
2918 return ret;
2919}
2920
2921template <>
2922template <>
2923OCTAVE_API ComplexMatrix
2924sparse_qr<SparseComplexMatrix>::min2norm_solve<MArray<Complex>, ComplexMatrix>
2925 (const SparseComplexMatrix& a, const MArray<Complex>& b, octave_idx_type& info, int order)
2926{
2927 info = -1;
2928
2929 octave_idx_type b_nc = b.cols ();
2930 octave_idx_type nc = a.cols ();
2931 ComplexMatrix x (nc, b_nc);
2932 cholmod_common cc;
2933
2934 cholmod_l_start (&cc);
2935 cholmod_sparse A = ocs2ccs (a);
2936 cholmod_dense B = ocd2ccd (b);
2937 cholmod_dense *X;
2938
2939 X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, &A, &B, &cc);
2940 spqr_error_handler (&cc);
2941
2942 Complex *xdata = x.rwdata ();
2943 for (octave_idx_type i = 0; i < nc * b_nc; i++)
2944 xdata[i] = reinterpret_cast<Complex *> (X->x)[i];
2945 info = 0;
2946
2947 if (octave_suitesparse_ptr_size_mismatch)
2948 {
2949 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
2950 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
2951 }
2952 cholmod_l_free_dense (&X, &cc);
2953 cholmod_l_finish (&cc);
2954
2955 return x;
2956}
2957
2958// FIXME: 2024/01/07: This template specialization does not appear to be
2959// reachable from current Octave code calling qrsolve from sparse-dmsolve.cc.
2960// It remains because liboctave may be used by other C++ code besides the
2961// Octave interpreter.
2962template <>
2963template <>
2964OCTAVE_API ComplexMatrix
2965sparse_qr<SparseComplexMatrix>::min2norm_solve<MArray<double>, ComplexMatrix>
2966 (const SparseComplexMatrix& a, const MArray<double>& b, octave_idx_type& info, int order)
2967{
2968 info = -1;
2969
2970 octave_idx_type b_nc = b.cols ();
2971 octave_idx_type nc = a.cols ();
2972 ComplexMatrix x (nc, b_nc);
2973 cholmod_common cc;
2974
2975 cholmod_l_start (&cc);
2976 cholmod_sparse A = ocs2ccs (a);
2977 cholmod_dense *B = ord2ccd (b, &cc);
2978 cholmod_dense *X;
2979
2980 X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, &A, B, &cc);
2981 spqr_error_handler (&cc);
2982
2983 Complex *xdata = x.rwdata ();
2984 for (octave_idx_type i = 0; i < nc * b_nc; i++)
2985 xdata[i] = reinterpret_cast<Complex *> (X->x)[i];
2986 info = 0;
2987
2988 if (octave_suitesparse_ptr_size_mismatch)
2989 {
2990 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
2991 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
2992 }
2993 cholmod_l_free_dense (&B, &cc);
2994 cholmod_l_free_dense (&X, &cc);
2995 cholmod_l_finish (&cc);
2996
2997 return x;
2998}
2999
3000template <>
3001template <>
3002OCTAVE_API SparseComplexMatrix
3003sparse_qr<SparseComplexMatrix>::min2norm_solve<SparseComplexMatrix, SparseComplexMatrix>
3004 (const SparseComplexMatrix& a, const SparseComplexMatrix& b, octave_idx_type& info, int order)
3005{
3006 info = -1;
3007
3008 cholmod_common cc;
3009
3010 cholmod_l_start (&cc);
3011 cholmod_sparse A = ocs2ccs (a);
3012 cholmod_sparse B = ocs2ccs (b);
3013 cholmod_sparse *X;
3014
3015 X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, &A, &B, &cc);
3016 spqr_error_handler (&cc);
3017
3018 SparseComplexMatrix ret = ccs2ocs (X, &cc);
3019 info = 0;
3020
3021 if (octave_suitesparse_ptr_size_mismatch)
3022 {
3023 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
3024 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
3025 delete [] reinterpret_cast<SuiteSparse_long *> (B.p);
3026 delete [] reinterpret_cast<SuiteSparse_long *> (B.i);
3027 }
3028 cholmod_l_free_sparse (&X, &cc);
3029 cholmod_l_finish (&cc);
3030
3031 return ret;
3032}
3033
3034// FIXME: 2024/01/07: This template specialization does not appear to be
3035// reachable from current Octave code calling qrsolve from sparse-dmsolve.cc.
3036// It remains because liboctave may be used by other C++ code besides the
3037// Octave interpreter.
3038template <>
3039template <>
3040OCTAVE_API SparseComplexMatrix
3041sparse_qr<SparseComplexMatrix>::min2norm_solve<SparseMatrix, SparseComplexMatrix>
3042 (const SparseComplexMatrix& a, const SparseMatrix& b, octave_idx_type& info, int order)
3043{
3044 info = -1;
3045
3046 cholmod_common cc;
3047
3048 cholmod_l_start (&cc);
3049 cholmod_sparse A = ocs2ccs (a);
3050 cholmod_sparse *B = ors2ccs (b, &cc);
3051 cholmod_sparse *X;
3052
3053 X = SuiteSparseQR_min2norm<Complex> (order, SPQR_DEFAULT_TOL, &A, B, &cc);
3054 spqr_error_handler (&cc);
3055
3056 SparseComplexMatrix ret = ccs2ocs (X, &cc);
3057 info = 0;
3058
3059 if (octave_suitesparse_ptr_size_mismatch)
3060 {
3061 delete [] reinterpret_cast<SuiteSparse_long *> (A.p);
3062 delete [] reinterpret_cast<SuiteSparse_long *> (A.i);
3063 }
3064 cholmod_l_free_sparse (&B, &cc);
3065 cholmod_l_free_sparse (&X, &cc);
3066 cholmod_l_finish (&cc);
3067
3068 return ret;
3069}
3070#endif
3071
3072// FIXME: Why is the "order" of the QR calculation as used in the
3073// CXSparse function sqr 3 for real matrices and 2 for complex? These
3074// values seem to be required but there was no explanation in David
3075// Bateman's original code.
3076
3077template <typename SPARSE_T>
3078class cxsparse_defaults
3079{
3080public:
3081 enum { order = -1 };
3082};
3083
3084template <>
3085class cxsparse_defaults<SparseMatrix>
3086{
3087public:
3088#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
3089 enum { order = SPQR_ORDERING_DEFAULT };
3090#elif defined (HAVE_CXSPARSE)
3091 enum { order = 3 };
3092#endif
3093};
3094
3095template <>
3096class cxsparse_defaults<SparseComplexMatrix>
3097{
3098public:
3099#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
3100 enum { order = SPQR_ORDERING_DEFAULT };
3101#elif defined (HAVE_CXSPARSE)
3102 enum { order = 2 };
3103#endif
3104};
3105
3106template <typename SPARSE_T>
3107template <typename RHS_T, typename RET_T>
3108RET_T
3109sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b,
3110 octave_idx_type& info)
3111{
3112#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
3113
3114 info = -1;
3115
3116 octave_idx_type nr = a.rows ();
3117 octave_idx_type nc = a.cols ();
3118
3119 octave_idx_type b_nc = b.cols ();
3120 octave_idx_type b_nr = b.rows ();
3121
3122 int order = cxsparse_defaults<SPARSE_T>::order;
3123
3124 if (nr < 0 || nc < 0 || b_nc < 0 || b_nr < 0)
3125 (*current_liboctave_error_handler)
3126 ("matrix dimension with negative size");
3127
3128 if ( nr != b_nr)
3129 (*current_liboctave_error_handler)
3130 ("matrix dimension mismatch in solution of minimum norm problem");
3131
3132 info = 0;
3133
3134 return min2norm_solve<RHS_T, RET_T> (a, b, info, order);
3135
3136#elif defined (HAVE_CXSPARSE)
3137
3138 info = -1;
3139
3140 octave_idx_type nr = a.rows ();
3141 octave_idx_type nc = a.cols ();
3142
3143 octave_idx_type b_nc = b.cols ();
3144 octave_idx_type b_nr = b.rows ();
3145
3146 int order = cxsparse_defaults<SPARSE_T>::order;
3147
3148 if (nr < 0 || nc < 0 || nr != b_nr)
3149 (*current_liboctave_error_handler)
3150 ("matrix dimension mismatch in solution of minimum norm problem");
3151
3152 if (nr == 0 || nc == 0 || b_nc == 0)
3153 {
3154 info = 0;
3155
3156 return RET_T (nc, b_nc, 0.0);
3157 }
3158 else if (nr >= nc)
3159 {
3160 sparse_qr<SPARSE_T> q (a, order);
3161
3162 return q.ok () ? q.tall_solve<RHS_T, RET_T> (b, info) : RET_T ();
3163 }
3164 else
3165 {
3166 sparse_qr<SPARSE_T> q (a.hermitian (), order);
3167
3168 return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T ();
3169 }
3170
3171#else
3172
3173 octave_unused_parameter (a);
3174 octave_unused_parameter (b);
3175 octave_unused_parameter (info);
3176
3177 return RET_T ();
3178
3179#endif
3180}
3181
3182//explicit instantiations of static member function solve
3183template
3184OCTAVE_API Matrix
3185sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix>
3186(const SparseMatrix& a, const MArray<double>& b, octave_idx_type& info);
3187
3188template
3189OCTAVE_API SparseMatrix
3190sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix>
3191(const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info);
3192
3193template
3194OCTAVE_API ComplexMatrix
3195sparse_qr<SparseMatrix>::solve<MArray<Complex>, ComplexMatrix>
3196(const SparseMatrix& a, const MArray<Complex>& b, octave_idx_type& info);
3197
3198template
3199OCTAVE_API SparseComplexMatrix
3200sparse_qr<SparseMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix>
3201(const SparseMatrix& a, const SparseComplexMatrix& b,
3202 octave_idx_type& info);
3203
3204template
3205OCTAVE_API ComplexMatrix
3206sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>, ComplexMatrix>
3207(const SparseComplexMatrix& a, const MArray<Complex>& b,
3208 octave_idx_type& info);
3209
3210template
3211OCTAVE_API SparseComplexMatrix
3212sparse_qr<SparseComplexMatrix>::solve<
3213SparseComplexMatrix, SparseComplexMatrix>
3214(const SparseComplexMatrix& a, const SparseComplexMatrix& b,
3215 octave_idx_type& info);
3216
3217template
3218OCTAVE_API ComplexMatrix
3219sparse_qr<SparseComplexMatrix>::solve<MArray<double>, ComplexMatrix>
3220(const SparseComplexMatrix& a, const MArray<double>& b,
3221 octave_idx_type& info);
3222
3223template
3224OCTAVE_API SparseComplexMatrix
3225sparse_qr<SparseComplexMatrix>::solve<SparseMatrix, SparseComplexMatrix>
3226(const SparseComplexMatrix& a, const SparseMatrix& b,
3227 octave_idx_type& info);
3228
3229//explicit instantiations of member function E_MAT
3230template
3231OCTAVE_API SparseMatrix
3232sparse_qr<SparseMatrix>::E_MAT () const;
3233
3234template
3235OCTAVE_API SparseMatrix
3236sparse_qr<SparseComplexMatrix>::E_MAT () const;
3237
3238template <typename SPARSE_T>
3239template <typename RHS_T, typename RET_T>
3240RET_T
3241sparse_qr<SPARSE_T>::tall_solve (const RHS_T& b, octave_idx_type& info) const
3242{
3243 return m_rep->template tall_solve<RHS_T, RET_T> (b, info);
3244}
3245
3246template <typename SPARSE_T>
3247template <typename RHS_T, typename RET_T>
3248RET_T
3249sparse_qr<SPARSE_T>::wide_solve (const RHS_T& b, octave_idx_type& info) const
3250{
3251 return m_rep->template wide_solve<RHS_T, RET_T> (b, info);
3252}
3253
3254// Explicitly instantiate all member functions
3255
3256template OCTAVE_API sparse_qr<SparseMatrix>::sparse_qr ();
3257template OCTAVE_API
3258sparse_qr<SparseMatrix>::sparse_qr (const SparseMatrix& a, int order);
3259template OCTAVE_API bool sparse_qr<SparseMatrix>::ok () const;
3260template OCTAVE_API ColumnVector sparse_qr<SparseMatrix>::E () const;
3261template OCTAVE_API SparseMatrix sparse_qr<SparseMatrix>::V () const;
3262template OCTAVE_API ColumnVector sparse_qr<SparseMatrix>::Pinv () const;
3263template OCTAVE_API ColumnVector sparse_qr<SparseMatrix>::P () const;
3264template OCTAVE_API SparseMatrix
3265sparse_qr<SparseMatrix>::R (bool econ) const;
3266template OCTAVE_API Matrix
3267sparse_qr<SparseMatrix>::C (const Matrix& b, bool econ) const;
3268template OCTAVE_API Matrix sparse_qr<SparseMatrix>::Q (bool econ) const;
3269
3270template OCTAVE_API sparse_qr<SparseComplexMatrix>::sparse_qr ();
3271template OCTAVE_API
3272sparse_qr<SparseComplexMatrix>::sparse_qr
3273(const SparseComplexMatrix& a, int order);
3274template OCTAVE_API bool sparse_qr<SparseComplexMatrix>::ok () const;
3275template OCTAVE_API ColumnVector
3276sparse_qr<SparseComplexMatrix>::E () const;
3277template OCTAVE_API SparseComplexMatrix
3278sparse_qr<SparseComplexMatrix>::V () const;
3279template OCTAVE_API ColumnVector
3280sparse_qr<SparseComplexMatrix>::Pinv () const;
3281template OCTAVE_API ColumnVector
3282sparse_qr<SparseComplexMatrix>::P () const;
3283template OCTAVE_API SparseComplexMatrix
3284sparse_qr<SparseComplexMatrix>::R (bool econ) const;
3285template OCTAVE_API ComplexMatrix
3286sparse_qr<SparseComplexMatrix>::C (const ComplexMatrix& b, bool econ) const;
3287template OCTAVE_API ComplexMatrix
3288sparse_qr<SparseComplexMatrix>::Q (bool econ) const;
3289
3290Matrix
3291qrsolve (const SparseMatrix& a, const MArray<double>& b,
3292 octave_idx_type& info)
3293{
3294 return sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix> (a, b,
3295 info);
3296}
3297
3298SparseMatrix
3299qrsolve (const SparseMatrix& a, const SparseMatrix& b,
3300 octave_idx_type& info)
3301{
3302 return sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix> (a, b,
3303 info);
3304}
3305
3306ComplexMatrix
3307qrsolve (const SparseMatrix& a, const MArray<Complex>& b,
3308 octave_idx_type& info)
3309{
3310 return sparse_qr<SparseMatrix>::solve<MArray<Complex>,
3311 ComplexMatrix> (a, b, info);
3312}
3313
3314SparseComplexMatrix
3315qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b,
3316 octave_idx_type& info)
3317{
3318 return sparse_qr<SparseMatrix>::solve<SparseComplexMatrix,
3319 SparseComplexMatrix> (a, b, info);
3320}
3321
3322ComplexMatrix
3323qrsolve (const SparseComplexMatrix& a, const MArray<double>& b,
3324 octave_idx_type& info)
3325{
3326 return sparse_qr<SparseComplexMatrix>::solve<MArray<double>,
3327 ComplexMatrix> (a, b, info);
3328}
3329
3330SparseComplexMatrix
3331qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b,
3332 octave_idx_type& info)
3333{
3334 return sparse_qr<SparseComplexMatrix>::solve<SparseMatrix,
3335 SparseComplexMatrix>
3336 (a, b, info);
3337}
3338
3339ComplexMatrix
3340qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b,
3341 octave_idx_type& info)
3342{
3343 return sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>,
3344 ComplexMatrix> (a, b, info);
3345}
3346
3347SparseComplexMatrix
3355
3356OCTAVE_END_NAMESPACE(math)
3357OCTAVE_END_NAMESPACE(octave)
C
#define C(a, b)
Definition Faddeeva.cc:256
Array::rows
octave_idx_type rows() const
Definition Array.h:463
Array::cols
octave_idx_type cols() const
Definition Array.h:473
Array::data
const T * data() const
Size of the specified dimension.
Definition Array.h:665
MArray
Template for N-dimensional array classes with like-type math operators.
Definition MArray.h:61
SparseComplexMatrix::hermitian
SparseComplexMatrix hermitian() const
Definition CSparse.cc:604
SparseMatrix::transpose
SparseMatrix transpose() const
Definition dSparse.h:125
Sparse::cols
octave_idx_type cols() const
Definition Sparse.h:349
Sparse::cidx
octave_idx_type * cidx()
Definition Sparse.h:593
Sparse::data
T * data()
Definition Sparse.h:571
Sparse::xdata
T * xdata()
Definition Sparse.h:573
Sparse::ridx
octave_idx_type * ridx()
Definition Sparse.h:580
Sparse::nnz
octave_idx_type nnz() const
Actual number of nonzero terms.
Definition Sparse.h:336
Sparse::rows
octave_idx_type rows() const
Definition Sparse.h:348
Sparse::xelem
T & xelem(octave_idx_type n)
Definition Sparse.h:392
Sparse::xcidx
octave_idx_type * xcidx()
Definition Sparse.h:599
Sparse::xridx
octave_idx_type * xridx()
Definition Sparse.h:586
octave_idx_type
sparse_qr
Definition sparse-qr.h:48
sparse_qr::E
ColumnVector E() const
Definition sparse-qr.cc:2740
sparse_qr::R
SPARSE_T R(bool econ=false) const
Definition sparse-qr.cc:2761
sparse_qr::E_MAT
SparseMatrix E_MAT() const
Definition sparse-qr.cc:2748
sparse_qr::Q
SPARSE_T::dense_matrix_type Q(bool econ=false) const
Definition sparse-qr.cc:2776
sparse_qr::C
SPARSE_T::dense_matrix_type C(const typename SPARSE_T::dense_matrix_type &b, bool econ=false) const
Definition sparse-qr.cc:2768
sparse_qr::ok
bool ok() const
Definition sparse-qr.cc:2712
sparse_qr::solve
static RET_T solve(const SPARSE_T &a, const RHS_T &b, octave_idx_type &info)
Definition sparse-qr.cc:3109
sparse_qr::V
SPARSE_T V() const
Definition sparse-qr.cc:2719
sparse_qr::Pinv
ColumnVector Pinv() const
Definition sparse-qr.cc:2726
sparse_qr::P
ColumnVector P() const
Definition sparse-qr.cc:2733
sparse_qr::sparse_qr
sparse_qr()
Definition sparse-qr.cc:2701
OCTAVE_BEGIN_NAMESPACE
OCTAVE_BEGIN_NAMESPACE(octave) static octave_value daspk_fcn
current_liboctave_error_handler
OCTAVE_NORETURN liboctave_error_handler current_liboctave_error_handler
Definition lo-error.c:41
B
F77_RET_T const F77_INT F77_CMPLX const F77_INT F77_CMPLX * B
Definition lo-lapack-proto.h:859
Q
F77_RET_T const F77_INT const F77_INT const F77_INT F77_DBLE const F77_INT F77_DBLE const F77_INT F77_DBLE * Q
Definition lo-lapack-proto.h:995
A
F77_RET_T const F77_INT F77_CMPLX * A
Definition lo-lapack-proto.h:858
N
F77_RET_T const F77_INT & N
Definition lo-lapack-proto.h:857
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
OCTAVE_LOCAL_BUFFER
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition oct-locbuf.h:44
to_suitesparse_intptr
suitesparse_integer * to_suitesparse_intptr(octave_idx_type *i)
Definition oct-sparse.cc:51
CXSPARSE_DNAME
#define CXSPARSE_DNAME(name)
Definition oct-sparse.h:159
suitesparse_integer
int suitesparse_integer
Definition oct-sparse.h:184
CXSPARSE_ZNAME
#define CXSPARSE_ZNAME(name)
Definition oct-sparse.h:160
from_suitesparse_long
octave_idx_type from_suitesparse_long(SuiteSparse_long x)
Definition oct-sparse.h:200
from_size_t
octave_idx_type from_size_t(std::size_t x)
Definition oct-sparse.h:211
free
void free(void *)
qrsolve
Matrix qrsolve(const SparseMatrix &a, const MArray< double > &b, octave_idx_type &info)
Definition sparse-qr.cc:3291