Mercurial > octave-libtiff

changeset 29259:66f162b6fa03

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 test for SPQR.
author Simon Hau <simon.hau79@gmail.com>
date Thu, 09 Jan 2020 21:45:54 +0100
parents fdfd874293f6
children 18ad3c01fc2a
files configure.ac libinterp/corefcn/qr.cc liboctave/numeric/sparse-qr.cc liboctave/numeric/sparse-qr.h liboctave/util/oct-sparse.cc liboctave/util/oct-sparse.h
diffstat 6 files changed, 1715 insertions(+), 129 deletions(-) [+]
line wrap: on
line diff
--- a/configure.ac	Mon Jan 04 14:52:06 2021 +0100
+++ b/configure.ac	Thu Jan 09 21:45:54 2020 +0100
@@ -2144,6 +2144,14 @@
   [], [don't use CHOLMOD library, disable some sparse matrix functionality])
 LIBS="$save_LIBS"
 
+### Check for SPQR library
+
+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 +2217,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 +3156,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
--- a/libinterp/corefcn/qr.cc	Mon Jan 04 14:52:06 2021 +0100
+++ b/libinterp/corefcn/qr.cc	Thu Jan 09 21:45:54 2020 +0100
@@ -79,6 +79,8 @@
     return octave::math::qr<T>::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{CXSparse} (e.g., 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<SparseComplexMatrix> 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<SparseComplexMatrix>::solve
+                     <MArray<Complex>, 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<SparseComplexMatrix>::solve
+                     <SparseComplexMatrix, SparseComplexMatrix>
+                     (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<SparseComplexMatrix>::solve
+                     <MArray<double>, 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<SparseComplexMatrix>::solve
+                     <SparseMatrix, SparseComplexMatrix>
+                     (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<SparseComplexMatrix>
+              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<SparseComplexMatrix>
+              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<SparseComplexMatrix>
+                  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<SparseComplexMatrix>
+                  q (arg.sparse_complex_matrix_value (), 0);
+                  retval = ovl (q.Q (economy), q.R (economy));
+                }
+              else
+                {
+                  octave::math::sparse_qr<SparseComplexMatrix>
+                  q (arg.sparse_complex_matrix_value (), 0);
+                  retval = ovl (q.R (economy));
+                }
+            }
         }
       else
         {
-          octave::math::sparse_qr<SparseMatrix> 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<SparseMatrix>::solve
+                                <SparseMatrix, SparseMatrix>
+                                (arg.sparse_matrix_value (),
+                                 args (1).sparse_matrix_value (), info));
+              else if (! args(1).issparse () && args(1).iscomplex ())
+                retval = ovl (octave::math::sparse_qr<SparseMatrix>::solve
+                                <MArray<Complex>, 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<SparseMatrix>::solve
+                                <MArray<double>, Matrix>
+                                (arg.sparse_matrix_value (),
+                                 args (1).matrix_value (), info));
+              else if (args(1).issparse () &&  args(1).iscomplex ())
+                retval = ovl (octave::math::sparse_qr<SparseMatrix>::solve
+                                <SparseComplexMatrix, SparseComplexMatrix>
+                                (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<SparseMatrix>
+              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<SparseMatrix>
+              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<SparseMatrix>
+                  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<SparseMatrix>
+                  q (arg.sparse_matrix_value (), 0);
+                  retval = ovl (q.Q (economy), q.R (economy));
+                }
+              else
+                {
+                  octave::math::sparse_qr<SparseMatrix>
+                  q (arg.sparse_matrix_value (), 0);
+                  retval = ovl (q.R (economy));
+                }
+            }
         }
     }
   else
@@ -443,8 +584,8 @@
                     octave::math::qr<FloatComplexMatrix> 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<ComplexMatrix> 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;
 
@@ -856,7 +997,7 @@
 ## The deactivated tests below can't be tested till rectangular back-subs is
 ## implemented for sparse matrices.
 
-%!testif HAVE_CXSPARSE
+%!testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -876,7 +1017,7 @@
 %! r = qr (a);
 %! assert (r'*r, a'*a, 1e-10);
 
-%!testif HAVE_CXSPARSE
+%!testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -885,7 +1026,7 @@
 %! [c,r] = qr (a, ones (n,1));
 %! assert (r\c, full (a)\ones (n,1), 10e-10);
 
