# HG changeset patch
# User simon.hau79@gmail.com
# Date 1578602754 -3600
# Thu Jan 09 21:45:54 2020 +0100
# Node ID a9ea9c972e0323a616708ba1f7677af05a2ba030
# Parent bd51beb6205e9d4160f6a36dbd791bcd3750169e
Implementation of qr(...,0), qr(...,"vector") and qr(...,"matrix") for sparse A.
Economy versions of C, Q and R are returned if qr(...,0) is called.
Column permutation of A is returned as (vector|matrix) if qr(...,("vector"|"matrix"))
is called.
Additional tests if HAVE_SPQR is set.
diff -r bd51beb6205e -r a9ea9c972e03 configure.ac
--- a/configure.ac Wed Jan 08 11:59:41 2020 -0500
+++ b/configure.ac Thu Jan 09 21:45:54 2020 +0100
@@ -1,11 +1,10 @@
dnl Process this file with autoconf to produce a configure script.
dnl
-########################################################################
-###
### Copyright (C) 1993-2020 The Octave Project Developers
###
-### See the file COPYRIGHT.md in the top-level directory of this
-### distribution or .
+### See the file COPYRIGHT.md in the top-level directory of this distribution
+### or .
+###
###
### This file is part of Octave.
###
@@ -22,8 +21,6 @@ dnl
### You should have received a copy of the GNU General Public License
### along with Octave; see the file COPYING. If not, see
### .
-###
-########################################################################
### Initialize Autoconf
AC_PREREQ([2.65])
@@ -2108,6 +2105,16 @@ OCTAVE_CHECK_LIB(cholmod, CHOLMOD,
[], [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,
@@ -2173,11 +2180,11 @@ dnl substitute some of them into other s
## 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)
@@ -3105,6 +3112,9 @@ Octave is now configured for $canonical_
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 bd51beb6205e -r a9ea9c972e03 libinterp/dldfcn/qr.cc
--- a/libinterp/dldfcn/qr.cc Wed Jan 08 11:59:41 2020 -0500
+++ b/libinterp/dldfcn/qr.cc Thu Jan 09 21:45:54 2020 +0100
@@ -1,27 +1,28 @@
-////////////////////////////////////////////////////////////////////////
-//
-// Copyright (C) 1996-2020 The Octave Project Developers
-//
-// See the file COPYRIGHT.md in the top-level directory of this
-// distribution or .
-//
-// 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
-// .
-//
-////////////////////////////////////////////////////////////////////////
+/*
+
+Copyright (C) 1996-2020 The Octave Project Developers
+
+See the file COPYRIGHT.md in the top-level directory of this distribution
+or .
+
+
+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
+.
+
+*/
#if defined (HAVE_CONFIG_H)
# include "config.h"
@@ -331,45 +332,160 @@ orthogonal basis of @code{span (A)}.
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
@@ -896,7 +1012,7 @@ orthogonal basis of @code{span (A)}.
%! assert (r\c, full (a)\b, 10e-10);
## Test under-determined systems!!
-%!#testif HAVE_CXSPARSE
+%!testif HAVE_CXSPARSE
%! n = 20; d = 0.2;
%! ## initialize generators to make behavior reproducible
%! rand ("state", 42);
@@ -946,7 +1062,7 @@ orthogonal basis of @code{span (A)}.
%! assert (r\c, full (a)\b, 10e-10);
## Test under-determined systems!!
-%!#testif HAVE_CXSPARSE
+%!testif HAVE_CXSPARSE
%! n = 20; d = 0.2;
%! ## initialize generators to make behavior reproducible
%! rand ("state", 42);
@@ -956,6 +1072,146 @@ orthogonal basis of @code{span (A)}.
%! [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 bd51beb6205e -r a9ea9c972e03 liboctave/numeric/sparse-qr.cc
--- a/liboctave/numeric/sparse-qr.cc Wed Jan 08 11:59:41 2020 -0500
+++ b/liboctave/numeric/sparse-qr.cc Thu Jan 09 21:45:54 2020 +0100
@@ -1,27 +1,28 @@
-////////////////////////////////////////////////////////////////////////
-//
-// Copyright (C) 2005-2020 The Octave Project Developers
-//
-// See the file COPYRIGHT.md in the top-level directory of this
-// distribution or .
-//
-// 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
-// .
-//
-////////////////////////////////////////////////////////////////////////
+/*
+
+Copyright (C) 2005-2020 The Octave Project Developers
+
+See the file COPYRIGHT.md in the top-level directory of this distribution
+or .
+
+
+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
+.
+
+*/
#if defined (HAVE_CONFIG_H)
# include "config.h"
@@ -47,15 +48,18 @@ namespace octave
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 +73,10 @@ namespace octave
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 +102,9 @@ namespace octave
bool ok (void) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+ return (H && Htau && HPinv && R_ && E_ && cc);
+#elif defined (HAVE_CXSPARSE)
return (N && S);
#else
return false;
@@ -108,21 +117,42 @@ namespace octave
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;
typename SPARSE_T::dense_matrix_type
+ C (const typename SPARSE_T::dense_matrix_type& b, bool econ) const;
+
+ typename SPARSE_T::dense_matrix_type
Q (void) const;
+ typename SPARSE_T::dense_matrix_type
+ Q (bool econ) const;
+
refcount count;
octave_idx_type nrows;
octave_idx_type ncols;
+ // octave_idx_type econ;
+
+#if defined (HAVE_SPQR)
typename cxsparse_types::symbolic_type *S;
typename cxsparse_types::numeric_type *N;
+ template
+ RET_T
+ 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
@@ -131,6 +161,21 @@ namespace octave
template
RET_T
wide_solve (const RHS_T& b, octave_idx_type& info) const;
+
+#if defined (HAVE_SPQR)
+
+ private:
+
+ cholmod_common *cc;
+ cholmod_sparse *R_; // R factor
+ // Column permutation for A. Fill-reducing ordering.
+ suitesparse_integer *E_;
+ cholmod_sparse *H; // Householder vectors
+ suitesparse_integer *HPinv;
+ cholmod_dense *Htau; // beta scalars
+
+#endif
+
};
template
@@ -157,7 +202,16 @@ namespace octave
ColumnVector
sparse_qr::sparse_qr_rep::P (void) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+ ColumnVector ret (nrows);
+
+ // FIXME: Is ret.xelem (HPinv[i]) = i + 1 correct?
+ for (octave_idx_type i = 0; i < nrows; i++)
+ ret.xelem (HPinv[i]) = i + 1;
+
+ return ret;
+
+#elif defined (HAVE_CXSPARSE)
ColumnVector ret (N->L->m);
@@ -173,20 +227,77 @@ namespace octave
#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 (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)
+ , R_ (nullptr), E_ (nullptr), H (nullptr), HPinv (nullptr), Htau (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 is not supported by SPQR");
+
+ cholmod_common Common;
+ cc = &Common;
+ CHOLMOD_NAME (start) (cc);
+ const cholmod_sparse A = ros2rcs (a);
+
+ BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+ SuiteSparseQR (order, SPQR_DEFAULT_TOL, A.nrow,
+ const_cast (&A), &R_, &E_, &H,
+ &HPinv, &Htau, cc);
+ spqr_error_handler (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 +321,36 @@ namespace octave
(*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) (&R_, cc);
+ CHOLMOD_NAME (free_sparse) (&H, cc);
+ CHOLMOD_NAME (free_dense) (&Htau, cc);
+ free (E_); // FIXME: use cholmod_l_free
+ free (HPinv);
+ CHOLMOD_NAME (finish) (cc);
+
+#elif defined (HAVE_CXSPARSE)
+
CXSPARSE_DNAME (_sfree) (S);
CXSPARSE_DNAME (_nfree) (N);
+
#endif
}
@@ -233,7 +358,11 @@ namespace octave
SparseMatrix
sparse_qr::sparse_qr_rep::V (void) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+ return rcs2ros (H);
+
+#elif defined (HAVE_CXSPARSE)
// Drop zeros from V and sort
// FIXME: Is the double transpose to sort necessary?
@@ -272,7 +401,31 @@ namespace octave
SparseMatrix
sparse_qr::sparse_qr_rep::R (bool econ) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+ octave_idx_type nr = static_cast (R_->nrow);
+ octave_idx_type nc = static_cast (R_->ncol);
+ octave_idx_type nz = static_cast (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 (R_->p));
+ octave_idx_type *Ri = to_octave_idx_type_ptr
+ (static_cast (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 (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 +465,45 @@ namespace octave
template <>
Matrix
- sparse_qr::sparse_qr_rep::C (const Matrix& b) const
+ sparse_qr::sparse_qr_rep::C (const Matrix &b) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+ octave_idx_type b_nr = b.rows ();
+ octave_idx_type b_nc = b.cols ();
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME(start) (cc1);
+ 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, H, Htau,HPinv,
+ const_cast(&B), cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+
+ CHOLMOD_NAME(finish) (cc1);
+
+ return ret;
+
+#elif defined (HAVE_CXSPARSE)
octave_idx_type b_nr = b.rows ();
octave_idx_type b_nc = b.cols ();
@@ -378,13 +567,156 @@ namespace octave
template <>
Matrix
+ sparse_qr::sparse_qr_rep::C (const Matrix &b, bool econ) const
+ {
+#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 ();
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME(start) (cc1);
+ 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, H, Htau,HPinv,
+ const_cast(&B), cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+
+ CHOLMOD_NAME(finish) (cc1);
+
+ return ret;
+
+#else
+
+ octave_unused_parameter (b);
+
+ return Matrix ();
+
+#endif
+ }
+
+ template <>
+ Matrix
+ sparse_qr::sparse_qr_rep::Q (bool econ) const
+ {
+#if defined (HAVE_SPQR)
+
+ octave_idx_type nc = (econ
+ ? (ncols > nrows ? nrows : ncols)
+ : nrows);
+ Matrix ret (nrows, nc);
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ cholmod_dense *q;
+
+ // I is nrows x nrows identity matrix
+ cholmod_dense *I = static_cast
+ (CHOLMOD_NAME (allocate_dense)
+ (nrows, nrows, nrows, CHOLMOD_REAL, cc1));
+
+ 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, H, Htau, HPinv, I, cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+ CHOLMOD_NAME (free_dense) (&I, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ return ret;
+
+#else
+
+ return Matrix ();
+
+#endif
+ }
+
+ template <>
+ Matrix
sparse_qr::sparse_qr_rep::Q (void) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+ Matrix ret (nrows, nrows);
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ cholmod_dense *q;
+
+ // I is nrows x nrows identity matrix
+ cholmod_dense *I = static_cast
+ (CHOLMOD_NAME (allocate_dense)
+ (nrows, nrows, nrows, CHOLMOD_REAL, cc1));
+
+ 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, H, Htau, HPinv, I, cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+ CHOLMOD_NAME (free_dense) (&I, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ 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");
@@ -424,7 +756,7 @@ namespace octave
}
for (octave_idx_type i = 0; i < nr; i++)
- vec[i+idx] = buf[i];
+ ret_vec[i+idx] = buf[i];
bvec[j] = 0.0;
}
@@ -447,7 +779,66 @@ namespace octave
{
info = -1;
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+ octave_idx_type b_nr = b.rows ();
+ octave_idx_type b_nc = b.cols ();
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME (start) (cc1);
+ Matrix x (ncols, b_nc); // X = E_'*(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, H, Htau, HPinv, const_cast (&B), cc1);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+ spqr_error_handler (cc1);
+
+ // convert R_ into CXSPARSE matrix R2
+ CXSPARSE_DNAME (_sparse) R2;
+ R2.n = ncols;
+ R2.m = ncols;
+ R2.nzmax = R_->nzmax;
+ R2.x = static_cast (R_->x);
+ R2.i = static_cast (R_->i);
+ R2.p = static_cast (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 (R_\(Q'*B(:,j)) and store result in QTB(:,j)
+ CXSPARSE_DNAME (_usolve)
+ (&R2, &(static_cast (QTB->x)[j * b_nr]));
+ // x(:,j) = E_' * (R_\(Q'*B(:,j))
+ CXSPARSE_DNAME (_ipvec)
+ (E_, &(static_cast (QTB->x)[j * b_nr]), &x_vec[j * ncols],
+ ncols);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+ }
+
+ CHOLMOD_NAME (free_dense) (&QTB, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ info = 0;
+
+ return x;
+
+#elif defined (HAVE_CXSPARSE)
octave_idx_type nr = nrows;
octave_idx_type nc = ncols;
@@ -509,7 +900,6 @@ namespace octave
(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 +952,6 @@ namespace octave
return x;
#else
-
octave_unused_parameter (b);
return Matrix ();
@@ -949,12 +1338,38 @@ namespace octave
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)
+ , R_ (nullptr), E_ (nullptr), H (nullptr), HPinv (nullptr),
+ Htau (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 is not supported by SPQR");
+
+ cholmod_common Common;
+ cc = &Common;
+
+ CHOLMOD_NAME (start) (cc);
+ const cholmod_sparse A = cos2ccs (a);
+
+ BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+ SuiteSparseQR (order, SPQR_DEFAULT_TOL, A.nrow,
+ const_cast(&A), &R_, &E_, &H,
+ &HPinv, &Htau, cc);
+ spqr_error_handler (cc);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+ }
+
+#elif defined (HAVE_CXSPARSE)
+ , S (nullptr), N (nullptr)
+ {
CXSPARSE_ZNAME () A;
A.nzmax = a.nnz ();
@@ -979,22 +1394,36 @@ namespace octave
(*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) (&R_, cc);
+ CHOLMOD_NAME (free_sparse) (&H, cc);
+ CHOLMOD_NAME (free_dense) (&Htau, cc);
+ free (E_); // FIXME: use cholmod_l_free
+ free (HPinv);
+ CHOLMOD_NAME (finish) (cc);
+
+#elif defined (HAVE_CXSPARSE)
+
CXSPARSE_ZNAME (_sfree) (S);
CXSPARSE_ZNAME (_nfree) (N);
+
#endif
}
@@ -1040,7 +1469,32 @@ namespace octave
SparseComplexMatrix
sparse_qr::sparse_qr_rep::R (bool econ) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+ octave_idx_type nr = static_cast (R_->nrow);
+ octave_idx_type nc = static_cast (R_->ncol);
+ octave_idx_type nz = static_cast (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 (R_->p));
+ octave_idx_type *Ri = to_octave_idx_type_ptr
+ (static_cast (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 (R_->x))[j];
+ }
+
+ return ret;
+
+#elif defined (HAVE_CXSPARSE)
+
// Drop zeros from R and sort
// FIXME: Is the double transpose to sort necessary?
@@ -1082,7 +1536,48 @@ namespace octave
ComplexMatrix
sparse_qr::sparse_qr_rep::C (const ComplexMatrix& b) const
{
-#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 ();
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME (start) (cc1);
+ 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, H, Htau, HPinv,
+ const_cast (&B),
+ cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+
+ CHOLMOD_NAME (finish) (cc1);
+
+ 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 +1638,104 @@ namespace octave
template <>
ComplexMatrix
+ sparse_qr::sparse_qr_rep::C
+ (const ComplexMatrix& b,bool econ) const
+ {
+#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 ();
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME (start) (cc1);
+ 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, H, Htau, HPinv,
+ const_cast (&B),
+ cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+
+ CHOLMOD_NAME (finish) (cc1);
+
+ return ret;
+
+#else
+
+ octave_unused_parameter (b);
+
+ return ComplexMatrix ();
+
+#endif
+ }
+
+
+ template <>
+ ComplexMatrix
sparse_qr::sparse_qr_rep::Q (void) const
{
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+ ComplexMatrix ret (nrows, nrows);
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME (start) (cc1);
+ cholmod_dense *q;
+
+ // I is nrows x nrows identity matrix
+ cholmod_dense *I = static_cast
+ (CHOLMOD_NAME (allocate_dense)
+ (nrows, nrows, nrows, CHOLMOD_COMPLEX, cc1));
+
+ 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, H, Htau, HPinv, I, cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+ CHOLMOD_NAME (free_dense) (&I, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ return ret;
+
+#elif defined (HAVE_CXSPARSE)
octave_idx_type nc = N->L->n;
octave_idx_type nr = nrows;
ComplexMatrix ret (nr, nr);
@@ -1206,6 +1796,57 @@ namespace octave
}
template <>
+ ComplexMatrix
+ sparse_qr::sparse_qr_rep::Q (bool econ) const
+ {
+#if defined (HAVE_SPQR)
+
+ octave_idx_type nc = (econ
+ ? (ncols > nrows ? nrows : ncols)
+ : nrows);
+ ComplexMatrix ret (nrows, nc);
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+ CHOLMOD_NAME (start) (cc1);
+ cholmod_dense *q;
+
+ // I is nrows x nrows identity matrix
+ cholmod_dense *I = static_cast
+ (CHOLMOD_NAME (allocate_dense)
+ (nrows, nrows, nrows, CHOLMOD_COMPLEX, cc1));
+
+ 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, H, Htau, HPinv, I, cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+ CHOLMOD_NAME (free_dense) (&I, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ return ret;
+
+#else
+
+ return ComplexMatrix ();
+
+#endif
+ }
+
+ template <>
template <>
SparseComplexMatrix
sparse_qr::sparse_qr_rep::tall_solveP ();
}
+
+ 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 +2867,25 @@ namespace octave
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 +2904,11 @@ namespace octave
cxsparse_defaults
{
public:
+#if defined (HAVE_SPQR)
+ enum { order = SPQR_ORDERING_DEFAULT };
+#elif defined (HAVE_CXSPARSE)
enum { order = 3 };
+#endif
};
template <>
@@ -2236,7 +2916,11 @@ namespace octave
cxsparse_defaults
{
public:
+#if defined (HAVE_SPQR)
+ enum { order = SPQR_ORDERING_DEFAULT };
+#elif defined (HAVE_CXSPARSE)
enum { order = 2 };
+#endif
};
template
@@ -2245,6 +2929,33 @@ namespace octave
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 octave::math::sparse_qr::min2norm_solve
+ (a, b, info, order);
+
+#elif defined (HAVE_CXSPARSE)
+
info = -1;
octave_idx_type nr = a.rows ();
@@ -2277,8 +2988,382 @@ namespace octave
return q.ok () ? q.wide_solve (b, info) : RET_T ();
}
+
+#endif
}
+#if defined (HAVE_SPQR)
+ //explicit instantiations of static member function solve
+ template
+ Matrix
+ sparse_qr::solve, Matrix>
+ (const SparseMatrix& a, const MArray& b, octave_idx_type& info);
+
+ template
+ SparseMatrix
+ sparse_qr::solve
+ (const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info);
+
+ template
+ ComplexMatrix
+ sparse_qr::solve>, ComplexMatrix>
+ (const SparseMatrix& a, const MArray>& b,
+ octave_idx_type& info);
+
+ template
+ SparseComplexMatrix
+ sparse_qr::solve
+ (const SparseMatrix& a, const SparseComplexMatrix& b,
+ octave_idx_type& info);
+
+ template
+ ComplexMatrix
+ sparse_qr::solve<
+ MArray>, ComplexMatrix>
+ (const SparseComplexMatrix& a, const MArray>& b,
+ octave_idx_type& info);
+
+ template
+ SparseComplexMatrix
+ sparse_qr::solve<
+ SparseComplexMatrix, SparseComplexMatrix>
+ (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+ octave_idx_type& info);
+
+ template
+ ComplexMatrix
+ sparse_qr::solve, ComplexMatrix>
+ (const SparseComplexMatrix& a, const MArray& b,
+ octave_idx_type& info);
+
+ template
+ SparseComplexMatrix
+ sparse_qr::solve
+ (const SparseComplexMatrix& a, const SparseMatrix& b,
+ octave_idx_type& info);
+
+ //explicit instantiations of member function E_MAT
+ template
+ SparseMatrix
+ sparse_qr::E_MAT (void) const;
+
+ template
+ SparseMatrix
+ sparse_qr::E_MAT (void) const;
+
+ //specializations of function min2norm_solve
+ template <>
+ template <>
+ 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 Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+ 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), cc1);
+ spqr_error_handler (cc1);
+ 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) (cc1);
+
+ return x;
+
+ }
+
+ template <>
+ template <>
+ SparseMatrix
+ sparse_qr::min2norm_solve
+ (const SparseMatrix& a, const SparseMatrix& b,
+ octave_idx_type& info, int order)
+ {
+ info = -1;
+ SparseMatrix x;
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+ 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,
+ cc1);
+ spqr_error_handler (cc1);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+ x = rcs2ros (X);
+ CHOLMOD_NAME (finish) (cc1);
+ info = 0;
+
+ return x;
+
+ }
+
+ template <>
+ template <>
+ 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 Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ cholmod_sparse *A = ros2ccs (a, cc1);
+ 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),
+ cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ info = 0;
+
+ return x;
+
+ }
+
+ template <>
+ template <>
+ SparseComplexMatrix
+ sparse_qr::min2norm_solve
+ (const SparseMatrix& a, const SparseComplexMatrix& b,
+ octave_idx_type& info, int order)
+ {
+ info = -1;
+
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ cholmod_sparse * A = ros2ccs (a, cc1);
+ 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),
+ cc1);
+ spqr_error_handler (cc1);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+ CHOLMOD_NAME (free_sparse) (&A, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ SparseComplexMatrix ret = ccs2cos(X);
+
+ info = 0;
+
+ return ret;
+
+ }
+
+ template <>
+ template <>
+ 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 Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ 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),
+ cc1);
+ spqr_error_handler (cc1);
+ 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) (cc1);
+
+ info = 0;
+
+ return x;
+
+ }
+
+ template <>
+ template <>
+ 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 Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ const cholmod_sparse A = cos2ccs (a);
+ cholmod_dense *B = rod2ccd (b, cc1);
+ cholmod_dense *X;
+
+ BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+ X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL,
+ const_cast (&A),
+ B, cc1);
+ spqr_error_handler (cc1);
+ 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, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ info = 0;
+
+ return x;
+
+ }
+
+ template <>
+ template <>
+ SparseComplexMatrix
+ sparse_qr::min2norm_solve
+ (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+ octave_idx_type& info, int order)
+ {
+ info = -1;
+
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ 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),
+ cc1);
+ spqr_error_handler (cc1);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+ CHOLMOD_NAME (finish) (cc1);
+
+ info = 0;
+
+ return ccs2cos (X);
+
+ }
+
+ template <>
+ template <>
+ SparseComplexMatrix
+ sparse_qr::min2norm_solve
+ (const SparseComplexMatrix& a, const SparseMatrix& b,
+ octave_idx_type& info,int order)
+ {
+ info = -1;
+
+ cholmod_common Common1;
+ cholmod_common *cc1 = &Common1;
+
+ CHOLMOD_NAME (start) (cc1);
+
+ const cholmod_sparse A = cos2ccs (a);
+ cholmod_sparse *B = ros2ccs (b,cc1);
+ cholmod_sparse *X;
+
+ BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+ X = SuiteSparseQR_min2norm (order, SPQR_DEFAULT_TOL,
+ const_cast (&A), B,
+ cc1);
+ spqr_error_handler (cc1);
+ END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+ CHOLMOD_NAME (finish) (cc1);
+
+ SparseComplexMatrix ret = ccs2cos(X);
+
+ CHOLMOD_NAME (free_sparse) (&B, cc1);
+ CHOLMOD_NAME (finish) (cc1);
+
+ info = 0;
+
+ return ret;
+
+ }
+
+#endif
+
template
template
RET_T
diff -r bd51beb6205e -r a9ea9c972e03 liboctave/numeric/sparse-qr.h
--- a/liboctave/numeric/sparse-qr.h Wed Jan 08 11:59:41 2020 -0500
+++ b/liboctave/numeric/sparse-qr.h Thu Jan 09 21:45:54 2020 +0100
@@ -1,27 +1,28 @@
-////////////////////////////////////////////////////////////////////////
-//
-// Copyright (C) 2005-2020 The Octave Project Developers
-//
-// See the file COPYRIGHT.md in the top-level directory of this
-// distribution or .
-//
-// 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
-// .
-//
-////////////////////////////////////////////////////////////////////////
+/*
+
+Copyright (C) 2005-2020 The Octave Project Developers
+
+See the file COPYRIGHT.md in the top-level directory of this distribution
+or .
+
+
+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
+.
+
+*/
#if ! defined (octave_sparse_qr_h)
#define octave_sparse_qr_h 1
@@ -54,7 +55,12 @@ namespace octave
sparse_qr (void);
+#if (HAVE_SPQR)
+ // order = 7 selects SPQR default ordering
+ sparse_qr (const SPARSE_T& a, int order = 7);
+#else
sparse_qr (const SPARSE_T& a, int order = 0);
+#endif
sparse_qr (const sparse_qr& a);
@@ -64,6 +70,12 @@ namespace octave
bool ok (void) const;
+ ColumnVector E (void) const;
+
+ //constructs permutation matrix from
+ //permutation vector rep -> E()
+ SparseMatrix E_MAT () const;
+
SPARSE_T V (void) const;
ColumnVector Pinv (void) const;
@@ -76,8 +88,15 @@ namespace octave
C (const typename SPARSE_T::dense_matrix_type& b) const;
typename SPARSE_T::dense_matrix_type
+ C (const typename SPARSE_T::dense_matrix_type& b,
+ bool econ) const;
+
+ typename SPARSE_T::dense_matrix_type
Q (void) const;
+ typename SPARSE_T::dense_matrix_type
+ Q (bool econ) const;
+
template
static RET_T
solve (const SPARSE_T& a, const RHS_T& b,
@@ -85,6 +104,11 @@ namespace octave
private:
+ template
+ static 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 bd51beb6205e -r a9ea9c972e03 liboctave/util/oct-sparse.cc
--- a/liboctave/util/oct-sparse.cc Wed Jan 08 11:59:41 2020 -0500
+++ b/liboctave/util/oct-sparse.cc Thu Jan 09 21:45:54 2020 +0100
@@ -1,27 +1,28 @@
-////////////////////////////////////////////////////////////////////////
-//
-// Copyright (C) 2017-2020 The Octave Project Developers
-//
-// See the file COPYRIGHT.md in the top-level directory of this
-// distribution or .
-//
-// 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
-// .
-//
-////////////////////////////////////////////////////////////////////////
+/*
+
+Copyright (C) 2017-2020 The Octave Project Developers
+
+See the file COPYRIGHT.md in the top-level directory of this distribution
+or .
+
+
+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
+.
+
+*/
#if defined (HAVE_CONFIG_H)
# include "config.h"
@@ -29,6 +30,8 @@
#include "lo-error.h"
#include "oct-sparse.h"
+#include "dSparse.h"
+#include "CSparse.h"
#if (defined (HAVE_AMD) || defined (HAVE_CCOLAMD) \
|| defined (HAVE_CHOLMOD) || defined (HAVE_COLAMD) \
@@ -78,6 +81,204 @@ namespace octave
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 bd51beb6205e -r a9ea9c972e03 liboctave/util/oct-sparse.h
--- a/liboctave/util/oct-sparse.h Wed Jan 08 11:59:41 2020 -0500
+++ b/liboctave/util/oct-sparse.h Thu Jan 09 21:45:54 2020 +0100
@@ -1,31 +1,34 @@
-////////////////////////////////////////////////////////////////////////
-//
-// Copyright (C) 2005-2020 The Octave Project Developers
-//
-// See the file COPYRIGHT.md in the top-level directory of this
-// distribution or .
-//
-// 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
-// .
-//
-////////////////////////////////////////////////////////////////////////
+/*
+
+Copyright (C) 2005-2020 The Octave Project Developers
+
+See the file COPYRIGHT.md in the top-level directory of this distribution
+or .
+
+
+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
+.
+
+*/
#if ! defined (octave_oct_sparse_h)
#define octave_oct_sparse_h 1
+#include "dSparse.h"
+#include "CSparse.h"
#include "octave-config.h"
#if defined (HAVE_SUITESPARSE_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 @@ namespace octave
extern 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 cholmod_sparse
+ ros2rcs (const SparseMatrix& a);
+
+ // Convert real sparse cholmod matrix to real sparse octave matrix.
+ // Returns a "shallow" copy of y.
+ extern SparseMatrix
+ rcs2ros (const cholmod_sparse *y);
+
+ // Convert real dense octave matrix to real dense cholmod matrix.
+ // Returns a "shallow" copy of a.
+ extern const cholmod_dense
+ rod2rcd (const MArray& a);
+
+ // Convert complex dense octave matrix to complex dense cholmod matrix.
+ // Returns a "shallow" copy of a.
+ extern const cholmod_dense
+ cod2ccd (const ComplexMatrix &a);
+
+ // Convert complex sparse octave matrix to complex sparse cholmod matrix.
+ // Returns a "shallow" copy of a.
+ extern const cholmod_sparse
+ cos2ccs (const SparseComplexMatrix &a);
+
+ // Convert complex sparse cholmod matrix to complex sparse octave matrix.
+ // Returns a "deep" copy of a.
+ extern 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 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 cholmod_dense *
+ rod2ccd (const MArray