Mercurial > octave
changeset 21176:791dcb32b657
revamp sparse QR factorizatino classes
* sparse-qr-inst.cc, sparse-qr.cc, sparse-qr.h: New files, adapted
from SparseCmplxQR.cc, SparseCmplxQR.h, SparseQR.cc, and SparseQR.h.
(sparse_qr): New template class.
* SparseCmplxQR.cc, SparseCmplxQR.h, SparseQR.cc, SparseQR.h: Delete.
* liboctave/numeric/module.mk: Update.
* dmperm.cc, qr.cc, CSparse.cc, dSparse.cc, sparse-dmsolve.cc:
Use new sparse_qr class.
* oct-sparse.h (CXSPARSE_DNAME, CXSPARSE_ZNAME):
Move macros here from SparseQR.h and SparseCmplxQR.h.
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 02 Feb 2016 12:34:06 -0500 |
| parents | f45f4f888db5 |
| children | a10f60e13243 |
| files | libinterp/dldfcn/dmperm.cc libinterp/dldfcn/qr.cc liboctave/array/CSparse.cc liboctave/array/dSparse.cc liboctave/numeric/SparseCmplxQR.cc liboctave/numeric/SparseCmplxQR.h liboctave/numeric/SparseQR.cc liboctave/numeric/SparseQR.h liboctave/numeric/module.mk liboctave/numeric/sparse-dmsolve.cc liboctave/numeric/sparse-qr-inst.cc liboctave/numeric/sparse-qr.cc liboctave/numeric/sparse-qr.h liboctave/util/oct-sparse.h |
| diffstat | 14 files changed, 2150 insertions(+), 1888 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/dldfcn/dmperm.cc Tue Feb 02 12:33:38 2016 -0500 +++ b/libinterp/dldfcn/dmperm.cc Tue Feb 02 12:34:06 2016 -0500 @@ -34,8 +34,7 @@ #include "oct-sparse.h" #include "ov-re-sparse.h" #include "ov-cx-sparse.h" -#include "SparseQR.h" -#include "SparseCmplxQR.h" +#include "sparse-qr.h" #if defined (ENABLE_64) # define CXSPARSE_NAME(name) cs_dl ## name
--- a/libinterp/dldfcn/qr.cc Tue Feb 02 12:33:38 2016 -0500 +++ b/libinterp/dldfcn/qr.cc Tue Feb 02 12:34:06 2016 -0500 @@ -34,8 +34,7 @@ #include "fCmplxQRP.h" #include "floatQR.h" #include "floatQRP.h" -#include "SparseQR.h" -#include "SparseCmplxQR.h" +#include "sparse-qr.h" #include "defun-dld.h" @@ -239,7 +238,7 @@ if (is_cmplx) { - SparseComplexQR q (arg.sparse_complex_matrix_value ()); + sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ()); if (have_b > 0) { @@ -255,7 +254,7 @@ } else { - SparseQR q (arg.sparse_matrix_value ()); + sparse_qr<SparseMatrix> q (arg.sparse_matrix_value ()); if (have_b > 0) {
--- a/liboctave/array/CSparse.cc Tue Feb 02 12:33:38 2016 -0500 +++ b/liboctave/array/CSparse.cc Tue Feb 02 12:34:06 2016 -0500 @@ -55,7 +55,7 @@ #include "oct-sparse.h" #include "sparse-util.h" #include "sparse-chol.h" -#include "SparseCmplxQR.h" +#include "sparse-qr.h" #include "Sparse-op-defs.h"
--- a/liboctave/array/dSparse.cc Tue Feb 02 12:33:38 2016 -0500 +++ b/liboctave/array/dSparse.cc Tue Feb 02 12:34:06 2016 -0500 @@ -50,7 +50,7 @@ #include "oct-sparse.h" #include "sparse-util.h" #include "sparse-chol.h" -#include "SparseQR.h" +#include "sparse-qr.h" #include "Sparse-op-defs.h"
--- a/liboctave/numeric/SparseCmplxQR.cc Tue Feb 02 12:33:38 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,736 +0,0 @@ -/* - -Copyright (C) 2005-2015 David Bateman - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif -#include <vector> - -#include "lo-error.h" -#include "SparseCmplxQR.h" -#include "oct-locbuf.h" - -SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep - (OCTAVE_UNUSED const SparseComplexMatrix& a, OCTAVE_UNUSED int order) - : count (1), nrows (0) -#ifdef HAVE_CXSPARSE - , S (0), N (0) -#endif -{ -#ifdef HAVE_CXSPARSE - CXSPARSE_ZNAME () A; - A.nzmax = a.nnz (); - A.m = a.rows (); - A.n = a.cols (); - nrows = A.m; - // Cast away const on A, with full knowledge that CSparse won't touch it - // Prevents the methods below making a copy of the data. - A.p = const_cast<octave_idx_type *>(a.cidx ()); - A.i = const_cast<octave_idx_type *>(a.ridx ()); - A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *> - (a.data ())); - A.nz = -1; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - S = CXSPARSE_ZNAME (_sqr) (order, &A, 1); - N = CXSPARSE_ZNAME (_qr) (&A, S); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - if (! N) - (*current_liboctave_error_handler) - ("SparseComplexQR: sparse matrix QR factorization filled"); - - count = 1; - -#else - (*current_liboctave_error_handler) - ("SparseComplexQR: support for CXSparse was unavailable or disabled when liboctave was built"); -#endif -} - -SparseComplexQR::SparseComplexQR_rep::~SparseComplexQR_rep (void) -{ -#ifdef HAVE_CXSPARSE - CXSPARSE_ZNAME (_sfree) (S); - CXSPARSE_ZNAME (_nfree) (N); -#endif -} - -SparseComplexMatrix -SparseComplexQR::SparseComplexQR_rep::V (void) const -{ -#ifdef HAVE_CXSPARSE - // Drop zeros from V and sort - // FIXME: Is the double transpose to sort necessary? - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_dropzeros) (N->L); - CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1); - CXSPARSE_ZNAME (_spfree) (N->L); - N->L = CXSPARSE_ZNAME (_transpose) (D, 1); - CXSPARSE_ZNAME (_spfree) (D); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - octave_idx_type nc = N->L->n; - octave_idx_type nz = N->L->nzmax; - SparseComplexMatrix ret (N->L->m, nc, nz); - for (octave_idx_type j = 0; j < nc+1; j++) - ret.xcidx (j) = N->L->p[j]; - for (octave_idx_type j = 0; j < nz; j++) - { - ret.xridx (j) = N->L->i[j]; - ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j]; - } - return ret; -#else - return SparseComplexMatrix (); -#endif -} - -ColumnVector -SparseComplexQR::SparseComplexQR_rep::Pinv (void) const -{ -#ifdef HAVE_CXSPARSE - ColumnVector ret(N->L->m); - for (octave_idx_type i = 0; i < N->L->m; i++) - ret.xelem (i) = S->pinv[i]; - return ret; -#else - return ColumnVector (); -#endif -} - -ColumnVector -SparseComplexQR::SparseComplexQR_rep::P (void) const -{ -#ifdef HAVE_CXSPARSE - ColumnVector ret(N->L->m); - for (octave_idx_type i = 0; i < N->L->m; i++) - ret.xelem (S->pinv[i]) = i; - return ret; -#else - return ColumnVector (); -#endif -} - -SparseComplexMatrix -SparseComplexQR::SparseComplexQR_rep::R (const bool econ) const -{ -#ifdef HAVE_CXSPARSE - // Drop zeros from R and sort - // FIXME: Is the double transpose to sort necessary? - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_dropzeros) (N->U); - CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1); - CXSPARSE_ZNAME (_spfree) (N->U); - N->U = CXSPARSE_ZNAME (_transpose) (D, 1); - CXSPARSE_ZNAME (_spfree) (D); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - octave_idx_type nc = N->U->n; - octave_idx_type nz = N->U->nzmax; - SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); - for (octave_idx_type j = 0; j < nc+1; j++) - ret.xcidx (j) = N->U->p[j]; - for (octave_idx_type j = 0; j < nz; j++) - { - ret.xridx (j) = N->U->i[j]; - ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j]; - } - return ret; -#else - return SparseComplexMatrix (); -#endif -} - -ComplexMatrix -SparseComplexQR::SparseComplexQR_rep::C (const ComplexMatrix &b) const -{ -#ifdef HAVE_CXSPARSE - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - octave_idx_type nc = N->L->n; - octave_idx_type nr = nrows; - const cs_complex_t *bvec = - reinterpret_cast<const cs_complex_t *>(b.fortran_vec ()); - ComplexMatrix ret(b_nr, b_nc); - Complex *vec = ret.fortran_vec (); - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) ("matrix dimension mismatch"); - - if (nr == 0 || nc == 0 || b_nc == 0) - ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); - else - { - OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); - for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) - { - octave_quit (); - volatile octave_idx_type nm = (nr < nc ? nr : nc); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_ipvec) - (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type i = 0; i < nm; i++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) - (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - for (octave_idx_type i = 0; i < b_nr; i++) - vec[i+idx] = buf[i]; - } - } - return ret; -#else - return ComplexMatrix (); -#endif -} - -ComplexMatrix -SparseComplexQR::SparseComplexQR_rep::Q (void) const -{ -#ifdef HAVE_CXSPARSE - octave_idx_type nc = N->L->n; - octave_idx_type nr = nrows; - ComplexMatrix ret(nr, nr); - Complex *vec = ret.fortran_vec (); - if (nr < 0 || nc < 0) - (*current_liboctave_error_handler) ("matrix dimension mismatch"); - - if (nr == 0 || nc == 0) - ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0)); - else - { - OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr); - for (octave_idx_type i = 0; i < nr; i++) - bvec[i] = cs_complex_t (0.0, 0.0); - OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); - for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) - { - octave_quit (); - bvec[j] = cs_complex_t (1.0, 0.0); - volatile octave_idx_type nm = (nr < nc ? nr : nc); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_ipvec) - (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type i = 0; i < nm; i++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) - (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - for (octave_idx_type i = 0; i < nr; i++) - vec[i+idx] = buf[i]; - bvec[j] = cs_complex_t (0.0, 0.0); - } - } - return ret.hermitian (); -#else - return ComplexMatrix (); -#endif -} - -ComplexMatrix -qrsolve (const SparseComplexMatrix& a, const MArray<double>& b, - octave_idx_type& info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nc = b.cols (); - octave_idx_type b_nr = b.rows (); - ComplexMatrix x; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); - else if (nr >= nc) - { - SparseComplexQR q (a, 2); - if (! q.ok ()) - return ComplexMatrix (); - x.resize (nc, b_nc); - cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2); - OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_ipvec) - (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_ZNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - info = 0; - } - else - { - SparseComplexMatrix at = a.hermitian (); - SparseComplexQR q (at, 2); - if (! q.ok ()) - return ComplexMatrix (); - x.resize (nc, b_nc); - cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); - OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); - OCTAVE_LOCAL_BUFFER (double, B, nr); - for (octave_idx_type i = 0; i < nr; i++) - B[i] = q.N ()->B[i]; - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) - (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); - CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - info = 0; - } - - return x; -#else - return ComplexMatrix (); -#endif -} - -SparseComplexMatrix -qrsolve (const SparseComplexMatrix&a, const SparseMatrix &b, - octave_idx_type &info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nc = b.cols (); - octave_idx_type b_nr = b.rows (); - SparseComplexMatrix x; - volatile octave_idx_type ii, x_nz; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = SparseComplexMatrix (nc, b_nc); - else if (nr >= nc) - { - SparseComplexQR q (a, 2); - if (! q.ok ()) - return SparseComplexMatrix (); - x = SparseComplexMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_ipvec) - (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_ZNAME (_ipvec) - (q.S ()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - Complex tmp = Xx[j]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - else - { - SparseComplexMatrix at = a.hermitian (); - SparseComplexQR q (at, 2); - if (! q.ok ()) - return SparseComplexMatrix (); - x = SparseComplexMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); - - OCTAVE_LOCAL_BUFFER (double, B, nr); - for (octave_idx_type i = 0; i < nr; i++) - B[i] = q.N ()->B[i]; - - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) - (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); - CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) - (q.S ()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - Complex tmp = Xx[j]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - - x.maybe_compress (); - return x; -#else - return SparseComplexMatrix (); -#endif -} - -ComplexMatrix -qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b, - octave_idx_type& info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nc = b.cols (); - octave_idx_type b_nr = b.rows (); - const cs_complex_t *bvec = - reinterpret_cast<const cs_complex_t *>(b.fortran_vec ()); - ComplexMatrix x; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); - else if (nr >= nc) - { - SparseComplexQR q (a, 2); - if (! q.ok ()) - return ComplexMatrix (); - x.resize (nc, b_nc); - cs_complex_t *vec = reinterpret_cast<cs_complex_t *> - (x.fortran_vec ()); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; - i++, idx+=nc, bidx+=b_nr) - { - octave_quit (); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_ZNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - info = 0; - } - else - { - SparseComplexMatrix at = a.hermitian (); - SparseComplexQR q (at, 2); - if (! q.ok ()) - return ComplexMatrix (); - x.resize (nc, b_nc); - cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); - OCTAVE_LOCAL_BUFFER (double, B, nr); - for (octave_idx_type i = 0; i < nr; i++) - B[i] = q.N ()->B[i]; - for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; - i++, idx+=nc, bidx+=b_nr) - { - octave_quit (); - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr); - CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - info = 0; - } - - return x; -#else - return ComplexMatrix (); -#endif -} - -SparseComplexMatrix -qrsolve (const SparseComplexMatrix&a, const SparseComplexMatrix &b, - octave_idx_type &info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nc = b.cols (); - octave_idx_type b_nr = b.rows (); - SparseComplexMatrix x; - volatile octave_idx_type ii, x_nz; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = SparseComplexMatrix (nc, b_nc); - else if (nr >= nc) - { - SparseComplexQR q (a, 2); - if (! q.ok ()) - return SparseComplexMatrix (); - x = SparseComplexMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_ipvec) - (q.S ()->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_ZNAME (_ipvec) - (q.S ()->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - Complex tmp = Xx[j]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - else - { - SparseComplexMatrix at = a.hermitian (); - SparseComplexQR q (at, 2); - if (! q.ok ()) - return SparseComplexMatrix (); - x = SparseComplexMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); - OCTAVE_LOCAL_BUFFER (double, B, nr); - for (octave_idx_type i = 0; i < nr; i++) - B[i] = q.N ()->B[i]; - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = cs_complex_t (0.0, 0.0); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) - (q.S ()->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); - CXSPARSE_ZNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_happly) (q.N ()->L, j, B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_ZNAME (_pvec) - (q.S ()->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - Complex tmp = Xx[j]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - - x.maybe_compress (); - return x; -#else - return SparseComplexMatrix (); -#endif -}
--- a/liboctave/numeric/SparseCmplxQR.h Tue Feb 02 12:33:38 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,173 +0,0 @@ -/* - -Copyright (C) 2005-2015 David Bateman - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#if ! defined (octave_SparseCmplxQR_h) -#define octave_SparseCmplxQR_h 1 - -#include <iosfwd> - -#include "dMatrix.h" -#include "CMatrix.h" -#include "dSparse.h" -#include "CSparse.h" -#include "oct-sparse.h" - -#if defined (ENABLE_64) -#define CXSPARSE_ZNAME(name) cs_cl ## name -#else -#define CXSPARSE_ZNAME(name) cs_ci ## name -#endif - -class -OCTAVE_API -SparseComplexQR -{ -protected: - class SparseComplexQR_rep - { - public: - SparseComplexQR_rep (const SparseComplexMatrix& a, int order); - - ~SparseComplexQR_rep (void); -#ifdef HAVE_CXSPARSE - bool ok (void) const { return (N && S); } -#else - bool ok (void) const { return false; } -#endif - SparseComplexMatrix V (void) const; - - ColumnVector Pinv (void) const; - - ColumnVector P (void) const; - - SparseComplexMatrix R (const bool econ) const; - - ComplexMatrix C (const ComplexMatrix &b) const; - - ComplexMatrix Q (void) const; - - octave_refcount<int> count; - - octave_idx_type nrows; -#ifdef HAVE_CXSPARSE - CXSPARSE_ZNAME (s) *S; - - CXSPARSE_ZNAME (n) *N; -#endif - private: - - // No copying! - - SparseComplexQR_rep (const SparseComplexQR_rep&); - - SparseComplexQR_rep operator = (const SparseComplexQR_rep&); - - }; -private: - SparseComplexQR_rep *rep; - -public: - SparseComplexQR (void) : - rep (new SparseComplexQR_rep (SparseComplexMatrix (), 0)) { } - - SparseComplexQR (const SparseComplexMatrix& a, int order = 0) : - rep (new SparseComplexQR_rep (a, order)) { } - - SparseComplexQR (const SparseComplexQR& a) : rep (a.rep) { rep->count++; } - - ~SparseComplexQR (void) - { - if (--rep->count == 0) - delete rep; - } - - SparseComplexQR& operator = (const SparseComplexQR& a) - { - if (this != &a) - { - if (--rep->count == 0) - delete rep; - - rep = a.rep; - rep->count++; - } - return *this; - } - - bool ok (void) const { return rep->ok (); } - - SparseComplexMatrix V (void) const { return rep->V (); } - - ColumnVector Pinv (void) const { return rep->P (); } - - ColumnVector P (void) const { return rep->P (); } - - SparseComplexMatrix R (const bool econ = false) const - { return rep->R(econ); } - - ComplexMatrix C (const ComplexMatrix &b) const { return rep->C(b); } - - ComplexMatrix Q (void) const { return rep->Q (); } - - friend ComplexMatrix qrsolve (const SparseComplexMatrix &a, - const MArray<double> &b, - octave_idx_type &info); - - friend SparseComplexMatrix qrsolve (const SparseComplexMatrix &a, - const SparseMatrix &b, - octave_idx_type &info); - - friend ComplexMatrix qrsolve (const SparseComplexMatrix &a, - const MArray<Complex> &b, - octave_idx_type &info); - - friend SparseComplexMatrix qrsolve (const SparseComplexMatrix &a, - const SparseComplexMatrix &b, - octave_idx_type &info); - -protected: -#ifdef HAVE_CXSPARSE - CXSPARSE_ZNAME (s) * S (void) { return rep->S; } - - CXSPARSE_ZNAME (n) * N (void) { return rep->N; } -#endif -}; - - -// Publish externally used friend functions. - -extern ComplexMatrix qrsolve (const SparseComplexMatrix &a, - const MArray<double> &b, - octave_idx_type &info); - -extern SparseComplexMatrix qrsolve (const SparseComplexMatrix &a, - const SparseMatrix &b, - octave_idx_type &info); - -extern ComplexMatrix qrsolve (const SparseComplexMatrix &a, - const MArray<Complex> &b, - octave_idx_type &info); - -extern SparseComplexMatrix qrsolve (const SparseComplexMatrix &a, - const SparseComplexMatrix &b, - octave_idx_type &info); -#endif
--- a/liboctave/numeric/SparseQR.cc Tue Feb 02 12:33:38 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,795 +0,0 @@ -/* - -Copyright (C) 2005-2015 David Bateman - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif -#include <vector> - -#include "lo-error.h" -#include "SparseQR.h" -#include "oct-locbuf.h" - -SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order) - : count (1), nrows (0) -#ifdef HAVE_CXSPARSE - , S (0), N (0) -#endif -{ -#ifdef HAVE_CXSPARSE - CXSPARSE_DNAME () A; - A.nzmax = a.nnz (); - A.m = a.rows (); - A.n = a.cols (); - nrows = A.m; - // Cast away const on A, with full knowledge that CSparse won't touch it - // Prevents the methods below making a copy of the data. - A.p = const_cast<octave_idx_type *>(a.cidx ()); - A.i = const_cast<octave_idx_type *>(a.ridx ()); - A.x = const_cast<double *>(a.data ()); - A.nz = -1; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - S = CXSPARSE_DNAME (_sqr) (order, &A, 1); - - N = CXSPARSE_DNAME (_qr) (&A, S); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - if (! N) - (*current_liboctave_error_handler) - ("SparseQR: sparse matrix QR factorization filled"); - - count = 1; - -#else - (*current_liboctave_error_handler) - ("SparseQR: support for CXSparse was unavailable or disabled when liboctave was built"); -#endif -} - -SparseQR::SparseQR_rep::~SparseQR_rep (void) -{ -#ifdef HAVE_CXSPARSE - CXSPARSE_DNAME (_sfree) (S); - CXSPARSE_DNAME (_nfree) (N); -#endif -} - -SparseMatrix -SparseQR::SparseQR_rep::V (void) const -{ -#ifdef HAVE_CXSPARSE - // Drop zeros from V and sort - // FIXME: Is the double transpose to sort necessary? - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_dropzeros) (N->L); - CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1); - CXSPARSE_DNAME (_spfree) (N->L); - N->L = CXSPARSE_DNAME (_transpose) (D, 1); - CXSPARSE_DNAME (_spfree) (D); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - octave_idx_type nc = N->L->n; - octave_idx_type nz = N->L->nzmax; - SparseMatrix ret (N->L->m, nc, nz); - for (octave_idx_type j = 0; j < nc+1; j++) - ret.xcidx (j) = N->L->p[j]; - for (octave_idx_type j = 0; j < nz; j++) - { - ret.xridx (j) = N->L->i[j]; - ret.xdata (j) = N->L->x[j]; - } - return ret; -#else - return SparseMatrix (); -#endif -} - -ColumnVector -SparseQR::SparseQR_rep::Pinv (void) const -{ -#ifdef HAVE_CXSPARSE - ColumnVector ret(N->L->m); - for (octave_idx_type i = 0; i < N->L->m; i++) - ret.xelem (i) = S->pinv[i]; - return ret; -#else - return ColumnVector (); -#endif -} - -ColumnVector -SparseQR::SparseQR_rep::P (void) const -{ -#ifdef HAVE_CXSPARSE - ColumnVector ret(N->L->m); - for (octave_idx_type i = 0; i < N->L->m; i++) - ret.xelem (S->pinv[i]) = i; - return ret; -#else - return ColumnVector (); -#endif -} - -SparseMatrix -SparseQR::SparseQR_rep::R (const bool econ) const -{ -#ifdef HAVE_CXSPARSE - // Drop zeros from R and sort - // FIXME: Is the double transpose to sort necessary? - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_dropzeros) (N->U); - CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1); - CXSPARSE_DNAME (_spfree) (N->U); - N->U = CXSPARSE_DNAME (_transpose) (D, 1); - CXSPARSE_DNAME (_spfree) (D); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - octave_idx_type nc = N->U->n; - octave_idx_type nz = N->U->nzmax; - - SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); - - for (octave_idx_type j = 0; j < nc+1; j++) - ret.xcidx (j) = N->U->p[j]; - for (octave_idx_type j = 0; j < nz; j++) - { - ret.xridx (j) = N->U->i[j]; - ret.xdata (j) = N->U->x[j]; - } - return ret; -#else - return SparseMatrix (); -#endif -} - -Matrix -SparseQR::SparseQR_rep::C (const Matrix &b) const -{ -#ifdef HAVE_CXSPARSE - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - octave_idx_type nc = N->L->n; - octave_idx_type nr = nrows; - const double *bvec = b.fortran_vec (); - Matrix ret (b_nr, b_nc); - double *vec = ret.fortran_vec (); - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) ("matrix dimension mismatch"); - - if (nr == 0 || nc == 0 || b_nc == 0) - ret = Matrix (nc, b_nc, 0.0); - else - { - OCTAVE_LOCAL_BUFFER (double, buf, S->m2); - for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) - { - octave_quit (); - for (octave_idx_type i = nr; i < S->m2; i++) - buf[i] = 0.; - volatile octave_idx_type nm = (nr < nc ? nr : nc); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (volatile octave_idx_type i = 0; i < nm; i++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - for (octave_idx_type i = 0; i < b_nr; i++) - vec[i+idx] = buf[i]; - } - } - return ret; -#else - return Matrix (); -#endif -} - -Matrix -SparseQR::SparseQR_rep::Q (void) const -{ -#ifdef HAVE_CXSPARSE - octave_idx_type nc = N->L->n; - octave_idx_type nr = nrows; - Matrix ret (nr, nr); - double *vec = ret.fortran_vec (); - if (nr < 0 || nc < 0) - (*current_liboctave_error_handler) ("matrix dimension mismatch"); - - if (nr == 0 || nc == 0) - ret = Matrix (nc, nr, 0.0); - else - { - OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1); - for (octave_idx_type i = 0; i < nr; i++) - bvec[i] = 0.; - OCTAVE_LOCAL_BUFFER (double, buf, S->m2); - for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) - { - octave_quit (); - bvec[j] = 1.0; - for (octave_idx_type i = nr; i < S->m2; i++) - buf[i] = 0.; - volatile octave_idx_type nm = (nr < nc ? nr : nc); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (volatile octave_idx_type i = 0; i < nm; i++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - for (octave_idx_type i = 0; i < nr; i++) - vec[i+idx] = buf[i]; - bvec[j] = 0.0; - } - } - return ret.transpose (); -#else - return Matrix (); -#endif -} - -Matrix -qrsolve (const SparseMatrix& a, const MArray<double>& b, octave_idx_type& info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nc = b.cols (); - octave_idx_type b_nr = b.rows (); - const double *bvec = b.fortran_vec (); - Matrix x; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = Matrix (nc, b_nc, 0.0); - else if (nr >= nc) - { - SparseQR q (a, 3); - if (! q.ok ()) - return Matrix (); - x.resize (nc, b_nc); - double *vec = x.fortran_vec (); - OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; - i++, idx+=nc, bidx+=b_nr) - { - octave_quit (); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, bvec + bidx, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, vec + idx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - info = 0; - } - else - { - SparseMatrix at = a.hermitian (); - SparseQR q (at, 3); - if (! q.ok ()) - return Matrix (); - x.resize (nc, b_nc); - double *vec = x.fortran_vec (); - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (double, buf, nbuf); - for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; - i++, idx+=nc, bidx+=b_nr) - { - octave_quit (); - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->q, bvec + bidx, buf, nr); - CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, vec + idx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - info = 0; - } - - return x; -#else - return Matrix (); -#endif -} - -SparseMatrix -qrsolve (const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - SparseMatrix x; - volatile octave_idx_type ii, x_nz; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = SparseMatrix (nc, b_nc); - else if (nr >= nc) - { - SparseQR q (a, 3); - if (! q.ok ()) - return SparseMatrix (); - x = SparseMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - double tmp = Xx[j]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - else - { - SparseMatrix at = a.hermitian (); - SparseQR q (at, 3); - if (! q.ok ()) - return SparseMatrix (); - x = SparseMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, buf, nbuf); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - Xx[j] = b.xelem (j,i); - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr); - CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - double tmp = Xx[j]; - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - - x.maybe_compress (); - return x; -#else - return SparseMatrix (); -#endif -} - -ComplexMatrix -qrsolve (const SparseMatrix& a, const MArray<Complex>& b, octave_idx_type& info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nc = b.cols (); - octave_idx_type b_nr = b.rows (); - ComplexMatrix x; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); - else if (nr >= nc) - { - SparseQR q (a, 3); - if (! q.ok ()) - return ComplexMatrix (); - x.resize (nc, b_nc); - Complex *vec = x.fortran_vec (); - OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - { - Complex c = b.xelem (j,i); - Xx[j] = std::real (c); - Xz[j] = std::imag (c); - } - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc); - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = 0.; - CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (octave_idx_type j = 0; j < nc; j++) - vec[j+idx] = Complex (Xx[j], Xz[j]); - } - info = 0; - } - else - { - SparseMatrix at = a.hermitian (); - SparseQR q (at, 3); - if (! q.ok ()) - return ComplexMatrix (); - x.resize (nc, b_nc); - Complex *vec = x.fortran_vec (); - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, buf, nbuf); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - { - Complex c = b.xelem (j,i); - Xx[j] = std::real (c); - Xz[j] = std::imag (c); - } - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr); - CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr); - CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (octave_idx_type j = 0; j < nc; j++) - vec[j+idx] = Complex (Xx[j], Xz[j]); - } - info = 0; - } - - return x; -#else - return ComplexMatrix (); -#endif -} - -SparseComplexMatrix -qrsolve (const SparseMatrix&a, const SparseComplexMatrix &b, - octave_idx_type &info) -{ - info = -1; -#ifdef HAVE_CXSPARSE - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - octave_idx_type b_nr = b.rows (); - octave_idx_type b_nc = b.cols (); - SparseComplexMatrix x; - volatile octave_idx_type ii, x_nz; - - if (nr < 0 || nc < 0 || nr != b_nr) - (*current_liboctave_error_handler) - ("matrix dimension mismatch in solution of minimum norm problem"); - - if (nr == 0 || nc == 0 || b_nc == 0) - x = SparseComplexMatrix (nc, b_nc); - else if (nr >= nc) - { - SparseQR q (a, 3); - if (! q.ok ()) - return SparseComplexMatrix (); - x = SparseComplexMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, buf, q.S ()->m2); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - { - Complex c = b.xelem (j,i); - Xx[j] = std::real (c); - Xz[j] = std::imag (c); - } - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xx, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (octave_idx_type j = nr; j < q.S ()->m2; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_ipvec) (q.S ()->pinv, Xz, buf, nr); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = 0; j < nc; j++) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_usolve) (q.N ()->U, buf); - CXSPARSE_DNAME (_ipvec) (q.S ()->q, buf, Xz, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - Complex tmp = Complex (Xx[j], Xz[j]); - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - else - { - SparseMatrix at = a.hermitian (); - SparseQR q (at, 3); - if (! q.ok ()) - return SparseComplexMatrix (); - x = SparseComplexMatrix (nc, b_nc, b.nnz ()); - x.xcidx (0) = 0; - x_nz = b.nnz (); - ii = 0; - volatile octave_idx_type nbuf = (nc > q.S ()->m2 ? nc : q.S ()->m2); - OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); - OCTAVE_LOCAL_BUFFER (double, buf, nbuf); - for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) - { - octave_quit (); - for (octave_idx_type j = 0; j < b_nr; j++) - { - Complex c = b.xelem (j,i); - Xx[j] = std::real (c); - Xz[j] = std::imag (c); - } - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->q, Xx, buf, nr); - CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xx, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (octave_idx_type j = nr; j < nbuf; j++) - buf[j] = 0.; - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->q, Xz, buf, nr); - CXSPARSE_DNAME (_utsolve) (q.N ()->U, buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - for (volatile octave_idx_type j = nr-1; j >= 0; j--) - { - octave_quit (); - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_happly) (q.N ()->L, j, q.N ()->B[j], buf); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - } - BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - CXSPARSE_DNAME (_pvec) (q.S ()->pinv, buf, Xz, nc); - END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; - - for (octave_idx_type j = 0; j < nc; j++) - { - Complex tmp = Complex (Xx[j], Xz[j]); - if (tmp != 0.0) - { - if (ii == x_nz) - { - // Resize the sparse matrix - octave_idx_type sz = x_nz * (b_nc - i) / b_nc; - sz = (sz > 10 ? sz : 10) + x_nz; - x.change_capacity (sz); - x_nz = sz; - } - x.xdata (ii) = tmp; - x.xridx (ii++) = j; - } - } - x.xcidx (i+1) = ii; - } - info = 0; - } - - x.maybe_compress (); - return x; -#else - return SparseComplexMatrix (); -#endif -}
--- a/liboctave/numeric/SparseQR.h Tue Feb 02 12:33:38 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,169 +0,0 @@ -/* - -Copyright (C) 2005-2015 David Bateman - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#if ! defined (octave_SparseQR_h) -#define octave_SparseQR_h 1 - -#include <iosfwd> - -#include "dMatrix.h" -#include "CMatrix.h" -#include "dSparse.h" -#include "CSparse.h" -#include "oct-sparse.h" - -#if defined (ENABLE_64) -#define CXSPARSE_DNAME(name) cs_dl ## name -#else -#define CXSPARSE_DNAME(name) cs_di ## name -#endif - -class -OCTAVE_API -SparseQR -{ -protected: - class SparseQR_rep - { - public: - SparseQR_rep (const SparseMatrix& a, int order); - - ~SparseQR_rep (void); -#ifdef HAVE_CXSPARSE - bool ok (void) const { return (N && S); } -#else - bool ok (void) const { return false; } -#endif - SparseMatrix V (void) const; - - ColumnVector Pinv (void) const; - - ColumnVector P (void) const; - - SparseMatrix R (const bool econ) const; - - Matrix C (const Matrix &b) const; - - Matrix Q (void) const; - - octave_refcount<int> count; - - octave_idx_type nrows; -#ifdef HAVE_CXSPARSE - CXSPARSE_DNAME (s) *S; - - CXSPARSE_DNAME (n) *N; -#endif - - private: - - // No copying! - - SparseQR_rep (const SparseQR_rep&); - - SparseQR_rep& operator = (const SparseQR_rep&); - }; - -private: - - SparseQR_rep *rep; - -public: - - SparseQR (void) : rep (new SparseQR_rep (SparseMatrix (), 0)) { } - - SparseQR (const SparseMatrix& a, int order = 0) : - rep (new SparseQR_rep (a, order)) { } - - SparseQR (const SparseQR& a) : rep (a.rep) { rep->count++; } - - ~SparseQR (void) - { - if (--rep->count == 0) - delete rep; - } - - SparseQR& operator = (const SparseQR& a) - { - if (this != &a) - { - if (--rep->count == 0) - delete rep; - - rep = a.rep; - rep->count++; - } - return *this; - } - - bool ok (void) const { return rep->ok (); } - - SparseMatrix V (void) const { return rep->V (); } - - ColumnVector Pinv (void) const { return rep->P (); } - - ColumnVector P (void) const { return rep->P (); } - - SparseMatrix R (const bool econ = false) const { return rep->R(econ); } - - Matrix C (const Matrix &b) const { return rep->C(b); } - - Matrix Q (void) const { return rep->Q (); } - - friend Matrix qrsolve (const SparseMatrix &a, const MArray<double> &b, - octave_idx_type &info); - - friend SparseMatrix qrsolve (const SparseMatrix &a, const SparseMatrix &b, - octave_idx_type &info); - - friend ComplexMatrix qrsolve (const SparseMatrix &a, const MArray<Complex> &b, - octave_idx_type &info); - - friend SparseComplexMatrix qrsolve (const SparseMatrix &a, - const SparseComplexMatrix &b, - octave_idx_type &info); - -protected: -#ifdef HAVE_CXSPARSE - CXSPARSE_DNAME (s) * S (void) { return rep->S; } - - CXSPARSE_DNAME (n) * N (void) { return rep->N; } -#endif -}; - - -// Publish externally used friend functions. - -extern Matrix qrsolve (const SparseMatrix &a, const MArray<double> &b, - octave_idx_type &info); - -extern SparseMatrix qrsolve (const SparseMatrix &a, const SparseMatrix &b, - octave_idx_type &info); - -extern ComplexMatrix qrsolve (const SparseMatrix &a, const MArray<Complex> &b, - octave_idx_type &info); - -extern SparseComplexMatrix qrsolve (const SparseMatrix &a, - const SparseComplexMatrix &b, - octave_idx_type &info); - -#endif
--- a/liboctave/numeric/module.mk Tue Feb 02 12:33:38 2016 -0500 +++ b/liboctave/numeric/module.mk Tue Feb 02 12:34:06 2016 -0500 @@ -81,8 +81,7 @@ liboctave/numeric/randpoisson.h \ liboctave/numeric/sparse-chol.h \ liboctave/numeric/sparse-lu.h \ - liboctave/numeric/SparseCmplxQR.h \ - liboctave/numeric/SparseQR.h + liboctave/numeric/sparse-qr.h NUMERIC_C_SRC = \ liboctave/numeric/randgamma.c \ @@ -142,11 +141,10 @@ liboctave/numeric/oct-spparms.cc \ liboctave/numeric/ODES.cc \ liboctave/numeric/Quad.cc \ - liboctave/numeric/SparseCmplxQR.cc \ - liboctave/numeric/SparseQR.cc \ liboctave/numeric/sparse-chol-inst.cc \ liboctave/numeric/sparse-lu-inst.cc \ -$(NUMERIC_C_SRC) + liboctave/numeric/sparse-qr-inst.cc \ + $(NUMERIC_C_SRC) LIBOCTAVE_TEMPLATE_SRC += \ liboctave/numeric/base-lu.cc \ @@ -155,6 +153,7 @@ liboctave/numeric/eigs-base.cc \ liboctave/numeric/sparse-chol.cc \ liboctave/numeric/sparse-lu.cc \ + liboctave/numeric/sparse-qr.cc \ liboctave/numeric/sparse-dmsolve.cc ## Special rules for sources which must be built before rest of compilation.
--- a/liboctave/numeric/sparse-dmsolve.cc Tue Feb 02 12:33:38 2016 -0500 +++ b/liboctave/numeric/sparse-dmsolve.cc Tue Feb 02 12:34:06 2016 -0500 @@ -28,8 +28,7 @@ #include "MArray.h" #include "MSparse.h" -#include "SparseQR.h" -#include "SparseCmplxQR.h" +#include "sparse-qr.h" #include "MatrixType.h" #include "oct-sort.h" #include "oct-locbuf.h"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/sparse-qr-inst.cc Tue Feb 02 12:34:06 2016 -0500 @@ -0,0 +1,32 @@ +/* + +Copyright (C) 2016 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "sparse-qr.h" +#include "sparse-qr.cc" + +template class sparse_qr<SparseMatrix>; + +template class sparse_qr<SparseComplexMatrix>;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/sparse-qr.cc Tue Feb 02 12:34:06 2016 -0500 @@ -0,0 +1,1975 @@ +/* + +Copyright (C) 2005-2015 David Bateman +Copyright (C) 2016 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "lo-error.h" +#include "sparse-qr.h" +#include "oct-locbuf.h" + +template <typename SPARSE_T> +class +cxsparse_types +{ +public: + typedef void symbolic_type; + typedef void numeric_type; +}; + +template <> +class +cxsparse_types<SparseMatrix> +{ +public: +#ifdef HAVE_CXSPARSE + typedef CXSPARSE_DNAME (s) symbolic_type; + typedef CXSPARSE_DNAME (n) numeric_type; +#endif +}; + +template <> +class +cxsparse_types<SparseComplexMatrix> +{ +public: +#ifdef HAVE_CXSPARSE + typedef CXSPARSE_ZNAME (s) symbolic_type; + typedef CXSPARSE_ZNAME (n) numeric_type; +#endif +}; + +template <typename SPARSE_T> +class sparse_qr<SPARSE_T>::sparse_qr_rep +{ +public: + + sparse_qr_rep (const SPARSE_T& a, int order); + + ~sparse_qr_rep (void); + + bool ok (void) const + { +#ifdef HAVE_CXSPARSE + return (N && S); +#else + return false; +#endif + } + + SPARSE_T V (void) const; + + ColumnVector Pinv (void) const; + + ColumnVector P (void) const; + + SPARSE_T R (bool econ) const; + + typename SPARSE_T::dense_matrix_type + C (const typename SPARSE_T::dense_matrix_type& b) const; + + typename SPARSE_T::dense_matrix_type + Q (void) const; + + octave_refcount<int> count; + + octave_idx_type nrows; + octave_idx_type ncols; + + typename cxsparse_types<SPARSE_T>::symbolic_type *S; + typename cxsparse_types<SPARSE_T>::numeric_type *N; + + template <typename RHS_T, typename RET_T> + RET_T + tall_solve (const RHS_T& b, octave_idx_type& info) const; + + template <typename RHS_T, typename RET_T> + RET_T + wide_solve (const RHS_T& b, octave_idx_type& info) const; + +private: + + // No copying! + + sparse_qr_rep (const sparse_qr_rep&); + + sparse_qr_rep& operator = (const sparse_qr_rep&); +}; + +template <typename SPARSE_T> +ColumnVector +sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const +{ +#ifdef HAVE_CXSPARSE + ColumnVector ret (N->L->m); + for (octave_idx_type i = 0; i < N->L->m; i++) + ret.xelem (i) = S->pinv[i]; + return ret; +#else + return ColumnVector (); +#endif +} + +template <typename SPARSE_T> +ColumnVector +sparse_qr<SPARSE_T>::sparse_qr_rep::P (void) const +{ +#ifdef HAVE_CXSPARSE + ColumnVector ret (N->L->m); + for (octave_idx_type i = 0; i < N->L->m; i++) + ret.xelem (S->pinv[i]) = i; + return ret; +#else + return ColumnVector (); +#endif +} + +// Specializations. + +// Real-valued matrices. + +template <> +sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep + (const SparseMatrix& a, int order) + : count (1), nrows (a.rows ()), ncols (a.columns ()) +#ifdef HAVE_CXSPARSE + , S (0), N (0) +#endif +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_DNAME () A; + A.nzmax = a.nnz (); + A.m = nrows; + A.n = ncols; + // Cast away const on A, with full knowledge that CSparse won't touch it + // Prevents the methods below making a copy of the data. + A.p = const_cast<octave_idx_type *>(a.cidx ()); + A.i = const_cast<octave_idx_type *>(a.ridx ()); + A.x = const_cast<double *>(a.data ()); + A.nz = -1; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + S = CXSPARSE_DNAME (_sqr) (order, &A, 1); + N = CXSPARSE_DNAME (_qr) (&A, S); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + if (! N) + (*current_liboctave_error_handler) + ("sparse_qr: sparse matrix QR factorization filled"); + + count = 1; + +#else + (*current_liboctave_error_handler) + ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built"); +#endif +} + +template <> +sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep (void) +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_DNAME (_sfree) (S); + CXSPARSE_DNAME (_nfree) (N); +#endif +} + +template <> +SparseMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::V (void) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from V and sort + // FIXME: Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_dropzeros) (N->L); + CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->L, 1); + CXSPARSE_DNAME (_spfree) (N->L); + N->L = CXSPARSE_DNAME (_transpose) (D, 1); + CXSPARSE_DNAME (_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->L->n; + octave_idx_type nz = N->L->nzmax; + SparseMatrix ret (N->L->m, nc, nz); + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->L->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->L->i[j]; + ret.xdata (j) = N->L->x[j]; + } + return ret; +#else + return SparseMatrix (); +#endif +} + +template <> +SparseMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from R and sort + // FIXME: Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_dropzeros) (N->U); + CXSPARSE_DNAME () *D = CXSPARSE_DNAME (_transpose) (N->U, 1); + CXSPARSE_DNAME (_spfree) (N->U); + N->U = CXSPARSE_DNAME (_transpose) (D, 1); + CXSPARSE_DNAME (_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->U->n; + octave_idx_type nz = N->U->nzmax; + + SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); + + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->U->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->U->i[j]; + ret.xdata (j) = N->U->x[j]; + } + return ret; +#else + return SparseMatrix (); +#endif +} + +template <> +Matrix +sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + const double *bvec = b.fortran_vec (); + Matrix ret (b_nr, b_nc); + double *vec = ret.fortran_vec (); + if (nr < 0 || nc < 0 || nr != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + if (nr == 0 || nc == 0 || b_nc == 0) + ret = Matrix (nc, b_nc, 0.0); + else + { + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) + { + octave_quit (); + for (octave_idx_type i = nr; i < S->m2; i++) + buf[i] = 0.; + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + idx, buf, b_nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type i = 0; i < nm; i++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < b_nr; i++) + vec[i+idx] = buf[i]; + } + } + return ret; +#else + return Matrix (); +#endif +} + +template <> +Matrix +sparse_qr<SparseMatrix>::sparse_qr_rep::Q (void) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + Matrix ret (nr, nr); + double *vec = ret.fortran_vec (); + if (nr < 0 || nc < 0) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + if (nr == 0 || nc == 0) + ret = Matrix (nc, nr, 0.0); + else + { + OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1); + for (octave_idx_type i = 0; i < nr; i++) + bvec[i] = 0.; + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) + { + octave_quit (); + bvec[j] = 1.0; + for (octave_idx_type i = nr; i < S->m2; i++) + buf[i] = 0.; + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type i = 0; i < nm; i++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, i, N->B[i], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < nr; i++) + vec[i+idx] = buf[i]; + bvec[j] = 0.0; + } + } + return ret.transpose (); +#else + return Matrix (); +#endif +} + +template <> +template <> +Matrix +sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix> + (const MArray<double>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + const double *bvec = b.data (); + + Matrix x (nc, b_nc); + double *vec = x.fortran_vec (); + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + octave_quit (); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_usolve) (N->U, buf); + CXSPARSE_DNAME (_ipvec) (S->q, buf, vec + idx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + + info = 0; + + return x; + +#else + + return Matrix (); + +#endif +} + +template <> +template <> +Matrix +sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>, Matrix> + (const MArray<double>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + const double *bvec = b.data (); + + Matrix x (nc, b_nc); + double *vec = x.fortran_vec (); + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (double, buf, nbuf); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + octave_quit (); + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->q, bvec + bidx, buf, nr); + CXSPARSE_DNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->pinv, buf, vec + idx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + + info = 0; + + return x; + +#else + + return Matrix (); + +#endif +} + +template <> +template <> +SparseMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix> + (const SparseMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + volatile octave_idx_type ii, x_nz; + + SparseMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_usolve) (N->U, buf); + CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + double tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + return x; + +#else + + return SparseMatrix (); + +#endif +} + +template <> +template <> +SparseMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseMatrix> + (const SparseMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + volatile octave_idx_type ii, x_nz; + + SparseMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, nbuf); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr); + CXSPARSE_DNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + double tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + x.maybe_compress (); + + return x; + +#else + + return SparseMatrix (); + +#endif +} + +template <> +template <> +ComplexMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix> + (const MArray<Complex>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + ComplexMatrix x (nc, b_nc); + Complex *vec = x.fortran_vec (); + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_usolve) (N->U, buf); + CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = 0.; + CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_usolve) (N->U, buf); + CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = 0; j < nc; j++) + vec[j+idx] = Complex (Xx[j], Xz[j]); + } + + info = 0; + + return x; + +#else + + return ComplexMatrix (); + +#endif +} + +template <> +template <> +ComplexMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix> + (const MArray<Complex>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + ComplexMatrix x (nc, b_nc); + Complex *vec = x.fortran_vec (); + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, nbuf); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr); + CXSPARSE_DNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr); + CXSPARSE_DNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = 0; j < nc; j++) + vec[j+idx] = Complex (Xx[j], Xz[j]); + } + + info = 0; + + return x; + +#else + + return ComplexMatrix (); + +#endif +} + +// Complex-valued matrices. + +template <> +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::sparse_qr_rep + (const SparseComplexMatrix& a, int order) + : count (1), nrows (a.rows ()), ncols (a.columns ()) +#ifdef HAVE_CXSPARSE + , S (0), N (0) +#endif +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_ZNAME () A; + A.nzmax = a.nnz (); + A.m = nrows; + A.n = ncols; + // Cast away const on A, with full knowledge that CSparse won't touch it + // Prevents the methods below making a copy of the data. + A.p = const_cast<octave_idx_type *>(a.cidx ()); + A.i = const_cast<octave_idx_type *>(a.ridx ()); + A.x = const_cast<cs_complex_t *>(reinterpret_cast<const cs_complex_t *> (a.data ())); + A.nz = -1; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + S = CXSPARSE_ZNAME (_sqr) (order, &A, 1); + N = CXSPARSE_ZNAME (_qr) (&A, S); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + if (! N) + (*current_liboctave_error_handler) + ("sparse_qr: sparse matrix QR factorization filled"); + + count = 1; + +#else + (*current_liboctave_error_handler) + ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built"); +#endif +} + +template <> +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::~sparse_qr_rep (void) +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_ZNAME (_sfree) (S); + CXSPARSE_ZNAME (_nfree) (N); +#endif +} + +template <> +SparseComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::V (void) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from V and sort + // FIXME: Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_dropzeros) (N->L); + CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->L, 1); + CXSPARSE_ZNAME (_spfree) (N->L); + N->L = CXSPARSE_ZNAME (_transpose) (D, 1); + CXSPARSE_ZNAME (_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->L->n; + octave_idx_type nz = N->L->nzmax; + SparseComplexMatrix ret (N->L->m, nc, nz); + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->L->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->L->i[j]; + ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j]; + } + return ret; +#else + return SparseComplexMatrix (); +#endif +} + +template <> +SparseComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from R and sort + // FIXME: Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_dropzeros) (N->U); + CXSPARSE_ZNAME () *D = CXSPARSE_ZNAME (_transpose) (N->U, 1); + CXSPARSE_ZNAME (_spfree) (N->U); + N->U = CXSPARSE_ZNAME (_transpose) (D, 1); + CXSPARSE_ZNAME (_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->U->n; + octave_idx_type nz = N->U->nzmax; + + SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); + + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->U->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->U->i[j]; + ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j]; + } + return ret; +#else + return SparseComplexMatrix (); +#endif +} + +template <> +ComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C (const ComplexMatrix& b) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>(b.fortran_vec ()); + ComplexMatrix ret (b_nr, b_nc); + Complex *vec = ret.fortran_vec (); + if (nr < 0 || nc < 0 || nr != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + if (nr == 0 || nc == 0 || b_nc == 0) + ret = ComplexMatrix (nc, b_nc, Complex (0.0, 0.0)); + else + { + OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) + { + octave_quit (); + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + idx, reinterpret_cast<cs_complex_t *>(buf), b_nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type i = 0; i < nm; i++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < b_nr; i++) + vec[i+idx] = buf[i]; + } + } + return ret; +#else + return ComplexMatrix (); +#endif +} + +template <> +ComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (void) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + ComplexMatrix ret (nr, nr); + Complex *vec = ret.fortran_vec (); + if (nr < 0 || nc < 0) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + if (nr == 0 || nc == 0) + ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0)); + else + { + OCTAVE_LOCAL_BUFFER (cs_complex_t, bvec, nr); + for (octave_idx_type i = 0; i < nr; i++) + bvec[i] = cs_complex_t (0.0, 0.0); + OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) + { + octave_quit (); + bvec[j] = cs_complex_t (1.0, 0.0); + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type i = 0; i < nm; i++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < nr; i++) + vec[i+idx] = buf[i]; + bvec[j] = cs_complex_t (0.0, 0.0); + } + } + return ret.hermitian (); +#else + return ComplexMatrix (); +#endif +} + +template <> +template <> +SparseComplexMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix> + (const SparseComplexMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + volatile octave_idx_type ii, x_nz; + + SparseComplexMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, Xx, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_usolve) (N->U, buf); + CXSPARSE_DNAME (_ipvec) (S->q, buf, Xx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_ipvec) (S->pinv, Xz, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_usolve) (N->U, buf); + CXSPARSE_DNAME (_ipvec) (S->q, buf, Xz, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Complex (Xx[j], Xz[j]); + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + return x; + +#else + + return SparseComplexMatrix (); + +#endif +} + +template <> +template <> +SparseComplexMatrix +sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix> + (const SparseComplexMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + + volatile octave_idx_type ii, x_nz; + + SparseComplexMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, nbuf); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->q, Xx, buf, nr); + CXSPARSE_DNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = 0.; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->q, Xz, buf, nr); + CXSPARSE_DNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (_pvec) (S->pinv, buf, Xz, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Complex (Xx[j], Xz[j]); + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + x.maybe_compress (); + + return x; + +#else + + return SparseComplexMatrix (); + +#endif +} + +template <> +template <> +ComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>, ComplexMatrix> + (const MArray<double>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + ComplexMatrix x (nc, b_nc); + cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); + OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_ipvec) (S->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_usolve) (N->U, buf); + CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + + info = 0; + + return x; + +#else + + return ComplexMatrix (); + +#endif +} + +template <> +template <> +ComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<double>, ComplexMatrix> + (const MArray<double>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseComplexMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + ComplexMatrix x (nc, b_nc); + cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); + + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); + OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); + OCTAVE_LOCAL_BUFFER (double, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = N->B[i]; + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); + CXSPARSE_ZNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + + info = 0; + + return x; + +#else + + return ComplexMatrix (); + +#endif +} + +template <> +template <> +SparseComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseComplexMatrix> + (const SparseMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + volatile octave_idx_type ii, x_nz; + + SparseComplexMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_ipvec) (S->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_usolve) (N->U, buf); + CXSPARSE_ZNAME (_ipvec) (S->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + x.maybe_compress (); + + return x; + +#else + + return SparseComplexMatrix (); + +#endif +} + +template <> +template <> +SparseComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseMatrix, SparseComplexMatrix> + (const SparseMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseComplexMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + volatile octave_idx_type ii, x_nz; + + SparseComplexMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); + + OCTAVE_LOCAL_BUFFER (double, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = N->B[i]; + + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); + CXSPARSE_ZNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + x.maybe_compress (); + + return x; + +#else + + return SparseComplexMatrix (); + +#endif +} + +template <> +template <> +ComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix> + (const MArray<Complex>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>(b.fortran_vec ()); + + ComplexMatrix x (nc, b_nc); + cs_complex_t *vec = reinterpret_cast<cs_complex_t *> + (x.fortran_vec ()); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + octave_quit (); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_ipvec) (S->pinv, bvec + bidx, buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_usolve) (N->U, buf); + CXSPARSE_ZNAME (_ipvec) (S->q, buf, vec + idx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + + info = 0; + + return x; + +#else + + return ComplexMatrix (); + +#endif +} + +template <> +template <> +ComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix> + (const MArray<Complex>& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseComplexMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + const cs_complex_t *bvec = reinterpret_cast<const cs_complex_t *>(b.fortran_vec ()); + + ComplexMatrix x (nc, b_nc); + cs_complex_t *vec = reinterpret_cast<cs_complex_t *> (x.fortran_vec ()); + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); + OCTAVE_LOCAL_BUFFER (double, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = N->B[i]; + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + octave_quit (); + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->q, bvec + bidx, buf, nr); + CXSPARSE_ZNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->pinv, buf, vec + idx, nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + info = 0; + + return x; + +#else + + return ComplexMatrix (); + +#endif +} + +template <> +template <> +SparseComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix> + (const SparseComplexMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + octave_idx_type nr = nrows; + octave_idx_type nc = ncols; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + volatile octave_idx_type ii, x_nz; + + SparseComplexMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, S->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < S->m2; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_ipvec) (S->pinv, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, N->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_usolve) (N->U, buf); + CXSPARSE_ZNAME (_ipvec) (S->q, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + x.maybe_compress (); + + return x; + +#else + + return SparseComplexMatrix (); + +#endif +} + +template <> +template <> +SparseComplexMatrix +sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix> + (const SparseComplexMatrix& b, octave_idx_type& info) const +{ + info = -1; + +#ifdef HAVE_CXSPARSE + + // These are swapped because the original matrix was transposed in + // sparse_qr<SparseComplexMatrix>::solve. + + octave_idx_type nr = ncols; + octave_idx_type nc = nrows; + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + volatile octave_idx_type ii, x_nz; + + SparseComplexMatrix x (nc, b_nc, b.nnz ()); + x.xcidx (0) = 0; + x_nz = b.nnz (); + ii = 0; + volatile octave_idx_type nbuf = (nc > S->m2 ? nc : S->m2); + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, nbuf); + OCTAVE_LOCAL_BUFFER (double, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = N->B[i]; + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + octave_quit (); + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem (j,i); + for (octave_idx_type j = nr; j < nbuf; j++) + buf[j] = cs_complex_t (0.0, 0.0); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->q, reinterpret_cast<cs_complex_t *>(Xx), buf, nr); + CXSPARSE_ZNAME (_utsolve) (N->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + octave_quit (); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_happly) (N->L, j, B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (_pvec) (S->pinv, buf, reinterpret_cast<cs_complex_t *>(Xx), nc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata (ii) = tmp; + x.xridx (ii++) = j; + } + } + x.xcidx (i+1) = ii; + } + + info = 0; + + x.maybe_compress (); + + return x; + +#else + + return SparseComplexMatrix (); + +#endif +} + +template <typename SPARSE_T> +sparse_qr<SPARSE_T>::sparse_qr (void) + : rep (new sparse_qr_rep (SPARSE_T (), 0)) +{ } + +template <typename SPARSE_T> +sparse_qr<SPARSE_T>::sparse_qr (const SPARSE_T& a, int order) + : rep (new sparse_qr_rep (a, order)) +{ } + +template <typename SPARSE_T> +sparse_qr<SPARSE_T>::sparse_qr (const sparse_qr<SPARSE_T>& a) + : rep (a.rep) +{ + rep->count++; +} + +template <typename SPARSE_T> +sparse_qr<SPARSE_T>::~sparse_qr (void) +{ + if (--rep->count == 0) + delete rep; +} + +template <typename SPARSE_T> +sparse_qr<SPARSE_T>& +sparse_qr<SPARSE_T>::operator = (const sparse_qr<SPARSE_T>& a) +{ + if (this != &a) + { + if (--rep->count == 0) + delete rep; + + rep = a.rep; + rep->count++; + } + + return *this; +} + +template <typename SPARSE_T> +bool +sparse_qr<SPARSE_T>::ok (void) const +{ + return rep->ok (); +} + +template <typename SPARSE_T> +SPARSE_T +sparse_qr<SPARSE_T>::V (void) const +{ + return rep->V (); +} + +template <typename SPARSE_T> +ColumnVector +sparse_qr<SPARSE_T>::Pinv (void) const +{ + return rep->P (); +} + +template <typename SPARSE_T> +ColumnVector +sparse_qr<SPARSE_T>::P (void) const +{ + return rep->P (); +} + +template <typename SPARSE_T> +SPARSE_T +sparse_qr<SPARSE_T>::R (bool econ) const +{ + return rep->R (econ); +} + +template <typename SPARSE_T> +typename SPARSE_T::dense_matrix_type +sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b) const +{ + return rep->C (b); +} + +template <typename SPARSE_T> +typename SPARSE_T::dense_matrix_type +sparse_qr<SPARSE_T>::Q (void) const +{ + return rep->Q (); +} + +// FIXME: Why is the "order" of the QR calculation as used in the +// CXSparse function sqr 3 for real matrices and 2 for complex? These +// values seem to be required but there was no explanation in David +// Bateman's original code. + +template <typename SPARSE_T> +class +cxsparse_defaults +{ +public: + enum { order = -1 }; +}; + +template <> +class +cxsparse_defaults<SparseMatrix> +{ +public: + enum { order = 3 }; +}; + +template <> +class +cxsparse_defaults<SparseComplexMatrix> +{ +public: + enum { order = 2 }; +}; + +template <typename SPARSE_T> +template <typename RHS_T, typename RET_T> +RET_T +sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b, + octave_idx_type& info) +{ + info = -1; + + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); + + octave_idx_type b_nc = b.cols (); + octave_idx_type b_nr = b.rows (); + + int order = cxsparse_defaults<SPARSE_T>::order; + + if (nr < 0 || nc < 0 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + + if (nr == 0 || nc == 0 || b_nc == 0) + { + info = 0; + + return RET_T (nc, b_nc, 0.0); + } + else if (nr >= nc) + { + sparse_qr<SPARSE_T> q (a, order); + + return q.ok () ? q.tall_solve<RHS_T, RET_T> (b, info) : RET_T (); + } + else + { + sparse_qr<SPARSE_T> q (a.hermitian (), order); + + return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T (); + } +} + +template <typename SPARSE_T> +template <typename RHS_T, typename RET_T> +RET_T +sparse_qr<SPARSE_T>::tall_solve (const RHS_T& b, octave_idx_type& info) const +{ + return rep->tall_solve<RHS_T, RET_T> (b, info); +} + +template <typename SPARSE_T> +template <typename RHS_T, typename RET_T> +RET_T +sparse_qr<SPARSE_T>::wide_solve (const RHS_T& b, octave_idx_type& info) const +{ + return rep->wide_solve<RHS_T, RET_T> (b, info); +} + +Matrix +qrsolve (const SparseMatrix& a, const MArray<double>& b, + octave_idx_type& info) +{ + return sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix> (a, b, info); +} + +SparseMatrix +qrsolve (const SparseMatrix& a, const SparseMatrix& b, + octave_idx_type& info) +{ + return sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix> (a, b, info); +} + +ComplexMatrix +qrsolve (const SparseMatrix& a, const MArray<Complex>& b, + octave_idx_type& info) +{ + return sparse_qr<SparseMatrix>::solve<MArray<Complex>, ComplexMatrix> (a, b, info); +} + +SparseComplexMatrix +qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info) +{ + return sparse_qr<SparseMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix> (a, b, info); +} + +ComplexMatrix +qrsolve (const SparseComplexMatrix& a, const MArray<double>& b, + octave_idx_type& info) +{ + return sparse_qr<SparseComplexMatrix>::solve<MArray<double>, ComplexMatrix> (a, b, info); +} + +SparseComplexMatrix +qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b, + octave_idx_type& info) +{ + return sparse_qr<SparseComplexMatrix>::solve<SparseMatrix, SparseComplexMatrix> (a, b, info); +} + +ComplexMatrix +qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b, + octave_idx_type& info) +{ + return sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>, ComplexMatrix> (a, b, info); +} + +SparseComplexMatrix +qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info) +{ + return sparse_qr<SparseComplexMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix> (a, b, info); +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/numeric/sparse-qr.h Tue Feb 02 12:34:06 2016 -0500 @@ -0,0 +1,122 @@ +/* + +Copyright (C) 2005-2015 David Bateman +Copyright (C) 2016 John W. Eaton + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#if ! defined (octave_sparse_qr_h) +#define octave_sparse_qr_h 1 + +#include "dMatrix.h" +#include "CMatrix.h" +#include "dSparse.h" +#include "CSparse.h" +#include "oct-sparse.h" + +template <typename SPARSE_T> +class +sparse_qr +{ +public: + + sparse_qr (void); + + sparse_qr (const SPARSE_T& a, int order = 0); + + sparse_qr (const sparse_qr& a); + + ~sparse_qr (void); + + sparse_qr& operator = (const sparse_qr& a); + + bool ok (void) const; + + SPARSE_T V (void) const; + + ColumnVector Pinv (void) const; + + ColumnVector P (void) const; + + SPARSE_T R (bool econ = false) const; + + typename SPARSE_T::dense_matrix_type + C (const typename SPARSE_T::dense_matrix_type& b) const; + + typename SPARSE_T::dense_matrix_type + Q (void) const; + + template <typename RHS_T, typename RET_T> + static RET_T + solve (const SPARSE_T& a, const RHS_T& b, + octave_idx_type& info); + +private: + + class sparse_qr_rep; + + sparse_qr_rep *rep; + + template <typename RHS_T, typename RET_T> + RET_T + tall_solve (const RHS_T& b, octave_idx_type& info) const; + + template <typename RHS_T, typename RET_T> + RET_T + wide_solve (const RHS_T& b, octave_idx_type& info) const; +}; + +// Provide qrsolve for backward compatibility. + +extern Matrix +qrsolve (const SparseMatrix& a, const MArray<double>& b, + octave_idx_type& info); + +extern SparseMatrix +qrsolve (const SparseMatrix& a, const SparseMatrix& b, + octave_idx_type& info); + +extern ComplexMatrix +qrsolve (const SparseMatrix& a, const MArray<Complex>& b, + octave_idx_type& info); + +extern SparseComplexMatrix +qrsolve (const SparseMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info); + +extern ComplexMatrix +qrsolve (const SparseComplexMatrix& a, const MArray<double>& b, + octave_idx_type& info); + +extern SparseComplexMatrix +qrsolve (const SparseComplexMatrix& a, const SparseMatrix& b, + octave_idx_type& info); + +extern ComplexMatrix +qrsolve (const SparseComplexMatrix& a, const MArray<Complex>& b, + octave_idx_type& info); + +extern SparseComplexMatrix +qrsolve (const SparseComplexMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info); + +typedef sparse_qr<SparseMatrix> SparseQR; +typedef sparse_qr<SparseComplexMatrix> SparseComplexQR; + +#endif
--- a/liboctave/util/oct-sparse.h Tue Feb 02 12:33:38 2016 -0500 +++ b/liboctave/util/oct-sparse.h Tue Feb 02 12:34:06 2016 -0500 @@ -107,4 +107,14 @@ # endif #endif +#if defined (HAVE_CXSPARSE) +# if defined (ENABLE_64) +# define CXSPARSE_DNAME(name) cs_dl ## name +# define CXSPARSE_ZNAME(name) cs_cl ## name +# else +# define CXSPARSE_DNAME(name) cs_di ## name +# define CXSPARSE_ZNAME(name) cs_ci ## name +# endif #endif + +#endif