-%!testif HAVE_CXSPARSE
+%!testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -896,7 +1037,7 @@
 %! assert (r\c, full (a)\b, 10e-10);
 
 ## Test under-determined systems!!
-%!#testif HAVE_CXSPARSE
+%!#testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -906,7 +1047,7 @@
 %! [c,r] = qr (a, b);
 %! assert (r\c, full (a)\b, 10e-10);
 
-%!testif HAVE_CXSPARSE
+%!testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -926,7 +1067,7 @@
 %! r = qr (a);
 %! assert (r'*r, a'*a, 1e-10);
 
-%!testif HAVE_CXSPARSE
+%!testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -935,7 +1076,7 @@
 %! [c,r] = qr (a, ones (n,1));
 %! assert (r\c, full (a)\ones (n,1), 10e-10);
 
-%!testif HAVE_CXSPARSE
+%!testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! rand ("state", 42);
@@ -946,16 +1087,156 @@
 %! assert (r\c, full (a)\b, 10e-10);
 
 ## Test under-determined systems!!
-%!#testif HAVE_CXSPARSE
+%!#testif ; __have_feature__ ("SPQR") || __have_feature__ ("CXSPARSE")
 %! n = 20;  d = 0.2;
 %! ## initialize generators to make behavior reproducible
 %! 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
--- a/liboctave/numeric/sparse-qr.cc	Mon Jan 04 14:52:06 2021 +0100
+++ b/liboctave/numeric/sparse-qr.cc	Thu Jan 09 21:45:54 2020 +0100
@@ -44,24 +44,19 @@
 {
   namespace math
   {
+#if defined (HAVE_CXSPARSE)
     template <typename SPARSE_T>
     class
     cxsparse_types
-    {
-    };
+      { };
 
     template <>
     class
     cxsparse_types<SparseMatrix>
     {
     public:
-#if defined (HAVE_CXSPARSE)
       typedef CXSPARSE_DNAME (s) symbolic_type;
       typedef CXSPARSE_DNAME (n) numeric_type;
-#else
-      typedef void symbolic_type;
-      typedef void numeric_type;
-#endif
     };
 
     template <>
@@ -69,14 +64,10 @@
     cxsparse_types<SparseComplexMatrix>
     {
     public:
-#if defined (HAVE_CXSPARSE)
       typedef CXSPARSE_ZNAME (s) symbolic_type;
       typedef CXSPARSE_ZNAME (n) numeric_type;
-#else
-      typedef void symbolic_type;
-      typedef void numeric_type;
+    };
 #endif
-    };
 
     template <typename SPARSE_T>
     class sparse_qr<SPARSE_T>::sparse_qr_rep
@@ -95,7 +86,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 +101,56 @@
 
       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<octave_idx_type> count;
 
       octave_idx_type nrows;
       octave_idx_type ncols;
 
+#if defined (HAVE_SPQR)
+
+      template <typename RHS_T, typename RET_T>
+      RET_T solve (const RHS_T& b, octave_idx_type& info) const;
+
+#elif defined (HAVE_CXSPARSE)
+
       typename cxsparse_types<SPARSE_T>::symbolic_type *S;
       typename cxsparse_types<SPARSE_T>::numeric_type *N;
 
+#endif
+
+      template <typename RHS_T, typename RET_T>
+      RET_T tall_solve (const RHS_T& b, octave_idx_type& info);
+
       template <typename RHS_T, typename RET_T>
-      RET_T
-      tall_solve (const RHS_T& b, octave_idx_type& info) const;
-
-      template <typename RHS_T, typename RET_T>
-      RET_T
-      wide_solve (const RHS_T& b, octave_idx_type& info) const;
+      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 <typename SPARSE_T>
@@ -157,7 +177,17 @@
     ColumnVector
     sparse_qr<SPARSE_T>::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 +203,76 @@
 #endif
     }
 
