# HG changeset patch # User Simon Hau # Date 1578602754 -3600 # Thu Jan 09 21:45:54 2020 +0100 # Node ID 385b07c4ee70f43fbdc0472ffa0fc606978594dd # Parent 4ab10bf273405bc9425cf0a54f4f5f2c058bfe7f Replacement of CXSPARSE by SPQR for QR factorization (bug #57033). * libinterp/corefcn/qr.cc (Fqr): Handle sparse input. Return economy versions of C, Q and R for qr(..., 0). Return column permutation of A as (vector|matrix) for qr(..., ("vector"|"matrix")). Add documentation for sparse input A. Additional tests if HAVE_SPQR. * liboctave/numeric/oct-sparse.h, liboctave/numeric/oct-sparse.cc (ros2rcs, cos2ccs, rod2ccd, rod2rcd, cod2ccd, rcs2ros, ccs2cos, ros2ccs): New functions to convert between Octave and SuiteSparse matrix types. (spqr_error_handler): New function to handle SPQR errors. * liboctave/numeric/sparse-qr.h (octave::math::sparse_qr): Add new member functions. * liboctave/numeric/sparse-qr.cc: Add interface to SPQR functios. * configure.ac: Add configure tests for SPQR. diff -r 4ab10bf27340 -r 385b07c4ee70 configure.ac --- a/configure.ac Mon Jan 04 10:23:28 2021 +0100 +++ b/configure.ac Thu Jan 09 21:45:54 2020 +0100 @@ -2144,6 +2144,16 @@ [], [don't use CHOLMOD library, disable some sparse matrix functionality]) LIBS="$save_LIBS" +### Check for SPQR library +### FIXME: We actually need to check for SuiteSparseQR_min2norm. But +### the macro cannot handle templated functions. + +OCTAVE_CHECK_LIB(spqr, SPQR, + [SPQR library not found. This will result in some lack of functionality for sparse matrices.], + [suitesparse/SuiteSparseQR.hpp], + [SuiteSparseQR_C], + [C++], [don't use SPQR library, disable some sparse matrix functionality]) + ### Check for CXSparse library OCTAVE_CHECK_LIB(cxsparse, CXSparse, @@ -2209,11 +2219,11 @@ ## Order matters, at least on some systems (Cygwin, for example). -SPARSE_XCPPFLAGS="$CHOLMOD_CPPFLAGS $UMFPACK_CPPFLAGS $AMD_CPPFLAGS $CAMD_CPPFLAGS $COLAMD_CPPFLAGS $CCOLAMD_CPPFLAGS $CXSPARSE_CPPFLAGS" - -SPARSE_XLDFLAGS="$CHOLMOD_LDFLAGS $UMFPACK_LDFLAGS $AMD_LDFLAGS $CAMD_LDFLAGS $COLAMD_LDFLAGS $CCOLAMD_LDFLAGS $CXSPARSE_LDFLAGS" - -SPARSE_XLIBS="$CHOLMOD_LIBS $UMFPACK_LIBS $AMD_LIBS $CAMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $CXSPARSE_LIBS $SUITESPARSECONFIG_LIBS" +SPARSE_XCPPFLAGS="$CHOLMOD_CPPFLAGS $UMFPACK_CPPFLAGS $AMD_CPPFLAGS $CAMD_CPPFLAGS $COLAMD_CPPFLAGS $CCOLAMD_CPPFLAGS $CXSPARSE_CPPFLAGS $SPQR_CPPFLAGS" + +SPARSE_XLDFLAGS="$CHOLMOD_LDFLAGS $UMFPACK_LDFLAGS $AMD_LDFLAGS $CAMD_LDFLAGS $COLAMD_LDFLAGS $CCOLAMD_LDFLAGS $CXSPARSE_LDFLAGS $SPQR_LDFLAGS" + +SPARSE_XLIBS="$CHOLMOD_LIBS $UMFPACK_LIBS $AMD_LIBS $CAMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $CXSPARSE_LIBS $SUITESPARSECONFIG_LIBS $SPQR_LIBS" AC_SUBST(SPARSE_XCPPFLAGS) AC_SUBST(SPARSE_XLDFLAGS) @@ -3148,6 +3158,9 @@ Sndfile CPPFLAGS: $SNDFILE_CPPFLAGS Sndfile LDFLAGS: $SNDFILE_LDFLAGS Sndfile libraries: $SNDFILE_LIBS + SPQR CPPFLAGS: $SPQR_CPPFLAGS + SPQR LDFLAGS: $SPQR_LDFLAGS + SPQR libraries: $SPQR_LIBS SuiteSparse config libraries: $SUITESPARSECONFIG_LIBS SUNDIALS IDA CPPFLAGS: $SUNDIALS_IDA_CPPFLAGS SUNDIALS IDA LDFLAGS: $SUNDIALS_IDA_LDFLAGS diff -r 4ab10bf27340 -r 385b07c4ee70 libinterp/corefcn/qr.cc --- a/libinterp/corefcn/qr.cc Mon Jan 04 10:23:28 2021 +0100 +++ b/libinterp/corefcn/qr.cc Thu Jan 09 21:45:54 2020 +0100 @@ -79,6 +79,8 @@ return octave::math::qr::std; } +// dense X +// // [Q, R] = qr (X): form Q unitary and R upper triangular such // that Q * R = X // @@ -95,13 +97,21 @@ // // qr (X) alone returns the output of the LAPACK routine dgeqrf, such // that R = triu (qr (X)) - +// +// sparse X +// +// X = qr (A, B): if M < N, X is the minimum 2-norm solution of +// A\B. If M >= N, X is the least squares +// approximation of A\B. X is calculated by +// SPQR-function SuiteSparseQR_min2norm. +// DEFUN (qr, args, nargout, doc: /* -*- texinfo -*- @deftypefn {} {[@var{Q}, @var{R}] =} qr (@var{A}) -@deftypefnx {} {[@var{Q}, @var{R}, @var{P}] =} qr (@var{A}) # non-sparse A +@deftypefnx {} {[@var{Q}, @var{R}, @var{P}] =} qr (@var{A}) @deftypefnx {} {@var{X} =} qr (@var{A}) # non-sparse A @deftypefnx {} {@var{R} =} qr (@var{A}) # sparse A +@deftypefnx {} {@var{X} =} qr (@var{A}, @var{B}) # sparse A @deftypefnx {} {[@var{C}, @var{R}] =} qr (@var{A}, @var{B}) @deftypefnx {} {[@dots{}] =} qr (@dots{}, 0) @deftypefnx {} {[@dots{}] =} qr (@dots{}, "vector") @@ -164,8 +174,8 @@ @var{A} is full. (Note: unlike most commands, the single return value is not the first return value when multiple values are requested.) -If the matrix @var{A} is full, and a third output @var{P} is requested, then -@code{qr} calculates the permuted QR@tie{}factorization +If a third output @var{P} is requested, then @code{qr} calculates the permuted +QR@tie{}factorization @tex $QR = AP$ where $Q$ is an orthogonal matrix, $R$ is upper triangular, and $P$ is a permutation matrix. @@ -181,9 +191,14 @@ @var{P} is a permutation matrix. @end ifnottex -The permuted QR@tie{}factorization has the additional property that the -diagonal entries of @var{R} are ordered by decreasing magnitude. In other -words, @code{abs (diag (@var{R}))} will be ordered from largest to smallest. +If @var{A} is dense, the permuted QR@tie{}factorization has the additional +property that the diagonal entries of @var{R} are ordered by decreasing +magnitude. In other words, @code{abs (diag (@var{R}))} will be ordered +from largest to smallest. + +If @var{A} is sparse, @var{P} is a fill-reducing ordering of the columns +of @var{A}. In that case, the diagonal entries of @var{R} are not ordered by +decreasing magnitude. For example, given the matrix @code{@var{A} = [1, 2; 3, 4]}, @@ -213,27 +228,38 @@ @end group @end example -If the input matrix @var{A} is sparse then the sparse QR@tie{}factorization -is computed using @sc{CSparse}. Because the matrix @var{Q} is, in general, a -full matrix, it is recommended to request only one return value @var{R}. In -that case, the computation avoids the construction of @var{Q} and returns -@var{R} such that @code{@var{R} = chol (@var{A}' * @var{A})}. +If the input matrix @var{A} is sparse, the sparse QR@tie{}factorization +is computed by using @sc{SPQR} (or @sc{CSparse} if @sc{SPQR} is not available). +Because the matrix @var{Q} is, in general, a full matrix, it is recommended to +request only one return value @var{R}. In that case, the computation avoids +the construction of @var{Q} and returns a sparse @var{R} such that +@code{@var{R} = chol (@var{A}' * @var{A})}. -If an additional matrix @var{B} is supplied and two return values are -requested, then @code{qr} returns @var{C}, where +If @var{A} is dense, an additional matrix @var{B} is supplied and two +return values are requested, then @code{qr} returns @var{C}, where @code{@var{C} = @var{Q}' * @var{B}}. This allows the least squares approximation of @code{@var{A} \ @var{B}} to be calculated as @example @group [@var{C}, @var{R}] = qr (@var{A}, @var{B}) -@var{x} = @var{R} \ @var{C} +@var{X} = @var{R} \ @var{C} @end group @end example +If @var{A} is a sparse MxN matrix and an additional matrix @var{B} is +supplied, one or two return values are possible. If one return value @var{X} +is requested and M < N, then @var{X} is the minimum 2-norm solution of +@w{@code{@var{A} \ @var{B}}}. If M >= N, @var{X} is the least squares +approximation @w{of @code{@var{A} \ @var{B}}}. If two return values are +requested, @var{C} and @var{R} have the same meaning as in the dense case +(@var{C} is dense and @var{R} is sparse). +The version with one return parameter should be preferred because +it uses less memory and can handle rank-deficient matrices better. + If the final argument is the string @qcode{"vector"} then @var{P} is a -permutation vector (of the columns of @var{A}) instead of a permutation matrix. -In this case, the defining relationship is +permutation vector (of the columns of @var{A}) instead of a permutation +matrix. In this case, the defining relationship is: @example @var{Q} * @var{R} = @var{A}(:, @var{P}) @@ -243,12 +269,13 @@ explicitly specified by using a final argument of @qcode{"matrix"}. If the final argument is the scalar 0 an @qcode{"economy"} factorization is -returned. When the original matrix @var{A} has size MxN and M > N then the +returned. If the original matrix @var{A} has size MxN and M > N, then the @qcode{"economy"} factorization will calculate just N rows in @var{R} and N -columns in @var{Q} and omit the zeros in @var{R}. If M @leq{} N there is no +columns in @var{Q} and omit the zeros in @var{R}. If M @leq{} N, there is no difference between the economy and standard factorizations. When calculating -an @qcode{"economy"} factorization the output @var{P} is always a vector -rather than a matrix. +an @qcode{"economy"} factorization and @var{A} is dense, the output @var{P} is +always a vector rather than a matrix. If @var{A} is sparse, output +@var{P} is a sparse permutation matrix. Background: The QR factorization has applications in the solution of least squares problems @@ -331,45 +358,159 @@ if (arg.issparse ()) { - if (nargout > 2) - error ("qr: Permutation output is not supported for sparse input"); + if (nargout > 3) + error ("qr: too many output arguments"); if (is_cmplx) { - octave::math::sparse_qr q (arg.sparse_complex_matrix_value ()); + if (have_b && nargout == 1) + { + octave_idx_type info; - if (have_b) + if (! args(1).issparse () && args(1).iscomplex ()) + retval = ovl + (octave::math::sparse_qr::solve + , ComplexMatrix> + (arg.sparse_complex_matrix_value (), + args(1).complex_matrix_value (), info)); + else if (args(1).issparse () && args(1).iscomplex ()) + retval = ovl + (octave::math::sparse_qr::solve + + (arg.sparse_complex_matrix_value (), + args(1).sparse_complex_matrix_value (), info)); + else if (! args(1).issparse () && ! args(1).iscomplex ()) + retval = ovl + (octave::math::sparse_qr::solve + , ComplexMatrix> + (arg.sparse_complex_matrix_value (), + args(1).matrix_value (), info)); + else if (args(1).issparse () && ! args(1).iscomplex ()) + retval = ovl + (octave::math::sparse_qr::solve + + (arg.sparse_complex_matrix_value (), + args(1).sparse_matrix_value (), info)); + else + error ("qr: b is not valid"); + } + else if (have_b && nargout == 2) { - retval = ovl (q.C (args(1).complex_matrix_value ()), + octave::math::sparse_qr + q (arg.sparse_complex_matrix_value (), 0); + retval = ovl (q.C (args(1).complex_matrix_value (), economy), q.R (economy)); - if (arg.rows () < arg.columns ()) - warning ("qr: non minimum norm solution for under-determined " - "problem %" OCTAVE_IDX_TYPE_FORMAT - "x%" OCTAVE_IDX_TYPE_FORMAT, - arg.rows (), arg.columns ()); + } + else if (have_b && nargout == 3) + { + octave::math::sparse_qr + q (arg.sparse_complex_matrix_value ()); + if (vector_p) + retval = ovl (q.C (args(1).complex_matrix_value (), economy), + q.R (economy), q.E ()); + else + retval = ovl (q.C (args(1).complex_matrix_value (), economy), + q.R (economy), q.E_MAT ()); } - else if (nargout > 1) - retval = ovl (q.Q (), q.R (economy)); else - retval = ovl (q.R (economy)); + { + if (nargout > 2) + { + octave::math::sparse_qr + q (arg.sparse_complex_matrix_value ()); + if (vector_p) + retval = ovl (q.Q (economy), q.R (economy), q.E ()); + else + retval = ovl (q.Q (economy), q.R (economy), + q.E_MAT ()); + } + else if (nargout > 1) + { + octave::math::sparse_qr + q (arg.sparse_complex_matrix_value (), 0); + retval = ovl (q.Q (economy), q.R (economy)); + } + else + { + octave::math::sparse_qr + q (arg.sparse_complex_matrix_value (), 0); + retval = ovl (q.R (economy)); + } + } } else { - octave::math::sparse_qr q (arg.sparse_matrix_value ()); - - if (have_b) + if (have_b && nargout == 1) + { + octave_idx_type info; + if (args(1).issparse () && ! args(1).iscomplex ()) + retval = ovl (octave::math::sparse_qr::solve + + (arg.sparse_matrix_value (), + args (1).sparse_matrix_value (), info)); + else if (! args(1).issparse () && args(1).iscomplex ()) + retval = ovl (octave::math::sparse_qr::solve + , ComplexMatrix> + (arg.sparse_matrix_value (), + args (1).complex_matrix_value (), info)); + else if (! args(1).issparse () && ! args(1).iscomplex ()) + retval = ovl (octave::math::sparse_qr::solve + , Matrix> + (arg.sparse_matrix_value (), + args (1).matrix_value (), info)); + else if (args(1).issparse () && args(1).iscomplex ()) + retval = ovl (octave::math::sparse_qr::solve + + (arg.sparse_matrix_value (), + args(1).sparse_complex_matrix_value (), + info)); + else + error ("qr: b is not valid"); + } + else if (have_b && nargout == 2) + { + octave::math::sparse_qr + q (arg.sparse_matrix_value (), 0); + retval = ovl (q.C (args(1).matrix_value (), economy), + q.R (economy)); + } + else if (have_b && nargout == 3) { - retval = ovl (q.C (args(1).matrix_value ()), q.R (economy)); - if (arg.rows () < arg.columns ()) - warning ("qr: non minimum norm solution for under-determined " - "problem %" OCTAVE_IDX_TYPE_FORMAT - "x%" OCTAVE_IDX_TYPE_FORMAT, - arg.rows (), arg.columns ()); + octave::math::sparse_qr + q (arg.sparse_matrix_value ()); + if (vector_p) + retval = ovl (q.C (args(1).matrix_value (), economy), + q.R (economy), q.E ()); + else + retval = ovl (q.C (args(1).matrix_value (), economy), + q.R (economy), q.E_MAT ()); } - else if (nargout > 1) - retval = ovl (q.Q (), q.R (economy)); + else - retval = ovl (q.R (economy)); + { + if (nargout > 2) + { + octave::math::sparse_qr + q (arg.sparse_matrix_value ()); + if (vector_p) + retval = ovl (q.Q (economy), q.R (economy), q.E ()); + else + retval = ovl (q.Q (economy), q.R (economy), + q.E_MAT ()); + } + else if (nargout > 1) + { + octave::math::sparse_qr + q (arg.sparse_matrix_value (), 0); + retval = ovl (q.Q (economy), q.R (economy)); + } + else + { + octave::math::sparse_qr + q (arg.sparse_matrix_value (), 0); + retval = ovl (q.R (economy)); + } + } } } else @@ -443,8 +584,8 @@ octave::math::qr fact (m, type); retval = ovl (fact.Q (), get_qr_r (fact)); if (have_b) - retval (0) = conj (fact.Q ().transpose ()) - * args(1).float_complex_matrix_value (); + retval(0) = conj (fact.Q ().transpose ()) + * args(1).float_complex_matrix_value (); } break; @@ -528,8 +669,8 @@ octave::math::qr fact (m, type); retval = ovl (fact.Q (), get_qr_r (fact)); if (have_b) - retval (0) = conj (fact.Q ().transpose ()) - * args(1).complex_matrix_value (); + retval(0) = conj (fact.Q ().transpose ()) + * args(1).complex_matrix_value (); } break; @@ -952,10 +1093,150 @@ %! rand ("state", 42); %! randn ("state", 42); %! a = 1i*sprandn (n,n+1,d) + speye (n,n+1); -%! b = randn (n,2); -%! [c,r] = qr (a, b); +%! b = randn (n, 2); +%! [c, r] = qr (a, b); +%! assert (r\c, full (a)\b, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = randn (m, 2); +%! [c, r] = qr (a, b); %! assert (r\c, full (a)\b, 10e-10); +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = sprandn (m, 2, d); +%! [c, r] = qr (a, b, 0); +%! [c2, r2] = qr (full (a), full (b), 0); +%! assert (r\c, r2\c2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = randn (m, 2); +%! [c, r, p] = qr (a, b, "matrix"); +%! x = p * (r\c); +%! [c2, r2] = qr (full (a), b); +%! x2 = r2 \ c2; +%! assert (x, x2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! [q, r, p] = qr (a, "matrix"); +%! assert (q * r, a * p, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = randn (m, 2); +%! x = qr (a, b); +%! [c2, r2] = qr (full (a), b); +%! assert (x, r2\c2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = i * randn (m, 2); +%! x = qr (a, b); +%! [c2, r2] = qr (full (a), b); +%! assert (x, r2\c2, 10e-10); + +%!#testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = i * randn (m, 2); +%! [c, r] = qr (a, b); +%! [c2, r2] = qr (full (a), b); +%! assert (r\c, r2\c2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = sprandn (m, n, d); +%! b = i * randn (m, 2); +%! [c, r, p] = qr (a, b, "matrix"); +%! x = p * (r\c); +%! [c2, r2] = qr (full (a), b); +%! x2 = r2 \ c2; +%! assert (x, x2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = i * sprandn (m, n, d); +%! b = sprandn (m, 2, d); +%! [c, r] = qr (a, b, 0); +%! [c2, r2] = qr (full (a), full (b), 0); +%! assert (r\c, r2\c2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = i * sprandn (m, n, d); +%! b = randn (m, 2); +%! [c, r, p] = qr (a, b, "matrix"); +%! x = p * (r\c); +%! [c2, r2] = qr (full (a), b); +%! x2 = r2 \ c2; +%! assert(x, x2, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = i * sprandn (m, n, d); +%! [q, r, p] = qr (a, "matrix"); +%! assert(q * r, a * p, 10e-10); + +%!testif HAVE_SPQR +%! n = 12; m = 20; d = 0.2; +%! ## initialize generators to make behavior reproducible +%! rand ("state", 42); +%! randn ("state", 42); +%! a = i * sprandn (m, n, d); +%! b = randn (m, 2); +%! x = qr (a, b); +%! [c2, r2] = qr (full (a), b); +%! assert (x, r2\c2, 10e-10); + +%!testif HAVE_SPQR +%! a = sparse (5, 6); +%! a(3,1) = 0.8; +%! a(2,2) = 1.4; +%! a(1,6) = -0.5; +%! r = qr (a); +%! assert (r'*r, a'*a, 10e-10); */ static diff -r 4ab10bf27340 -r 385b07c4ee70 liboctave/numeric/sparse-qr.cc --- a/liboctave/numeric/sparse-qr.cc Mon Jan 04 10:23:28 2021 +0100 +++ b/liboctave/numeric/sparse-qr.cc Thu Jan 09 21:45:54 2020 +0100 @@ -47,15 +47,18 @@ template class cxsparse_types - { - }; + { }; template <> class cxsparse_types { public: -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + typedef CXSPARSE_DNAME (s) symbolic_type; + typedef CXSPARSE_DNAME (n) numeric_type; + typedef Matrix dense_matrix_type; +#elif defined (HAVE_CXSPARSE) typedef CXSPARSE_DNAME (s) symbolic_type; typedef CXSPARSE_DNAME (n) numeric_type; #else @@ -69,7 +72,10 @@ cxsparse_types { public: -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + typedef CXSPARSE_ZNAME (s) symbolic_type; + typedef CXSPARSE_ZNAME (n) numeric_type; +#elif defined (HAVE_CXSPARSE) typedef CXSPARSE_ZNAME (s) symbolic_type; typedef CXSPARSE_ZNAME (n) numeric_type; #else @@ -95,7 +101,9 @@ bool ok (void) const { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + return (m_H && m_Htau && m_HPinv && m_R && m_E && &m_cc); +#elif defined (HAVE_CXSPARSE) return (N && S); #else return false; @@ -108,29 +116,62 @@ ColumnVector P (void) const; + ColumnVector E (void) const; + SPARSE_T R (bool econ) const; typename SPARSE_T::dense_matrix_type - C (const typename SPARSE_T::dense_matrix_type& b) const; + C (const typename SPARSE_T::dense_matrix_type& b); typename SPARSE_T::dense_matrix_type - Q (void) const; + C (const typename SPARSE_T::dense_matrix_type& b, bool econ); + + typename SPARSE_T::dense_matrix_type Q (void); + + typename SPARSE_T::dense_matrix_type Q (bool econ); refcount count; octave_idx_type nrows; octave_idx_type ncols; +#if defined (HAVE_SPQR) + typename cxsparse_types::symbolic_type *S; typename cxsparse_types::numeric_type *N; template RET_T - tall_solve (const RHS_T& b, octave_idx_type& info) const; + solve (const RHS_T& b, octave_idx_type& info) const; + +#elif defined (HAVE_CXSPARSE) + + typename cxsparse_types::symbolic_type *S; + typename cxsparse_types::numeric_type *N; + +#endif + + template + RET_T + tall_solve (const RHS_T& b, octave_idx_type& info); template RET_T wide_solve (const RHS_T& b, octave_idx_type& info) const; + +#if defined (HAVE_SPQR) + + private: + + cholmod_common m_cc; + cholmod_sparse *m_R; // R factor + // Column permutation for A. Fill-reducing ordering. + suitesparse_integer *m_E; + cholmod_sparse *m_H; // Householder vectors + cholmod_dense *m_Htau; // beta scalars + suitesparse_integer *m_HPinv; + +#endif }; template @@ -157,7 +198,17 @@ ColumnVector sparse_qr::sparse_qr_rep::P (void) const { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + ColumnVector ret (nrows); + + // FIXME: Is ret.xelem (m_HPinv[i]) = i + 1 correct? + for (octave_idx_type i = 0; i < nrows; i++) + ret.xelem (m_HPinv[i]) = i + 1; + + return ret; + +#elif defined (HAVE_CXSPARSE) ColumnVector ret (N->L->m); @@ -173,20 +224,76 @@ #endif } + template + ColumnVector + sparse_qr::sparse_qr_rep::E (void) const + { +#if defined (HAVE_SPQR) + + ColumnVector ret (ncols); + + for (octave_idx_type i = 0; i < ncols; i++) + ret(i) = static_cast (m_E[i]) + 1; + + return ret; + +#else + + return ColumnVector (); + +#endif + } + // Specializations. // Real-valued matrices. + // Arguments for parameter order (taken from SuiteSparseQR documentation). + // 0: fixed ordering 0 (no permutation of columns) + // 1: natural ordering 1 (only singleton columns are permuted to the left of + // the matrix) + // 2: colamd + // 3: + // 4: CHOLMOD best-effort (COLAMD, METIS,...) + // 5: AMD(a'*a) + // 6: metis(a'*a) + // 7: SuiteSparseQR default ordering + // 8: try COLAMD, AMD, and METIS; pick best + // 9: try COLAMD and AMD; pick best + //FIXME: What is order = 3? template <> sparse_qr::sparse_qr_rep::sparse_qr_rep (const SparseMatrix& a, int order) : count (1), nrows (a.rows ()), ncols (a.columns ()) -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + , m_cc (), m_R (nullptr), m_E (nullptr), m_H (nullptr), m_Htau (nullptr), + m_HPinv (nullptr) + { + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); + + if (nr <= 0 || nc <= 0) + (*current_liboctave_error_handler) + ("matrix dimension with negative or zero size"); + + if (order < 0 || order > 9) + (*current_liboctave_error_handler) + ("ordering %d is not supported by SPQR", order); + + CHOLMOD_NAME (start) (&m_cc); + const cholmod_sparse A = ros2rcs (a); + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + SuiteSparseQR (order, SPQR_DEFAULT_TOL, A.nrow, + const_cast (&A), &m_R, &m_E, &m_H, + &m_HPinv, &m_Htau, &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + +#elif defined (HAVE_CXSPARSE) , S (nullptr), N (nullptr) -#endif { -#if defined (HAVE_CXSPARSE) - CXSPARSE_DNAME () A; A.nzmax = a.nnz (); @@ -210,22 +317,36 @@ (*current_liboctave_error_handler) ("sparse_qr: sparse matrix QR factorization filled"); + } + #else + { octave_unused_parameter (order); (*current_liboctave_error_handler) - ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built"); + ("sparse_qr: support for SPQR or CXSparse was unavailable or disabled when liboctave was built"); + } #endif - } template <> sparse_qr::sparse_qr_rep::~sparse_qr_rep (void) { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + CHOLMOD_NAME (free_sparse) (&m_R, &m_cc); + CHOLMOD_NAME (free_sparse) (&m_H, &m_cc); + CHOLMOD_NAME (free_dense) (&m_Htau, &m_cc); + free (m_E); // FIXME: use cholmod_l_free + free (m_HPinv); + CHOLMOD_NAME (finish) (&m_cc); + +#elif defined (HAVE_CXSPARSE) + CXSPARSE_DNAME (_sfree) (S); CXSPARSE_DNAME (_nfree) (N); + #endif } @@ -233,7 +354,11 @@ SparseMatrix sparse_qr::sparse_qr_rep::V (void) const { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + return rcs2ros (m_H); + +#elif defined (HAVE_CXSPARSE) // Drop zeros from V and sort // FIXME: Is the double transpose to sort necessary? @@ -272,7 +397,31 @@ SparseMatrix sparse_qr::sparse_qr_rep::R (bool econ) const { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + octave_idx_type nr = static_cast (m_R->nrow); + octave_idx_type nc = static_cast (m_R->ncol); + octave_idx_type nz = static_cast (m_R->nzmax); + + // FIXME: Does this work if econ = true? + SparseMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz); + octave_idx_type *Rp = to_octave_idx_type_ptr + (static_cast (m_R->p)); + octave_idx_type *Ri = to_octave_idx_type_ptr + (static_cast (m_R->i)); + + for (octave_idx_type j = 0; j < nc + 1; j++) + ret.xcidx (j) = Rp[j]; + + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = Ri[j]; + ret.xdata (j) = (static_cast (m_R->x))[j]; + } + + return ret; + +#elif defined (HAVE_CXSPARSE) // Drop zeros from R and sort // FIXME: Is the double transpose to sort necessary? @@ -312,9 +461,40 @@ template <> Matrix - sparse_qr::sparse_qr_rep::C (const Matrix& b) const + sparse_qr::sparse_qr_rep::C (const Matrix &b) { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + Matrix ret (b_nr, b_nc); + + if (nrows != b_nr) + (*current_liboctave_error_handler) + ("sparse_qr: matrix dimension mismatch"); + else if (b_nc <= 0 || b_nr <= 0) + (*current_liboctave_error_handler) + ("sparse_qr: matrix dimension with negative or zero size"); + + cholmod_dense *QTB; // Q' * B + const cholmod_dense B = rod2rcd (b); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + QTB = SuiteSparseQR_qmult (SPQR_QTX, m_H, m_Htau, m_HPinv, + const_cast(&B), &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + // copy QTB into ret + double *QTB_x = static_cast (QTB->x); + double *ret_vec = static_cast (ret.fortran_vec ()); + for (octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < b_nr; i++) + ret_vec[j * b_nr + i] = QTB_x[j * b_nr + i]; + + CHOLMOD_NAME (free_dense) (&QTB, &m_cc); + + return ret; + +#elif defined (HAVE_CXSPARSE) octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); @@ -378,13 +558,93 @@ template <> Matrix - sparse_qr::sparse_qr_rep::Q (void) const + sparse_qr::sparse_qr_rep::C (const Matrix &b, bool econ) { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + octave_idx_type nr = (econ + ? (ncols > nrows ? nrows : ncols) + : nrows); + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + Matrix ret (nr, b_nc); + + if (nrows != b_nr) + (*current_liboctave_error_handler) + ("sparse_qr: matrix dimension mismatch"); + else if (b_nc <= 0 || b_nr <= 0) + (*current_liboctave_error_handler) + ("sparse_qr: matrix dimension with negative or zero size"); + + cholmod_dense *QTB; // Q' * B + const cholmod_dense B = rod2rcd (b); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + QTB = SuiteSparseQR_qmult (SPQR_QTX, m_H, m_Htau, m_HPinv, + const_cast(&B), &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + // copy QTB into ret + double *QTB_x = static_cast (QTB->x); + double *ret_vec = static_cast (ret.fortran_vec ()); + for (octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + ret_vec[j * nr + i] = QTB_x[j * b_nr + i]; + + CHOLMOD_NAME (free_dense) (&QTB, &m_cc); + + return ret; + +#else + + octave_unused_parameter (b); + + return Matrix (); + +#endif + } + + template <> + Matrix + sparse_qr::sparse_qr_rep::Q (void) + { +#if defined (HAVE_SPQR) + + Matrix ret (nrows, nrows); + cholmod_dense *q; + + // I is nrows x nrows identity matrix + cholmod_dense *I = static_cast + (CHOLMOD_NAME (allocate_dense) + (nrows, nrows, nrows, CHOLMOD_REAL, &m_cc)); + + for (octave_idx_type i = 0; i < nrows * nrows; i++) + (static_cast (I->x))[i] = 0.0; + + for (octave_idx_type i = 0; i < nrows; i++) + (static_cast (I->x))[i * nrows + i] = 1.0; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + q = SuiteSparseQR_qmult (SPQR_QX, m_H, m_Htau, m_HPinv, I, &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + double *q_x = static_cast (q->x); + double *ret_vec = const_cast (ret.fortran_vec ()); + for (octave_idx_type j = 0; j < nrows; j++) + for (octave_idx_type i = 0; i < nrows; i++) + ret_vec[j * nrows + i] = q_x[j * nrows + i]; + + CHOLMOD_NAME (free_dense) (&q, &m_cc); + CHOLMOD_NAME (free_dense) (&I, &m_cc); + + return ret; + +#elif defined (HAVE_CXSPARSE) + octave_idx_type nc = N->L->n; octave_idx_type nr = nrows; Matrix ret (nr, nr); - double *vec = ret.fortran_vec (); + double *ret_vec = ret.fortran_vec (); if (nr < 0 || nc < 0) (*current_liboctave_error_handler) ("matrix dimension mismatch"); @@ -400,7 +660,7 @@ OCTAVE_LOCAL_BUFFER (double, buf, S->m2); - for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) + for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx += nr) { octave_quit (); @@ -424,7 +684,7 @@ } for (octave_idx_type i = 0; i < nr; i++) - vec[i+idx] = buf[i]; + ret_vec[i+idx] = buf[i]; bvec[j] = 0.0; } @@ -440,14 +700,116 @@ } template <> + Matrix + sparse_qr::sparse_qr_rep::Q (bool econ) + { +#if defined (HAVE_SPQR) + + octave_idx_type nc = (econ + ? (ncols > nrows ? nrows : ncols) + : nrows); + Matrix ret (nrows, nc); + cholmod_dense *q; + + // I is nrows x nrows identity matrix + cholmod_dense *I = static_cast + (CHOLMOD_NAME (allocate_dense) + (nrows, nrows, nrows, CHOLMOD_REAL, &m_cc)); + + for (octave_idx_type i = 0; i < nrows * nrows; i++) + (static_cast (I->x))[i] = 0.0; + + for (octave_idx_type i = 0; i < nrows; i++) + (static_cast (I->x))[i * nrows + i] = 1.0; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + q = SuiteSparseQR_qmult (SPQR_QX, m_H, m_Htau, m_HPinv, I, &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + double *q_x = static_cast (q->x); + double *ret_vec = const_cast (ret.fortran_vec ()); + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nrows; i++) + ret_vec[j * nrows + i] = q_x[j * nrows + i]; + + CHOLMOD_NAME (free_dense) (&q, &m_cc); + CHOLMOD_NAME (free_dense) (&I, &m_cc); + + return ret; + +#else + + return Matrix (); + +#endif + } + + template <> template <> Matrix sparse_qr::sparse_qr_rep::tall_solve, Matrix> - (const MArray& b, octave_idx_type& info) const + (const MArray& b, octave_idx_type& info) { info = -1; -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + Matrix x (ncols, b_nc); // X = m_E'*(m_R\(Q'*B)) + + if (nrows <= 0 || ncols <= 0 || b_nc <= 0 || b_nr <= 0) + (*current_liboctave_error_handler) + ("matrix dimension with negative or zero size"); + + if (nrows < 0 || ncols < 0 || nrows != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + cholmod_dense *QTB; // Q' * B + const cholmod_dense B = rod2rcd (b); + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + // FIXME: Process b column by column as in the CXSPARSE version below. + // This avoids a large dense matrix Q' * B in memory. + QTB = SuiteSparseQR_qmult + (SPQR_QTX, m_H, m_Htau, m_HPinv, const_cast (&B), + &m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + spqr_error_handler (&m_cc); + + // convert m_R into CXSPARSE matrix R2 + CXSPARSE_DNAME (_sparse) R2; + R2.n = ncols; + R2.m = ncols; + R2.nzmax = m_R->nzmax; + R2.x = static_cast (m_R->x); + R2.i = static_cast (m_R->i); + R2.p = static_cast (m_R->p); + R2.nz = -1; + double *x_vec = const_cast (x.fortran_vec ()); + for (volatile octave_idx_type j = 0; j < b_nc; j++) + { + // fill x(:,j) + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + // solve (m_R\(Q'*B(:,j)) and store result in QTB(:,j) + CXSPARSE_DNAME (_usolve) + (&R2, &(static_cast (QTB->x)[j * b_nr])); + // x(:,j) = m_E' * (m_R\(Q'*B(:,j)) + CXSPARSE_DNAME (_ipvec) + (m_E, &(static_cast (QTB->x)[j * b_nr]), &x_vec[j * ncols], + ncols); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + + CHOLMOD_NAME (free_dense) (&QTB, &m_cc); + + info = 0; + + return x; + +#elif defined (HAVE_CXSPARSE) octave_idx_type nr = nrows; octave_idx_type nc = ncols; @@ -509,7 +871,6 @@ (const MArray& b, octave_idx_type& info) const { info = -1; - #if defined (HAVE_CXSPARSE) // These are swapped because the original matrix was transposed in @@ -562,7 +923,6 @@ return x; #else - octave_unused_parameter (b); return Matrix (); @@ -574,7 +934,7 @@ template <> SparseMatrix sparse_qr::sparse_qr_rep::tall_solve - (const SparseMatrix& b, octave_idx_type& info) const + (const SparseMatrix& b, octave_idx_type& info) { info = -1; @@ -758,7 +1118,7 @@ template <> ComplexMatrix sparse_qr::sparse_qr_rep::tall_solve, ComplexMatrix> - (const MArray& b, octave_idx_type& info) const + (const MArray& b, octave_idx_type& info) { info = -1; @@ -949,12 +1309,35 @@ sparse_qr::sparse_qr_rep::sparse_qr_rep (const SparseComplexMatrix& a, int order) : count (1), nrows (a.rows ()), ncols (a.columns ()) -#if defined (HAVE_CXSPARSE) - , S (nullptr), N (nullptr) -#endif +#if defined (HAVE_SPQR) + , m_cc (), m_R (nullptr), m_E (nullptr), m_H (nullptr), + m_Htau (nullptr), m_HPinv (nullptr) { -#if defined (HAVE_CXSPARSE) - + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); + + if (nr <= 0 || nc <= 0) + (*current_liboctave_error_handler) + ("matrix dimension with negative or zero size"); + + if (order < 0 || order > 9) + (*current_liboctave_error_handler) + ("ordering %d is not supported by SPQR", order); + + CHOLMOD_NAME (start) (&m_cc); + const cholmod_sparse A = cos2ccs (a); + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + SuiteSparseQR (order, SPQR_DEFAULT_TOL, A.nrow, + const_cast(&A), &m_R, &m_E, &m_H, + &m_HPinv, &m_Htau, &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + +#elif defined (HAVE_CXSPARSE) + , S (nullptr), N (nullptr) + { CXSPARSE_ZNAME () A; A.nzmax = a.nnz (); @@ -979,22 +1362,36 @@ (*current_liboctave_error_handler) ("sparse_qr: sparse matrix QR factorization filled"); + } + #else + { octave_unused_parameter (order); (*current_liboctave_error_handler) ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built"); + } #endif - } template <> sparse_qr::sparse_qr_rep::~sparse_qr_rep (void) { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + CHOLMOD_NAME (free_sparse) (&m_R, &m_cc); + CHOLMOD_NAME (free_sparse) (&m_H, &m_cc); + CHOLMOD_NAME (free_dense) (&m_Htau, &m_cc); + free (m_E); // FIXME: use cholmod_l_free + free (m_HPinv); + CHOLMOD_NAME (finish) (&m_cc); + +#elif defined (HAVE_CXSPARSE) + CXSPARSE_ZNAME (_sfree) (S); CXSPARSE_ZNAME (_nfree) (N); + #endif } @@ -1040,7 +1437,32 @@ SparseComplexMatrix sparse_qr::sparse_qr_rep::R (bool econ) const { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + octave_idx_type nr = static_cast (m_R->nrow); + octave_idx_type nc = static_cast (m_R->ncol); + octave_idx_type nz = static_cast (m_R->nzmax); + + // FIXME: Does this work if econ = true? + SparseComplexMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz); + octave_idx_type *Rp = to_octave_idx_type_ptr + (static_cast (m_R->p)); + octave_idx_type *Ri = to_octave_idx_type_ptr + (static_cast (m_R->i)); + + for (octave_idx_type j = 0; j < nc + 1; j++) + ret.xcidx (j) = Rp[j]; + + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = Ri[j]; + ret.xdata (j) = (static_cast (m_R->x))[j]; + } + + return ret; + +#elif defined (HAVE_CXSPARSE) + // Drop zeros from R and sort // FIXME: Is the double transpose to sort necessary? @@ -1080,9 +1502,45 @@ template <> ComplexMatrix - sparse_qr::sparse_qr_rep::C (const ComplexMatrix& b) const + sparse_qr::sparse_qr_rep::C (const ComplexMatrix& b) { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + // FIXME: not tested + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + ComplexMatrix ret (b_nr, b_nc); + + if (nrows != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + if (b_nc <= 0 || b_nr <= 0) + (*current_liboctave_error_handler) + ("matrix dimension with negative or zero size"); + + cholmod_dense *QTB; // Q' * B + const cholmod_dense B = cod2ccd (b); + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + QTB = SuiteSparseQR_qmult (SPQR_QTX, m_H, m_Htau, m_HPinv, + const_cast (&B), + &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + // copy QTB into ret + Complex *QTB_x = static_cast (QTB->x); + Complex *ret_vec = static_cast (ret.fortran_vec ()); + for (octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < b_nr; i++) + ret_vec[j * b_nr + i] = QTB_x[j * b_nr + i]; + + CHOLMOD_NAME (free_dense) (&QTB, &m_cc); + + return ret; + +#elif defined (HAVE_CXSPARSE) + octave_idx_type b_nr = b.rows (); octave_idx_type b_nc = b.cols (); octave_idx_type nc = N->L->n; @@ -1143,9 +1601,97 @@ template <> ComplexMatrix - sparse_qr::sparse_qr_rep::Q (void) const + sparse_qr::sparse_qr_rep::C + (const ComplexMatrix& b, bool econ) { -#if defined (HAVE_CXSPARSE) +#if defined (HAVE_SPQR) + + // FIXME: not tested + octave_idx_type nr = (econ + ? (ncols > nrows ? nrows : ncols) + : nrows); + octave_idx_type b_nr = b.rows (); + octave_idx_type b_nc = b.cols (); + ComplexMatrix ret (nr, b_nc); + + if (nrows != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + + if (b_nc <= 0 || b_nr <= 0) + (*current_liboctave_error_handler) + ("matrix dimension with negative or zero size"); + + cholmod_dense *QTB; // Q' * B + const cholmod_dense B = cod2ccd (b); + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + QTB = SuiteSparseQR_qmult (SPQR_QTX, m_H, m_Htau, m_HPinv, + const_cast (&B), + &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + // copy QTB into ret + Complex *QTB_x = static_cast (QTB->x); + Complex *ret_vec = static_cast (ret.fortran_vec ()); + for (octave_idx_type j = 0; j < b_nc; j++) + for (octave_idx_type i = 0; i < nr; i++) + ret_vec[j * nr + i] = QTB_x[j * b_nr + i]; + + CHOLMOD_NAME (free_dense) (&QTB, &m_cc); + + return ret; + +#else + + octave_unused_parameter (b); + + return ComplexMatrix (); + +#endif + } + + + template <> + ComplexMatrix + sparse_qr::sparse_qr_rep::Q (void) + { +#if defined (HAVE_SPQR) + + ComplexMatrix ret (nrows, nrows); + cholmod_dense *q; + + // I is nrows x nrows identity matrix + cholmod_dense *I = static_cast + (CHOLMOD_NAME (allocate_dense) + (nrows, nrows, nrows, CHOLMOD_COMPLEX, &m_cc)); + + for (octave_idx_type i = 0; i < nrows * nrows; i++) + (static_cast (I->x))[i] = 0.0; + + for (octave_idx_type i = 0; i < nrows; i++) + (static_cast (I->x))[i * nrows + i] = 1.0; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + q = SuiteSparseQR_qmult (SPQR_QX, m_H, m_Htau, m_HPinv, I, + &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + Complex *q_x = static_cast (q->x); + Complex *ret_vec = const_cast (ret.fortran_vec ()); + + for (octave_idx_type j = 0; j < nrows; j++) + for (octave_idx_type i = 0; i < nrows; i++) + ret_vec[j * nrows + i] = q_x[j * nrows + i]; + + CHOLMOD_NAME (free_dense) (&q, &m_cc); + CHOLMOD_NAME (free_dense) (&I, &m_cc); + + return ret; + +#elif defined (HAVE_CXSPARSE) + octave_idx_type nc = N->L->n; octave_idx_type nr = nrows; ComplexMatrix ret (nr, nr); @@ -1206,11 +1752,59 @@ } template <> + ComplexMatrix + sparse_qr::sparse_qr_rep::Q (bool econ) + { +#if defined (HAVE_SPQR) + + octave_idx_type nc = (econ + ? (ncols > nrows ? nrows : ncols) + : nrows); + ComplexMatrix ret (nrows, nc); + cholmod_dense *q; + + // I is nrows x nrows identity matrix + cholmod_dense *I = static_cast + (CHOLMOD_NAME (allocate_dense) + (nrows, nrows, nrows, CHOLMOD_COMPLEX, &m_cc)); + + for (octave_idx_type i = 0; i < nrows * nrows; i++) + (static_cast (I->x))[i] = 0.0; + + for (octave_idx_type i = 0; i < nrows; i++) + (static_cast (I->x))[i * nrows + i] = 1.0; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + q = SuiteSparseQR_qmult (SPQR_QX, m_H, m_Htau, m_HPinv, I, + &m_cc); + spqr_error_handler (&m_cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + Complex *q_x = static_cast (q->x); + Complex *ret_vec = const_cast (ret.fortran_vec ()); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = 0; i < nrows; i++) + ret_vec[j * nrows + i] = q_x[j * nrows + i]; + + CHOLMOD_NAME (free_dense) (&q, &m_cc); + CHOLMOD_NAME (free_dense) (&I, &m_cc); + + return ret; + +#else + + return ComplexMatrix (); + +#endif + } + + template <> template <> SparseComplexMatrix sparse_qr::sparse_qr_rep::tall_solve - (const SparseComplexMatrix& b, octave_idx_type& info) const + (const SparseComplexMatrix& b, octave_idx_type& info) { info = -1; @@ -1448,7 +2042,7 @@ ComplexMatrix sparse_qr::sparse_qr_rep::tall_solve, ComplexMatrix> - (const MArray& b, octave_idx_type& info) const + (const MArray& b, octave_idx_type& info) { info = -1; @@ -1590,7 +2184,7 @@ SparseComplexMatrix sparse_qr::sparse_qr_rep::tall_solve - (const SparseMatrix& b, octave_idx_type& info) const + (const SparseMatrix& b, octave_idx_type& info) { info = -1; @@ -1790,7 +2384,7 @@ ComplexMatrix sparse_qr::sparse_qr_rep::tall_solve, ComplexMatrix> - (const MArray& b, octave_idx_type& info) const + (const MArray& b, octave_idx_type& info) { info = -1; @@ -1928,8 +2522,9 @@ template <> template <> SparseComplexMatrix - sparse_qr::sparse_qr_rep::tall_solve - (const SparseComplexMatrix& b, octave_idx_type& info) const + sparse_qr::sparse_qr_rep::tall_solve + (const SparseComplexMatrix& b, octave_idx_type& info) { info = -1; @@ -2023,7 +2618,8 @@ template <> template <> SparseComplexMatrix - sparse_qr::sparse_qr_rep::wide_solve + sparse_qr::sparse_qr_rep::wide_solve (const SparseComplexMatrix& b, octave_idx_type& info) const { info = -1; @@ -2188,6 +2784,27 @@ { return rep->P (); } + + template + ColumnVector + sparse_qr::E (void) const + { + return rep->E(); + } + + + template + SparseMatrix + sparse_qr::E_MAT (void) const + { + ColumnVector perm = rep-> E (); + octave_idx_type nrows = perm.rows (); + SparseMatrix ret (nrows,nrows,nrows); + for (octave_idx_type i = 0; i < nrows; i++) + ret(perm(i) - 1,i) = 1.0; + return ret; + } + template SPARSE_T @@ -2205,11 +2822,25 @@ template typename SPARSE_T::dense_matrix_type + sparse_qr::C (const typename SPARSE_T::dense_matrix_type& b, + bool econ) const + { + return rep->C (b,econ); + } + + template + typename SPARSE_T::dense_matrix_type sparse_qr::Q (void) const { return rep->Q (); } + template + typename SPARSE_T::dense_matrix_type + sparse_qr::Q (bool econ) const + { + return rep->Q (econ); + } // 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 @@ -2228,7 +2859,11 @@ cxsparse_defaults { public: +#if defined (HAVE_SPQR) + enum { order = SPQR_ORDERING_DEFAULT }; +#elif defined (HAVE_CXSPARSE) enum { order = 3 }; +#endif }; template <> @@ -2236,7 +2871,11 @@ cxsparse_defaults { public: +#if defined (HAVE_SPQR) + enum { order = SPQR_ORDERING_DEFAULT }; +#elif defined (HAVE_CXSPARSE) enum { order = 2 }; +#endif }; template @@ -2245,6 +2884,32 @@ sparse_qr::solve (const SPARSE_T& a, const RHS_T& b, octave_idx_type& info) { +#if defined (HAVE_SPQR) + + 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::order; + + if (nr <= 0 || nc <= 0 || b_nc <= 0 || b_nr <= 0) + (*current_liboctave_error_handler) + ("matrix dimension with negative or zero size"); + + if ( nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + + info = 0; + + return min2norm_solve (a, b, info, order); + +#elif defined (HAVE_CXSPARSE) + info = -1; octave_idx_type nr = a.rows (); @@ -2277,8 +2942,371 @@ return q.ok () ? q.wide_solve (b, info) : RET_T (); } + +#endif } +#if defined (HAVE_SPQR) + //explicit instantiations of static member function solve + template + OCTAVE_API Matrix + sparse_qr::solve, Matrix> + (const SparseMatrix& a, const MArray& b, octave_idx_type& info); + + template + OCTAVE_API SparseMatrix + sparse_qr::solve + (const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info); + + template + OCTAVE_API ComplexMatrix + sparse_qr::solve, ComplexMatrix> + (const SparseMatrix& a, const MArray& b, octave_idx_type& info); + + template + OCTAVE_API SparseComplexMatrix + sparse_qr::solve + (const SparseMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info); + + template + OCTAVE_API ComplexMatrix + sparse_qr::solve, ComplexMatrix> + (const SparseComplexMatrix& a, const MArray& b, + octave_idx_type& info); + + template + OCTAVE_API SparseComplexMatrix + sparse_qr::solve< + SparseComplexMatrix, SparseComplexMatrix> + (const SparseComplexMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info); + + template + OCTAVE_API ComplexMatrix + sparse_qr::solve, ComplexMatrix> + (const SparseComplexMatrix& a, const MArray& b, + octave_idx_type& info); + + template + OCTAVE_API SparseComplexMatrix + sparse_qr::solve + (const SparseComplexMatrix& a, const SparseMatrix& b, + octave_idx_type& info); + + //explicit instantiations of member function E_MAT + template + OCTAVE_API SparseMatrix + sparse_qr::E_MAT (void) const; + + template + OCTAVE_API SparseMatrix + sparse_qr::E_MAT (void) const; + + //specializations of function min2norm_solve + template <> + template <> + OCTAVE_API Matrix + sparse_qr::min2norm_solve, Matrix> + (const SparseMatrix& a, const MArray& b, + octave_idx_type& info, int order) + { + info = -1; + octave_idx_type b_nc = b.cols (); + octave_idx_type nc = a.cols (); + Matrix x (nc, b_nc); + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + const cholmod_sparse A = ros2rcs (a); + const cholmod_dense B = rod2rcd (b); + cholmod_dense *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, + const_cast (&A), + const_cast (&B), &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + double *vec = x.fortran_vec (); + for (volatile octave_idx_type i = 0; i < nc * b_nc; i++) + vec[i] = static_cast (X->x)[i]; + + info = 0; + CHOLMOD_NAME (finish) (&cc); + + return x; + + } + + template <> + template <> + OCTAVE_API SparseMatrix + sparse_qr::min2norm_solve + (const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info, + int order) + { + info = -1; + SparseMatrix x; + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + const cholmod_sparse A = ros2rcs(a); + cholmod_sparse B = ros2rcs(b); + cholmod_sparse *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, + const_cast(&A), &B, + &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + x = rcs2ros (X); + CHOLMOD_NAME (finish) (&cc); + info = 0; + + return x; + + } + + template <> + template <> + OCTAVE_API ComplexMatrix + sparse_qr::min2norm_solve, ComplexMatrix> + (const SparseMatrix& a, const MArray& b, + octave_idx_type& info, int order) + { + info = -1; + + octave_idx_type b_nc = b.cols (); + octave_idx_type nc = a.cols (); + + ComplexMatrix x (nc, b_nc); + + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + + cholmod_sparse *A = ros2ccs (a, &cc); + const cholmod_dense B = cod2ccd (b); + cholmod_dense *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, A, + const_cast (&B), + &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + Complex *vec = x.fortran_vec (); + for (volatile octave_idx_type i = 0; i < nc * b_nc; i++) + vec[i] = static_cast (X->x)[i]; + + CHOLMOD_NAME (free_sparse) (&A, &cc); + CHOLMOD_NAME (finish) (&cc); + + info = 0; + + return x; + + } + + template <> + template <> + OCTAVE_API SparseComplexMatrix + sparse_qr::min2norm_solve + (const SparseMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info, int order) + { + info = -1; + + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + + cholmod_sparse * A = ros2ccs (a, &cc); + const cholmod_sparse B = cos2ccs (b); + cholmod_sparse *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, A, + const_cast (&B), + &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + CHOLMOD_NAME (free_sparse) (&A, &cc); + CHOLMOD_NAME (finish) (&cc); + + SparseComplexMatrix ret = ccs2cos(X); + + info = 0; + + return ret; + + } + + template <> + template <> + OCTAVE_API ComplexMatrix + sparse_qr::min2norm_solve, + ComplexMatrix> + (const SparseComplexMatrix& a, const MArray& b, + octave_idx_type& info,int order) + { + info = -1; + octave_idx_type b_nc = b.cols (); + octave_idx_type nc = a.cols (); + ComplexMatrix x (nc, b_nc); + + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + + const cholmod_sparse A = cos2ccs (a); + const cholmod_dense B = cod2ccd (b); + cholmod_dense *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, + const_cast (&A), + const_cast (&B), + &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + Complex *vec = x.fortran_vec (); + for (volatile octave_idx_type i = 0; i < nc * b_nc; i++) + vec[i] = static_cast (X->x)[i]; + + CHOLMOD_NAME (finish) (&cc); + + info = 0; + + return x; + + } + + template <> + template <> + OCTAVE_API ComplexMatrix + sparse_qr::min2norm_solve, + ComplexMatrix> + (const SparseComplexMatrix& a, const MArray& b, + octave_idx_type& info, int order) + { + info = -1; + + octave_idx_type b_nc = b.cols (); + octave_idx_type nc = a.cols (); + ComplexMatrix x (nc, b_nc); + + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + + const cholmod_sparse A = cos2ccs (a); + cholmod_dense *B = rod2ccd (b, &cc); + cholmod_dense *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, + const_cast (&A), + B, &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + Complex *vec = x.fortran_vec (); + + for (volatile octave_idx_type i = 0; i < nc * b_nc; i++) + vec[i] = static_cast (X->x)[i]; + + CHOLMOD_NAME (free_dense) (&B, &cc); + CHOLMOD_NAME (finish) (&cc); + + info = 0; + + return x; + + } + + template <> + template <> + OCTAVE_API SparseComplexMatrix + sparse_qr::min2norm_solve + (const SparseComplexMatrix& a, const SparseComplexMatrix& b, + octave_idx_type& info, int order) + { + info = -1; + + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + + const cholmod_sparse A = cos2ccs (a); + const cholmod_sparse B = cos2ccs (b); + cholmod_sparse *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, + const_cast (&A), + const_cast (&B), + &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + CHOLMOD_NAME (finish) (&cc); + + info = 0; + + return ccs2cos (X); + + } + + template <> + template <> + OCTAVE_API SparseComplexMatrix + sparse_qr::min2norm_solve + (const SparseComplexMatrix& a, const SparseMatrix& b, + octave_idx_type& info,int order) + { + info = -1; + + cholmod_common cc; + + CHOLMOD_NAME (start) (&cc); + + const cholmod_sparse A = cos2ccs (a); + cholmod_sparse *B = ros2ccs (b, &cc); + cholmod_sparse *X; + + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL, + const_cast (&A), B, + &cc); + spqr_error_handler (&cc); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + CHOLMOD_NAME (finish) (&cc); + + SparseComplexMatrix ret = ccs2cos(X); + + CHOLMOD_NAME (free_sparse) (&B, &cc); + CHOLMOD_NAME (finish) (&cc); + + info = 0; + + return ret; + + } + +#endif + template template RET_T @@ -2304,6 +3332,7 @@ sparse_qr::sparse_qr (const sparse_qr& a); template OCTAVE_API sparse_qr::~sparse_qr (void); template OCTAVE_API bool sparse_qr::ok (void) const; + template OCTAVE_API ColumnVector sparse_qr::E (void) const; template OCTAVE_API SparseMatrix sparse_qr::V (void) const; template OCTAVE_API ColumnVector sparse_qr::Pinv (void) const; template OCTAVE_API ColumnVector sparse_qr::P (void) const; @@ -2311,7 +3340,10 @@ sparse_qr::R (bool econ) const; template OCTAVE_API Matrix sparse_qr::C (const Matrix& b) const; + template OCTAVE_API Matrix + sparse_qr::C (const Matrix& b, bool econ) const; template OCTAVE_API Matrix sparse_qr::Q (void) const; + template OCTAVE_API Matrix sparse_qr::Q (bool econ) const; template OCTAVE_API sparse_qr::sparse_qr (void); template OCTAVE_API @@ -2322,6 +3354,8 @@ (const sparse_qr& a); template OCTAVE_API sparse_qr::~sparse_qr (void); template OCTAVE_API bool sparse_qr::ok (void) const; + template OCTAVE_API ColumnVector + sparse_qr::E (void) const; template OCTAVE_API SparseComplexMatrix sparse_qr::V (void) const; template OCTAVE_API ColumnVector @@ -2333,7 +3367,11 @@ template OCTAVE_API ComplexMatrix sparse_qr::C (const ComplexMatrix& b) const; template OCTAVE_API ComplexMatrix + sparse_qr::C (const ComplexMatrix& b, bool econ) const; + template OCTAVE_API ComplexMatrix sparse_qr::Q (void) const; + template OCTAVE_API ComplexMatrix + sparse_qr::Q (bool econ) const; Matrix qrsolve (const SparseMatrix& a, const MArray& b, diff -r 4ab10bf27340 -r 385b07c4ee70 liboctave/numeric/sparse-qr.h --- a/liboctave/numeric/sparse-qr.h Mon Jan 04 10:23:28 2021 +0100 +++ b/liboctave/numeric/sparse-qr.h Thu Jan 09 21:45:54 2020 +0100 @@ -54,7 +54,12 @@ OCTAVE_API sparse_qr (void); +#if (HAVE_SPQR) + // order = 7 selects SPQR default ordering + OCTAVE_API sparse_qr (const SPARSE_T& a, int order = 7); +#else OCTAVE_API sparse_qr (const SPARSE_T& a, int order = 0); +#endif OCTAVE_API sparse_qr (const sparse_qr& a); @@ -64,6 +69,11 @@ OCTAVE_API bool ok (void) const; + OCTAVE_API ColumnVector E (void) const; + + // constructs permutation matrix from permutation vector rep -> E() + OCTAVE_API SparseMatrix E_MAT () const; + OCTAVE_API SPARSE_T V (void) const; OCTAVE_API ColumnVector Pinv (void) const; @@ -76,8 +86,14 @@ C (const typename SPARSE_T::dense_matrix_type& b) const; OCTAVE_API typename SPARSE_T::dense_matrix_type + C (const typename SPARSE_T::dense_matrix_type& b, bool econ) const; + + OCTAVE_API typename SPARSE_T::dense_matrix_type Q (void) const; + OCTAVE_API typename SPARSE_T::dense_matrix_type + Q (bool econ) const; + template static OCTAVE_API RET_T solve (const SPARSE_T& a, const RHS_T& b, @@ -85,6 +101,11 @@ private: + template + static OCTAVE_API RET_T + min2norm_solve (const SPARSE_T& a, const RHS_T& b, + octave_idx_type& info, int order); + class sparse_qr_rep; sparse_qr_rep *rep; diff -r 4ab10bf27340 -r 385b07c4ee70 liboctave/util/oct-sparse.cc --- a/liboctave/util/oct-sparse.cc Mon Jan 04 10:23:28 2021 +0100 +++ b/liboctave/util/oct-sparse.cc Thu Jan 09 21:45:54 2020 +0100 @@ -78,6 +78,202 @@ return reinterpret_cast (i); } + + const cholmod_sparse + ros2rcs (const SparseMatrix &a) + { + cholmod_sparse A; + + A.ncol = a.cols (); + A.nrow = a.rows (); +# if defined (OCTAVE_ENABLE_64) + A.itype = CHOLMOD_LONG; +# else + A.itype = CHOLMOD_INT; +# endif + A.nzmax = a.nnz (); + A.sorted = 0; + A.packed = 1; + A.stype = 0; + A.xtype = CHOLMOD_REAL; + A.dtype = CHOLMOD_DOUBLE; + A.nz = NULL; + A.z = NULL; + A.p = const_cast (to_suitesparse_intptr (a.cidx ())); + A.i = const_cast (to_suitesparse_intptr (a.ridx ())); + A.x = const_cast (a.data ()); + + return A; + } + + const cholmod_sparse + cos2ccs (const SparseComplexMatrix &a) + { + cholmod_sparse A; + + A.ncol = a.cols (); + A.nrow = a.rows (); +# if defined (OCTAVE_ENABLE_64) + A.itype = CHOLMOD_LONG; +# else + A.itype = CHOLMOD_INT; +# endif + A.nzmax = a.nnz (); + A.sorted = 0; + A.packed = 1; + A.stype = 0; + A.xtype = CHOLMOD_COMPLEX; + A.dtype = CHOLMOD_DOUBLE; + A.nz = NULL; + A.z = NULL; + A.p = const_cast (to_suitesparse_intptr (a.cidx ())); + A.i = const_cast (to_suitesparse_intptr (a.ridx ())); + A.x = const_cast + (reinterpret_cast (a.data ())); + + return A; + } + + cholmod_dense* + rod2ccd (const MArray &a, cholmod_common* cc1) + { + cholmod_dense *A = static_cast + (CHOLMOD_NAME (allocate_dense) + (a.rows (), a.cols (), a.rows(), CHOLMOD_COMPLEX, + cc1)); + + const double *a_x = a.data (); + + for (octave_idx_type j = 0; j < a.cols() * a.rows() ; j++) + (static_cast (A->x))[j] = Complex (a_x[j], 0.0); + + return A; + } + + const cholmod_dense + rod2rcd (const MArray &a) + { + cholmod_dense A; + + A.ncol = a.cols (); + A.nrow = a.rows (); + A.nzmax = a.cols() * a.rows(); + A.xtype = CHOLMOD_REAL; + A.dtype = CHOLMOD_DOUBLE; + A.z = NULL; + A.d = a.rows(); + A.x = const_cast (a.data ()); + + return A; + } + + const cholmod_dense + cod2ccd (const ComplexMatrix &a) + { + cholmod_dense A; + + A.ncol = a.cols (); + A.nrow = a.rows (); + A.nzmax = a.cols () * a.rows (); + A.xtype = CHOLMOD_COMPLEX; + A.dtype = CHOLMOD_DOUBLE; + A.z = NULL; + A.d = a.rows(); + A.x = const_cast + (reinterpret_cast (a.data ())); + + return A; + } + + SparseMatrix + rcs2ros (const cholmod_sparse* y) + { + SparseMatrix ret (static_cast (y->nrow), + static_cast (y->ncol), + static_cast (y->nzmax)); + + octave_idx_type nz = static_cast (y->nzmax); + + for (octave_idx_type j = 0; j < static_cast (y->ncol) + 1; + j++) + ret.xcidx (j) = (static_cast (y->p))[j]; + + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = (static_cast (y->i))[j]; + ret.xdata (j) = (static_cast (y->x))[j]; + } + + return ret; + } + + SparseComplexMatrix + ccs2cos (const cholmod_sparse* a) + { + SparseComplexMatrix ret (static_cast (a->nrow), + static_cast (a->ncol), + static_cast (a->nzmax)); + + octave_idx_type nz = static_cast (a->nzmax); + + for (octave_idx_type j = 0; j < static_cast (a->ncol) + 1; + j++) + ret.xcidx (j) = (static_cast (a->p))[j]; + + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = (static_cast (a->i))[j]; + ret.xdata (j) = (static_cast (a->x))[j]; + } + + return ret; + } + + cholmod_sparse* + ros2ccs (const SparseMatrix& a, cholmod_common* cc1) + { + cholmod_sparse *A = static_cast + (CHOLMOD_NAME(allocate_sparse) + (a.rows (), a.cols (), a.nnz (), 0, 1, 0, + CHOLMOD_COMPLEX, cc1)); + + for (octave_idx_type j = 0; + j < static_cast (a.cols ()) + 1; j++) + static_cast (A->p)[j] = a.cidx(j); + + const double *a_x = a.data (); + for (octave_idx_type j = 0; j < a.nnz (); j++) + { + (static_cast (A->x))[j] = Complex (a_x[j], 0.0); + (static_cast (A->i))[j] = a.ridx(j); + } + return A; + } + + void + spqr_error_handler (const cholmod_common* cc) + { + if (cc->status >= 0) + return; + + switch (cc->status) + { + case CHOLMOD_OUT_OF_MEMORY: + (*current_liboctave_error_handler) + ("sparse_qr: sparse matrix QR factorization failed" + " - out of memory"); + case CHOLMOD_TOO_LARGE: + (*current_liboctave_error_handler) + ("sparse_qr: sparse matrix QR factorization failed" + " - integer overflow occurred"); + default: + (*current_liboctave_error_handler) + ("sparse_qr: sparse matrix QR factorization failed"); + } + + // FIXME: Free memory? + // FIXME: Can cc-status > 0 (CHOLMOD_NOT_POSDEF, CHOLMOD_DSMALL) occur? + } } #endif diff -r 4ab10bf27340 -r 385b07c4ee70 liboctave/util/oct-sparse.h --- a/liboctave/util/oct-sparse.h Mon Jan 04 10:23:28 2021 +0100 +++ b/liboctave/util/oct-sparse.h Thu Jan 09 21:45:54 2020 +0100 @@ -28,6 +28,9 @@ #include "octave-config.h" +#include "dSparse.h" +#include "CSparse.h" + #if defined (HAVE_SUITESPARSE_AMD_H) # include #elif defined (HAVE_UFSPARSE_AMD_H) @@ -88,6 +91,10 @@ # include #endif +#if defined (HAVE_SUITESPARSE_SUITESPARSEQR_HPP) +# include +#endif + // Cope with new SuiteSparse versions #if defined (SUITESPARSE_VERSION) @@ -182,6 +189,51 @@ extern OCTAVE_API const octave_idx_type* to_octave_idx_type_ptr (const suitesparse_integer *i); + + // Convert real sparse octave matrix to real sparse cholmod matrix. + // Returns a "shallow" copy of a. + extern const OCTAVE_API cholmod_sparse + ros2rcs (const SparseMatrix& a); + + // Convert real sparse cholmod matrix to real sparse octave matrix. + // Returns a "shallow" copy of y. + extern OCTAVE_API SparseMatrix + rcs2ros (const cholmod_sparse *y); + + // Convert real dense octave matrix to real dense cholmod matrix. + // Returns a "shallow" copy of a. + extern const OCTAVE_API cholmod_dense + rod2rcd (const MArray& a); + + // Convert complex dense octave matrix to complex dense cholmod matrix. + // Returns a "shallow" copy of a. + extern const OCTAVE_API cholmod_dense + cod2ccd (const ComplexMatrix &a); + + // Convert complex sparse octave matrix to complex sparse cholmod matrix. + // Returns a "shallow" copy of a. + extern const OCTAVE_API cholmod_sparse + cos2ccs (const SparseComplexMatrix &a); + + // Convert complex sparse cholmod matrix to complex sparse octave matrix. + // Returns a "deep" copy of a. + extern OCTAVE_API SparseComplexMatrix + ccs2cos (const cholmod_sparse *a); + + // Convert real sparse octave matrix to complex sparse cholmod matrix. + // Returns a "deep" copy of a. + // FIXME: const return type not necessary since deep copy is returned. + extern OCTAVE_API cholmod_sparse * + ros2ccs (const SparseMatrix& a, cholmod_common *cc1); + + // Convert real dense octave matrix to complex dense cholmod matrix. + // Returns a "deep" copy of a. + // FIXME: const return type not necessary since deep copy is returned. + extern OCTAVE_API cholmod_dense * + rod2ccd (const MArray &a, cholmod_common *cc1); + + extern OCTAVE_API void + spqr_error_handler (const cholmod_common *cc); } #endif