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