+    template <typename SPARSE_T>
+    ColumnVector
+    sparse_qr<SPARSE_T>::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<octave_idx_type> (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<SparseMatrix>::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<double> (order, SPQR_DEFAULT_TOL, A.nrow,
+                             const_cast<cholmod_sparse*> (&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 +296,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<SparseMatrix>::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 +333,11 @@
     SparseMatrix
     sparse_qr<SparseMatrix>::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 +376,31 @@
     SparseMatrix
     sparse_qr<SparseMatrix>::sparse_qr_rep::R (bool econ) const
     {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+      octave_idx_type nr = static_cast<octave_idx_type> (m_R->nrow);
+      octave_idx_type nc = static_cast<octave_idx_type> (m_R->ncol);
+      octave_idx_type nz = static_cast<octave_idx_type> (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<suitesparse_integer*> (m_R->p));
+      octave_idx_type *Ri = to_octave_idx_type_ptr
+                              (static_cast<suitesparse_integer*> (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<double*> (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 +440,40 @@
 
     template <>
     Matrix
-    sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b) const
+    sparse_qr<SparseMatrix>::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<double> (SPQR_QTX, m_H, m_Htau, m_HPinv,
+                                         const_cast<cholmod_dense*>(&B), &m_cc);
+      spqr_error_handler (&m_cc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      // copy QTB into ret
+      double *QTB_x = static_cast<double*> (QTB->x);
+      double *ret_vec = static_cast<double*> (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 +537,93 @@
 
     template <>
     Matrix
-    sparse_qr<SparseMatrix>::sparse_qr_rep::Q (void) const
+    sparse_qr<SparseMatrix>::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<double> (SPQR_QTX, m_H, m_Htau, m_HPinv,
+                                         const_cast<cholmod_dense*>(&B), &m_cc);
+      spqr_error_handler (&m_cc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      // copy QTB into ret
+      double *QTB_x = static_cast<double*> (QTB->x);
+      double *ret_vec = static_cast<double*> (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<SparseMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                             (nrows, nrows, nrows, CHOLMOD_REAL, &m_cc));
+
+      for (octave_idx_type i = 0; i < nrows * nrows; i++)
+        (static_cast<double*> (I->x))[i] = 0.0;
+
+      for (octave_idx_type i = 0; i < nrows; i++)
+       (static_cast<double*> (I->x))[i * nrows + i] = 1.0;
+
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      q = SuiteSparseQR_qmult<double> (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<double*> (q->x);
+      double *ret_vec = const_cast<double*> (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 +639,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 +663,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 +679,116 @@
     }
 
     template <>
+    Matrix
+    sparse_qr<SparseMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                             (nrows, nrows, nrows, CHOLMOD_REAL, &m_cc));
+
+      for (octave_idx_type i = 0; i < nrows * nrows; i++)
+        (static_cast<double*> (I->x))[i] = 0.0;
+
+      for (octave_idx_type i = 0; i < nrows; i++)
+       (static_cast<double*> (I->x))[i * nrows + i] = 1.0;
+
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      q = SuiteSparseQR_qmult<double> (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<double*> (q->x);
+      double *ret_vec = const_cast<double*> (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<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix>
-    (const MArray<double>& b, octave_idx_type& info) const
+    (const MArray<double>& b, octave_idx_type& info)
     {
       info = -1;
 
-#if defined (HAVE_CXSPARSE)
+#if (defined (HAVE_SPQR) && defined (HAVE_CXSPARSE))
+
+      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<double>
+             (SPQR_QTX, m_H, m_Htau, m_HPinv, const_cast<cholmod_dense*> (&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<double*> (m_R->x);
+      R2.i = static_cast<suitesparse_integer*> (m_R->i);
+      R2.p = static_cast<suitesparse_integer*> (m_R->p);
+      R2.nz = -1;
+      double *x_vec = const_cast<double*> (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<double*> (QTB->x)[j * b_nr]));
+          // x(:,j) = m_E' * (m_R\(Q'*B(:,j))
+          CXSPARSE_DNAME (_ipvec)
+            (m_E, &(static_cast<double*> (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 +850,6 @@
     (const MArray<double>& b, octave_idx_type& info) const
     {
       info = -1;
-
 #if defined (HAVE_CXSPARSE)
 
       // These are swapped because the original matrix was transposed in
@@ -562,7 +902,6 @@
       return x;
 
 #else
-
       octave_unused_parameter (b);
 
       return Matrix ();
@@ -574,7 +913,7 @@
     template <>
     SparseMatrix
     sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix>
-    (const SparseMatrix& b, octave_idx_type& info) const
+    (const SparseMatrix& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -758,7 +1097,7 @@
     template <>
     ComplexMatrix
     sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix>
-    (const MArray<Complex>& b, octave_idx_type& info) const
+    (const MArray<Complex>& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -949,12 +1288,35 @@
     sparse_qr<SparseComplexMatrix>::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<Complex> (order, SPQR_DEFAULT_TOL, A.nrow,
+                              const_cast<cholmod_sparse*>(&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 +1341,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<SparseComplexMatrix>::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 +1416,32 @@
     SparseComplexMatrix
     sparse_qr<SparseComplexMatrix>::sparse_qr_rep::R (bool econ) const
     {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+      octave_idx_type nr = static_cast<octave_idx_type> (m_R->nrow);
+      octave_idx_type nc = static_cast<octave_idx_type> (m_R->ncol);
+      octave_idx_type nz = static_cast<octave_idx_type> (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<suitesparse_integer*> (m_R->p));
+      octave_idx_type *Ri = to_octave_idx_type_ptr
+                              (static_cast<suitesparse_integer*> (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<Complex*> (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 +1481,45 @@
 
     template <>
     ComplexMatrix
-    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::C (const ComplexMatrix& b) const
+    sparse_qr<SparseComplexMatrix>::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<Complex> (SPQR_QTX, m_H, m_Htau, m_HPinv,
+                                          const_cast<cholmod_dense*> (&B),
+                                          &m_cc);
+      spqr_error_handler (&m_cc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      // copy QTB into ret
+      Complex *QTB_x = static_cast<Complex*> (QTB->x);
+      Complex *ret_vec = static_cast<Complex*> (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 +1580,97 @@
 
     template <>
     ComplexMatrix
-    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::Q (void) const
+    sparse_qr<SparseComplexMatrix>::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<Complex> (SPQR_QTX, m_H, m_Htau, m_HPinv,
+                                          const_cast<cholmod_dense*> (&B),
+                                          &m_cc);
+      spqr_error_handler (&m_cc);
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+
+      // copy QTB into ret
+      Complex *QTB_x = static_cast<Complex*> (QTB->x);
+      Complex *ret_vec = static_cast<Complex*> (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<SparseComplexMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                              (nrows, nrows, nrows, CHOLMOD_COMPLEX, &m_cc));
+
+      for (octave_idx_type i = 0; i < nrows * nrows; i++)
+        (static_cast<Complex*> (I->x))[i] = 0.0;
+
+      for (octave_idx_type i = 0; i < nrows; i++)
+        (static_cast<Complex*> (I->x))[i * nrows + i] = 1.0;
+
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      q = SuiteSparseQR_qmult<Complex> (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<Complex*> (q->x);
+      Complex *ret_vec = const_cast<Complex*> (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 +1731,59 @@
     }
 
     template <>
+    ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                              (nrows, nrows, nrows, CHOLMOD_COMPLEX, &m_cc));
+
+      for (octave_idx_type i = 0; i < nrows * nrows; i++)
+        (static_cast<Complex*> (I->x))[i] = 0.0;
+
+      for (octave_idx_type i = 0; i < nrows; i++)
+        (static_cast<Complex*> (I->x))[i * nrows + i] = 1.0;
+
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      q = SuiteSparseQR_qmult<Complex> (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<Complex*> (q->x);
+      Complex *ret_vec = const_cast<Complex*> (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<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
                                                        SparseComplexMatrix>
-      (const SparseComplexMatrix& b, octave_idx_type& info) const
+      (const SparseComplexMatrix& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -1448,7 +2021,7 @@
     ComplexMatrix
     sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<double>,
                                                               ComplexMatrix>
-      (const MArray<double>& b, octave_idx_type& info) const
+      (const MArray<double>& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -1590,7 +2163,7 @@
     SparseComplexMatrix
     sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseMatrix,
                                                               SparseComplexMatrix>
-      (const SparseMatrix& b, octave_idx_type& info) const
+      (const SparseMatrix& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -1790,7 +2363,7 @@
     ComplexMatrix
     sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>,
                                                               ComplexMatrix>
-      (const MArray<Complex>& b, octave_idx_type& info) const
+      (const MArray<Complex>& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -1928,8 +2501,9 @@
     template <>
     template <>
     SparseComplexMatrix
-    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix, SparseComplexMatrix>
-      (const SparseComplexMatrix& b, octave_idx_type& info) const
+    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
+                                                              SparseComplexMatrix>
+      (const SparseComplexMatrix& b, octave_idx_type& info)
     {
       info = -1;
 
@@ -2023,7 +2597,8 @@
     template <>
     template <>
     SparseComplexMatrix
-    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix, SparseComplexMatrix>
+    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::wide_solve<SparseComplexMatrix,
+                                                              SparseComplexMatrix>
       (const SparseComplexMatrix& b, octave_idx_type& info) const
     {
       info = -1;
@@ -2188,6 +2763,27 @@
     {
       return rep->P ();
     }
+    
+    template <typename SPARSE_T>
+    ColumnVector
+    sparse_qr<SPARSE_T>::E (void) const
+    {
+      return rep->E();
+    }
+
+
+    template <typename SPARSE_T>
+    SparseMatrix
+    sparse_qr<SPARSE_T>::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 <typename SPARSE_T>
     SPARSE_T
@@ -2205,11 +2801,25 @@
 
     template <typename SPARSE_T>
     typename SPARSE_T::dense_matrix_type
+    sparse_qr<SPARSE_T>::C (const typename SPARSE_T::dense_matrix_type& b,
+                            bool econ) const
+    {
+      return rep->C (b,econ);
+    }
+
+    template <typename SPARSE_T>
+    typename SPARSE_T::dense_matrix_type
     sparse_qr<SPARSE_T>::Q (void) const
     {
       return rep->Q ();
     }
 
+    template <typename SPARSE_T>
+    typename SPARSE_T::dense_matrix_type
+    sparse_qr<SPARSE_T>::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 +2838,11 @@
     cxsparse_defaults<SparseMatrix>
     {
     public:
+#if defined (HAVE_SPQR)
+      enum { order = SPQR_ORDERING_DEFAULT };
+#elif defined (HAVE_CXSPARSE)
       enum { order = 3 };
+#endif
     };
 
     template <>
@@ -2236,7 +2850,11 @@
     cxsparse_defaults<SparseComplexMatrix>
     {
     public:
+#if defined (HAVE_SPQR)
+      enum { order = SPQR_ORDERING_DEFAULT };
+#elif defined (HAVE_CXSPARSE)
       enum { order = 2 };
+#endif
     };
 
     template <typename SPARSE_T>
@@ -2245,6 +2863,32 @@
     sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b,
                                 octave_idx_type& info)
     {
+#if (defined (HAVE_SPQR) && defined (HAVE_CHOLMOD))
+
+      info = -1;
+
+      octave_idx_type nr = a.rows ();
+      octave_idx_type nc = a.cols ();
+
+      octave_idx_type b_nc = b.cols ();
+      octave_idx_type b_nr = b.rows ();
+
+      int order = cxsparse_defaults<SPARSE_T>::order;
+
+      if (nr <= 0 || nc <= 0 || 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<RHS_T, RET_T> (a, b, info, order);
+
+#elif defined (HAVE_CXSPARSE)
+
       info = -1;
 
       octave_idx_type nr = a.rows ();
@@ -2277,8 +2921,372 @@
 
           return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T ();
         }
+
+#endif
     }
 
+#if defined (HAVE_SPQR)
+    //explicit instantiations of static member function solve
+    template
+    OCTAVE_API Matrix
+    sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix>
+    (const SparseMatrix& a, const MArray<double>& b, octave_idx_type& info);
+
+    template
+    OCTAVE_API SparseMatrix
+    sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix>
+    (const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info);
+
+    template
+    OCTAVE_API ComplexMatrix
+    sparse_qr<SparseMatrix>::solve<MArray<Complex>, ComplexMatrix>
+    (const SparseMatrix& a, const MArray<Complex>& b, octave_idx_type& info);
+
+    template
+    OCTAVE_API SparseComplexMatrix
+    sparse_qr<SparseMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix>
+    (const SparseMatrix& a, const SparseComplexMatrix& b,
+     octave_idx_type& info);
+
+    template
+    OCTAVE_API ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<MArray<Complex>, ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<Complex>& b,
+     octave_idx_type& info);
+
+    template
+    OCTAVE_API SparseComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<
+    SparseComplexMatrix, SparseComplexMatrix>
+    (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+     octave_idx_type& info);
+
+    template
+    OCTAVE_API ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<MArray<double>, ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<double>& b,
+     octave_idx_type& info);
+
+    template
+    OCTAVE_API SparseComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<SparseMatrix, SparseComplexMatrix>
+    (const SparseComplexMatrix& a, const SparseMatrix& b,
+     octave_idx_type& info);
+
+    //explicit instantiations of member function E_MAT
+    template
+    OCTAVE_API SparseMatrix
+    sparse_qr<SparseMatrix>::E_MAT (void) const;
+
+    template
+    OCTAVE_API SparseMatrix
+    sparse_qr<SparseComplexMatrix>::E_MAT (void) const;
+
+#  if defined (HAVE_CHOLMOD)
+    //specializations of function min2norm_solve
+    template <>
+    template <>
+    OCTAVE_API Matrix
+    sparse_qr<SparseMatrix>::min2norm_solve<MArray<double>, Matrix>
+    (const SparseMatrix& a, const MArray<double>& 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<double> (order, SPQR_DEFAULT_TOL,
+                                          const_cast<cholmod_sparse*> (&A),
+                                          const_cast<cholmod_dense*> (&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<double*> (X->x)[i];
+
+      info = 0;
+      CHOLMOD_NAME (finish) (&cc);
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    OCTAVE_API SparseMatrix
+    sparse_qr<SparseMatrix>::min2norm_solve<SparseMatrix, SparseMatrix>
+    (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<double> (order, SPQR_DEFAULT_TOL,
+                                          const_cast<cholmod_sparse*>(&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<SparseMatrix>::min2norm_solve<MArray<Complex>, ComplexMatrix>
+    (const SparseMatrix& a, const MArray<Complex>& 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<Complex> (order, SPQR_DEFAULT_TOL, A,
+                                           const_cast<cholmod_dense*> (&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<Complex*> (X->x)[i];
+
+      CHOLMOD_NAME (free_sparse) (&A, &cc);
+      CHOLMOD_NAME (finish) (&cc);
+
+      info = 0;
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    OCTAVE_API SparseComplexMatrix
+    sparse_qr<SparseMatrix>::min2norm_solve<SparseComplexMatrix,
+                                            SparseComplexMatrix>
+    (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<Complex> (order, SPQR_DEFAULT_TOL, A,
+                                           const_cast<cholmod_sparse*> (&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<SparseComplexMatrix>::min2norm_solve<MArray<Complex>,
+                                                   ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<Complex>& 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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&A),
+                                           const_cast<cholmod_dense*> (&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<Complex*> (X->x)[i];
+
+      CHOLMOD_NAME (finish) (&cc);
+
+      info = 0;
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    OCTAVE_API ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::min2norm_solve<MArray<double>,
+                                                   ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<double>& 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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&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<Complex*> (X->x)[i];
+
+      CHOLMOD_NAME (free_dense) (&B, &cc);
+      CHOLMOD_NAME (finish) (&cc);
+
+      info = 0;
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    OCTAVE_API SparseComplexMatrix
+    sparse_qr<SparseComplexMatrix>::min2norm_solve<SparseComplexMatrix,
+                                                   SparseComplexMatrix>
+    (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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&A),
+                                           const_cast<cholmod_sparse*> (&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<SparseComplexMatrix>::min2norm_solve<SparseMatrix,
+                                                   SparseComplexMatrix>
+    (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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&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
+#endif
+
     template <typename SPARSE_T>
     template <typename RHS_T, typename RET_T>
     RET_T
@@ -2304,6 +3312,7 @@
     sparse_qr<SparseMatrix>::sparse_qr (const sparse_qr<SparseMatrix>& a);
     template OCTAVE_API sparse_qr<SparseMatrix>::~sparse_qr (void);
     template OCTAVE_API bool sparse_qr<SparseMatrix>::ok (void) const;
+    template OCTAVE_API ColumnVector sparse_qr<SparseMatrix>::E (void) const;
     template OCTAVE_API SparseMatrix sparse_qr<SparseMatrix>::V (void) const;
     template OCTAVE_API ColumnVector sparse_qr<SparseMatrix>::Pinv (void) const;
     template OCTAVE_API ColumnVector sparse_qr<SparseMatrix>::P (void) const;
@@ -2311,7 +3320,10 @@
     sparse_qr<SparseMatrix>::R (bool econ) const;
     template OCTAVE_API Matrix
     sparse_qr<SparseMatrix>::C (const Matrix& b) const;
+    template OCTAVE_API Matrix
+    sparse_qr<SparseMatrix>::C (const Matrix& b, bool econ) const;
     template OCTAVE_API Matrix sparse_qr<SparseMatrix>::Q (void) const;
+    template OCTAVE_API Matrix sparse_qr<SparseMatrix>::Q (bool econ) const;
 
     template OCTAVE_API sparse_qr<SparseComplexMatrix>::sparse_qr (void);
     template OCTAVE_API
@@ -2322,6 +3334,8 @@
     (const sparse_qr<SparseComplexMatrix>& a);
     template OCTAVE_API sparse_qr<SparseComplexMatrix>::~sparse_qr (void);
     template OCTAVE_API bool sparse_qr<SparseComplexMatrix>::ok (void) const;
+    template OCTAVE_API ColumnVector
+    sparse_qr<SparseComplexMatrix>::E (void) const;
     template OCTAVE_API SparseComplexMatrix
     sparse_qr<SparseComplexMatrix>::V (void) const;
     template OCTAVE_API ColumnVector
@@ -2333,7 +3347,11 @@
     template OCTAVE_API ComplexMatrix
     sparse_qr<SparseComplexMatrix>::C (const ComplexMatrix& b) const;
     template OCTAVE_API ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::C (const ComplexMatrix& b, bool econ) const;
+    template OCTAVE_API ComplexMatrix
     sparse_qr<SparseComplexMatrix>::Q (void) const;
+    template OCTAVE_API ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::Q (bool econ) const;
 
     Matrix
     qrsolve (const SparseMatrix& a, const MArray<double>& b,
--- a/liboctave/numeric/sparse-qr.h	Mon Jan 04 14:52:06 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 <typename RHS_T, typename RET_T>
       static OCTAVE_API RET_T
       solve (const SPARSE_T& a, const RHS_T& b,
@@ -85,6 +101,11 @@
 
     private:
 
+      template <typename RHS_T,typename RET_T>
+      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;
--- a/liboctave/util/oct-sparse.cc	Mon Jan 04 14:52:06 2021 +0100
+++ b/liboctave/util/oct-sparse.cc	Thu Jan 09 21:45:54 2020 +0100
@@ -78,6 +78,204 @@
 
     return reinterpret_cast<const octave_idx_type *> (i);
   }
+
+#  if defined (HAVE_CHOLMOD)
+  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<suitesparse_integer*> (to_suitesparse_intptr (a.cidx ()));
+    A.i = const_cast<suitesparse_integer*> (to_suitesparse_intptr (a.ridx ()));
+    A.x = const_cast<double*> (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<suitesparse_integer*> (to_suitesparse_intptr (a.cidx ()));
+    A.i = const_cast<suitesparse_integer*> (to_suitesparse_intptr (a.ridx ()));
+    A.x = const_cast<Complex*>
+                 (reinterpret_cast<const Complex*> (a.data ()));
+
+    return A;
+  }
+
+  cholmod_dense*
+  rod2ccd (const MArray<double> &a, cholmod_common* cc1)
+  {
+    cholmod_dense *A = static_cast<cholmod_dense*>
+                         (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<Complex*> (A->x))[j] = Complex (a_x[j], 0.0);
+
+    return A;
+  }
+
+  const cholmod_dense
+  rod2rcd (const MArray<double> &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<double*> (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<Complex*>
+            (reinterpret_cast<const Complex*> (a.data ()));
+
+    return A;
+  }
+
+  SparseMatrix
+  rcs2ros (const cholmod_sparse* y)
+  {
+    SparseMatrix ret (static_cast<octave_idx_type> (y->nrow),
+                      static_cast<octave_idx_type> (y->ncol),
+                      static_cast<octave_idx_type> (y->nzmax));
+
+    octave_idx_type nz = static_cast<octave_idx_type> (y->nzmax);
+
+    for (octave_idx_type j = 0; j < static_cast<octave_idx_type> (y->ncol) + 1;
+         j++)
+      ret.xcidx (j) = (static_cast<octave_idx_type*> (y->p))[j];
+
+    for (octave_idx_type j = 0; j < nz; j++)
+      {
+        ret.xridx (j) = (static_cast<octave_idx_type*> (y->i))[j];
+        ret.xdata (j) = (static_cast<double*> (y->x))[j];
+      }
+
+    return ret;
+  }
+
+  SparseComplexMatrix
+  ccs2cos (const cholmod_sparse* a)
+  {
+    SparseComplexMatrix ret (static_cast<octave_idx_type> (a->nrow),
+                             static_cast<octave_idx_type> (a->ncol),
+                             static_cast<octave_idx_type> (a->nzmax));
+
+    octave_idx_type nz = static_cast<octave_idx_type> (a->nzmax);
+
+    for (octave_idx_type j = 0; j < static_cast<octave_idx_type> (a->ncol) + 1;
+         j++)
+      ret.xcidx (j) = (static_cast<octave_idx_type*> (a->p))[j];
+
+    for (octave_idx_type j = 0; j < nz; j++)
+      {
+        ret.xridx (j) = (static_cast<octave_idx_type*> (a->i))[j];
+        ret.xdata (j) = (static_cast<Complex*> (a->x))[j];
+      }
+ 
+    return ret;
+  }
+
+  cholmod_sparse*
+  ros2ccs (const SparseMatrix& a, cholmod_common* cc1)
+  {
+    cholmod_sparse *A = static_cast<cholmod_sparse*>
+                          (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<octave_idx_type> (a.cols ()) + 1; j++)
+      static_cast<octave_idx_type*> (A->p)[j] = a.cidx(j);
+
+    const double *a_x = a.data ();
+    for (octave_idx_type j = 0; j < a.nnz (); j++)
+      {
+        (static_cast<Complex*> (A->x))[j] = Complex (a_x[j], 0.0);
+        (static_cast<octave_idx_type*> (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
 }
 
 #endif
--- a/liboctave/util/oct-sparse.h	Mon Jan 04 14:52:06 2021 +0100
+++ b/liboctave/util/oct-sparse.h	Thu Jan 09 21:45:54 2020 +0100
@@ -28,6 +28,11 @@
 
 #include "octave-config.h"
 
+#if defined (HAVE_CHOLMOD)
+#  include "dSparse.h"
+#  include "CSparse.h"
+#endif
+
 #if defined (HAVE_SUITESPARSE_AMD_H)
 #  include <suitesparse/amd.h>
 #elif defined (HAVE_UFSPARSE_AMD_H)
@@ -88,6 +93,10 @@
 #  include <umfpack.h>
 #endif
 
+#if defined (HAVE_SUITESPARSE_SUITESPARSEQR_HPP)
+#  include <suitesparse/SuiteSparseQR.hpp>
+#endif
+
 // Cope with new SuiteSparse versions
 
 #if defined (SUITESPARSE_VERSION)
@@ -161,7 +170,8 @@
 
 #if (defined (HAVE_AMD) || defined (HAVE_CCOLAMD)               \
      || defined (HAVE_CHOLMOD) || defined (HAVE_COLAMD)         \
-     || defined (HAVE_CXSPARSE) || defined (HAVE_UMFPACK))
+     || defined (HAVE_CXSPARSE) || defined (HAVE_SPQR)          \
+     || defined (HAVE_UMFPACK))
 
 namespace octave
 {
@@ -182,6 +192,53 @@
 
   extern OCTAVE_API const octave_idx_type*
   to_octave_idx_type_ptr (const suitesparse_integer *i);
+
+#  if defined (HAVE_CHOLMOD)
+  // 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<double>& 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<double> &a, cholmod_common *cc1);
+
+  extern OCTAVE_API void
+  spqr_error_handler (const cholmod_common *cc);
+#  endif
 }
 
 #endif