# HG changeset patch
# User Simon Hau <simon.hau79@gmail.com>
# Date 1578602754 -3600
#      Thu Jan 09 21:45:54 2020 +0100
# Node ID f2b2c2745c6263499c570e579df30b5fbba5cb7b
# Parent  baf16e6f498b866d1d62feeb1dd838cf5191bd45
Replacement of CXSPARSE by SPQR for QR factorization (bug #57033).

* qr.cc (Fqr): Handle sparse input. Return economy versions of C, Q and R for
  qr(..., 0). Retun column permutation of A as (vector|matrix) for
  qr(..., ("vector"|"matrix")). Additional tests if HAVE_SPQR.
* oct-sparse.[cc,h] (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.
* sparse-qr.[cc,h]: Add interface to SPQR functios.
* configure.ac: Add configure tests for SPQR.

diff --git a/configure.ac b/configure.ac
--- a/configure.ac
+++ b/configure.ac
@@ -2118,6 +2118,16 @@
   [], [don't use CHOLMOD library, disable some sparse matrix functionality])
 LIBS="$save_LIBS"
 
+### Check for SPQR library
+### FIXME: We actually need to check for SuiteSparseQR_min2norm<double>.  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,
@@ -2183,11 +2193,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)
@@ -3115,6 +3125,9 @@
   Sndfile CPPFLAGS:              $SNDFILE_CPPFLAGS
   Sndfile LDFLAGS:               $SNDFILE_LDFLAGS
   Sndfile libraries:             $SNDFILE_LIBS
+  SPQR CPPFLAGS:                 $SPQR_CPPFLAGS
+  SPQR LDFLAGS:                  $SPQR_LDFLAGS
+  SPQR libraries:                $SPQR_LIBS
   SuiteSparse config libraries:  $SUITESPARSECONFIG_LIBS
   SUNDIALS IDA CPPFLAGS:         $SUNDIALS_IDA_CPPFLAGS
   SUNDIALS IDA LDFLAGS:          $SUNDIALS_IDA_LDFLAGS
diff --git a/libinterp/corefcn/qr.cc b/libinterp/corefcn/qr.cc
--- a/libinterp/corefcn/qr.cc
+++ b/libinterp/corefcn/qr.cc
@@ -331,45 +331,161 @@
 
   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<std::complex<double>>, 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<std::complex<double>>, 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
@@ -896,7 +1012,7 @@
 %! 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 @@
 %! 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 @@
 %! [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 --git a/liboctave/numeric/sparse-qr.cc b/liboctave/numeric/sparse-qr.cc
--- a/liboctave/numeric/sparse-qr.cc
+++ b/liboctave/numeric/sparse-qr.cc
@@ -47,15 +47,18 @@
     template <typename SPARSE_T>
     class
     cxsparse_types
-    {
-    };
+      { };
 
     template <>
     class
     cxsparse_types<SparseMatrix>
     {
     public:
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+      typedef CXSPARSE_DNAME (s) symbolic_type;
+      typedef CXSPARSE_DNAME (n) numeric_type;
+      typedef Matrix dense_matrix_type;
+#elif defined (HAVE_CXSPARSE)
       typedef CXSPARSE_DNAME (s) symbolic_type;
       typedef CXSPARSE_DNAME (n) numeric_type;
 #else
@@ -69,7 +72,10 @@
     cxsparse_types<SparseComplexMatrix>
     {
     public:
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+      typedef CXSPARSE_ZNAME (s) symbolic_type;
+      typedef CXSPARSE_ZNAME (n) numeric_type;
+#elif defined (HAVE_CXSPARSE)
       typedef CXSPARSE_ZNAME (s) symbolic_type;
       typedef CXSPARSE_ZNAME (n) numeric_type;
 #else
@@ -95,7 +101,9 @@
 
       bool ok (void) const
       {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+        return (H && Htau && HPinv && R_ && E_ && cc);
+#elif defined (HAVE_CXSPARSE)
         return (N && S);
 #else
         return false;
@@ -108,21 +116,42 @@
 
       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<octave_idx_type> count;
 
       octave_idx_type nrows;
       octave_idx_type ncols;
+      // octave_idx_type econ;
+
+#if defined (HAVE_SPQR)
 
       typename cxsparse_types<SPARSE_T>::symbolic_type *S;
       typename cxsparse_types<SPARSE_T>::numeric_type *N;
+      template <typename RHS_T, typename RET_T>
+      RET_T
+      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
@@ -131,6 +160,21 @@
       template <typename RHS_T, typename RET_T>
       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 <typename SPARSE_T>
@@ -157,7 +201,16 @@
     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 (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 +226,77 @@
 #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> (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)
+      , 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<double> (order, SPQR_DEFAULT_TOL, A.nrow,
+                             const_cast<cholmod_sparse*> (&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 +320,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<SparseMatrix>::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 +357,11 @@
     SparseMatrix
     sparse_qr<SparseMatrix>::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 +400,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> (R_->nrow);
+      octave_idx_type nc = static_cast<octave_idx_type> (R_->ncol);
+      octave_idx_type nz = static_cast<octave_idx_type> (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*> (R_->p));
+      octave_idx_type *Ri = to_octave_idx_type_ptr
+                              (static_cast<suitesparse_integer*> (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*> (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 +464,45 @@
 
     template <>
     Matrix
-    sparse_qr<SparseMatrix>::sparse_qr_rep::C (const Matrix& b) const
+    sparse_qr<SparseMatrix>::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<double> (SPQR_QTX, H, Htau,HPinv,
+                                         const_cast<cholmod_dense*>(&B), cc1);
+      spqr_error_handler (cc1);
+      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, 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 +566,156 @@
 
     template <>
     Matrix
+    sparse_qr<SparseMatrix>::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<double> (SPQR_QTX, H, Htau,HPinv,
+                                         const_cast<cholmod_dense*>(&B), cc1);
+      spqr_error_handler (cc1);
+      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, cc1);
+
+      CHOLMOD_NAME(finish) (cc1);
+
+      return ret;
+
+#else
+
+      octave_unused_parameter (b);
+
+      return Matrix ();
+
+#endif
+    }
+
+    template <>
+    Matrix
+    sparse_qr<SparseMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                             (nrows, nrows, nrows, CHOLMOD_REAL, cc1));
+
+      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, H, Htau, HPinv, I, cc1);
+      spqr_error_handler (cc1);
+      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, cc1);
+      CHOLMOD_NAME (free_dense) (&I, cc1);
+      CHOLMOD_NAME (finish) (cc1);
+
+      return ret;
+
+#else
+
+      return Matrix ();
+
+#endif
+    }
+
+    template <>
+    Matrix
     sparse_qr<SparseMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                             (nrows, nrows, nrows, CHOLMOD_REAL, cc1));
+
+      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, H, Htau, HPinv, I, cc1);
+      spqr_error_handler (cc1);
+      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, 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 +755,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;
             }
@@ -447,7 +778,66 @@
     {
       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<double>
+             (SPQR_QTX, H, Htau, HPinv, const_cast<cholmod_dense*> (&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<double*> (R_->x);
+      R2.i = static_cast<suitesparse_integer*> (R_->i);
+      R2.p = static_cast<suitesparse_integer*> (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 (R_\(Q'*B(:,j)) and store result in QTB(:,j)
+          CXSPARSE_DNAME (_usolve)
+            (&R2, &(static_cast<double*> (QTB->x)[j * b_nr]));
+          // x(:,j) = E_' * (R_\(Q'*B(:,j))
+          CXSPARSE_DNAME (_ipvec)
+            (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, 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 +899,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 +951,6 @@
       return x;
 
 #else
-
       octave_unused_parameter (b);
 
       return Matrix ();
@@ -949,12 +1337,38 @@
     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)
+        , 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<Complex> (order, SPQR_DEFAULT_TOL, A.nrow,
+                              const_cast<cholmod_sparse*>(&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 +1393,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) (&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 +1468,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> (R_->nrow);
+      octave_idx_type nc = static_cast<octave_idx_type> (R_->ncol);
+      octave_idx_type nz = static_cast<octave_idx_type> (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*> (R_->p));
+      octave_idx_type *Ri = to_octave_idx_type_ptr
+                              (static_cast<suitesparse_integer*> (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*> (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 +1535,48 @@
     ComplexMatrix
     sparse_qr<SparseComplexMatrix>::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<Complex> (SPQR_QTX, H, Htau, HPinv,
+                                          const_cast<cholmod_dense*> (&B),
+                                          cc1);
+      spqr_error_handler (cc1);
+      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, 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 +1637,104 @@
 
     template <>
     ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::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<Complex> (SPQR_QTX, H, Htau, HPinv,
+                                          const_cast<cholmod_dense*> (&B),
+                                          cc1);
+      spqr_error_handler (cc1);
+      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, cc1);
+
+      CHOLMOD_NAME (finish) (cc1);
+
+      return ret;
+
+#else
+
+      octave_unused_parameter (b);
+
+      return ComplexMatrix ();
+
+#endif
+    }
+
+
+    template <>
+    ComplexMatrix
     sparse_qr<SparseComplexMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                              (nrows, nrows, nrows, CHOLMOD_COMPLEX, cc1));
+
+      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, H, Htau, HPinv, I, cc1);
+      spqr_error_handler (cc1);
+      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, 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 +1795,57 @@
     }
 
     template <>
+    ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::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_dense*>
+                           (CHOLMOD_NAME (allocate_dense)
+                              (nrows, nrows, nrows, CHOLMOD_COMPLEX, cc1));
+
+      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, H, Htau, HPinv, I, cc1);
+      spqr_error_handler (cc1);
+      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, cc1);
+      CHOLMOD_NAME (free_dense) (&I, cc1);
+      CHOLMOD_NAME (finish) (cc1);
+
+      return ret;
+
+#else
+
+      return ComplexMatrix ();
+
+#endif
+    }
+
+    template <>
     template <>
     SparseComplexMatrix
     sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseComplexMatrix,
@@ -2188,6 +2828,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 +2866,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 +2903,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 +2915,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 +2928,33 @@
     sparse_qr<SPARSE_T>::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<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 octave::math::sparse_qr<SPARSE_T>::min2norm_solve<RHS_T, RET_T>
+               (a, b, info, order);
+
+#elif defined (HAVE_CXSPARSE)
+
       info = -1;
 
       octave_idx_type nr = a.rows ();
@@ -2277,8 +2987,382 @@
 
           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
+    Matrix
+    sparse_qr<SparseMatrix>::solve<MArray<double>, Matrix>
+    (const SparseMatrix& a, const MArray<double>& b, octave_idx_type& info);
+
+    template
+    SparseMatrix
+    sparse_qr<SparseMatrix>::solve<SparseMatrix, SparseMatrix>
+    (const SparseMatrix& a, const SparseMatrix& b, octave_idx_type& info);
+
+    template
+    ComplexMatrix
+    sparse_qr<SparseMatrix>::solve<MArray<std::complex<double>>, ComplexMatrix>
+    (const SparseMatrix& a, const MArray<std::complex<double>>& b,
+     octave_idx_type& info);
+
+    template
+    SparseComplexMatrix
+    sparse_qr<SparseMatrix>::solve<SparseComplexMatrix, SparseComplexMatrix>
+    (const SparseMatrix& a, const SparseComplexMatrix& b,
+     octave_idx_type& info);
+
+    template
+    ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<
+    MArray<std::complex<double>>, ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<std::complex<double>>& b,
+     octave_idx_type& info);
+
+    template
+    SparseComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<
+    SparseComplexMatrix, SparseComplexMatrix>
+    (const SparseComplexMatrix& a, const SparseComplexMatrix& b,
+     octave_idx_type& info);
+
+    template
+    ComplexMatrix
+    sparse_qr<SparseComplexMatrix>::solve<MArray<double>, ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<double>& b,
+     octave_idx_type& info);
+
+    template
+    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
+    SparseMatrix
+    sparse_qr<SparseMatrix>::E_MAT (void) const;
+
+    template
+    SparseMatrix
+    sparse_qr<SparseComplexMatrix>::E_MAT (void) const;
+
+    //specializations of function min2norm_solve
+    template <>
+    template <>
+    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 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<double> (order, SPQR_DEFAULT_TOL,
+                                          const_cast<cholmod_sparse*> (&A),
+                                          const_cast<cholmod_dense*> (&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<double*> (X->x)[i];
+
+      info = 0;
+      CHOLMOD_NAME (finish) (cc1);
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    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 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<double> (order, SPQR_DEFAULT_TOL,
+                                          const_cast<cholmod_sparse*>(&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<SparseMatrix>::min2norm_solve<MArray<std::complex<double>>,
+                                            ComplexMatrix>
+    (const SparseMatrix& a,const MArray<std::complex<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 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<Complex> (order, SPQR_DEFAULT_TOL, A,
+                                           const_cast<cholmod_dense*> (&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<Complex*> (X->x)[i];
+
+      CHOLMOD_NAME (free_sparse) (&A, cc1);
+      CHOLMOD_NAME (finish) (cc1);
+
+      info = 0;
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    SparseComplexMatrix
+    sparse_qr<SparseMatrix>::min2norm_solve<SparseComplexMatrix,
+                                            SparseComplexMatrix>
+    (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<Complex> (order, SPQR_DEFAULT_TOL, A,
+                                           const_cast<cholmod_sparse*> (&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<SparseComplexMatrix>::min2norm_solve<MArray<std::complex<double>>,
+                                                   ComplexMatrix>
+    (const SparseComplexMatrix& a, const MArray<std::complex<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 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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&A),
+                                           const_cast<cholmod_dense*> (&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<Complex*> (X->x)[i];
+
+      CHOLMOD_NAME (finish) (cc1);
+
+      info = 0;
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    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 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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&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<Complex*> (X->x)[i];
+
+      CHOLMOD_NAME (free_dense) (&B, cc1);
+      CHOLMOD_NAME (finish) (cc1);
+
+      info = 0;
+
+      return x;
+
+    }
+
+    template <>
+    template <>
+    SparseComplexMatrix
+    sparse_qr<SparseComplexMatrix>::min2norm_solve<SparseComplexMatrix,
+                                                   SparseComplexMatrix>
+    (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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&A),
+                                           const_cast<cholmod_sparse*> (&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<SparseComplexMatrix>::min2norm_solve<SparseMatrix,
+                                                   SparseComplexMatrix>
+    (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<Complex> (order, SPQR_DEFAULT_TOL,
+                                           const_cast<cholmod_sparse*> (&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 <typename SPARSE_T>
     template <typename RHS_T, typename RET_T>
     RET_T
diff --git a/liboctave/numeric/sparse-qr.h b/liboctave/numeric/sparse-qr.h
--- a/liboctave/numeric/sparse-qr.h
+++ b/liboctave/numeric/sparse-qr.h
@@ -54,7 +54,12 @@
 
       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 +69,11 @@
 
       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 +86,14 @@
       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 <typename RHS_T, typename RET_T>
       static RET_T
       solve (const SPARSE_T& a, const RHS_T& b,
@@ -85,6 +101,11 @@
 
     private:
 
+      template <typename RHS_T,typename RET_T>
+      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 --git a/liboctave/util/oct-sparse.cc b/liboctave/util/oct-sparse.cc
--- a/liboctave/util/oct-sparse.cc
+++ b/liboctave/util/oct-sparse.cc
@@ -78,6 +78,204 @@
 
     return reinterpret_cast<const octave_idx_type *> (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<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<cs_complex_t *>
+                 (reinterpret_cast<const cs_complex_t *> (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<cs_complex_t *>
+            (reinterpret_cast<const cs_complex_t *> (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
diff --git a/liboctave/util/oct-sparse.h b/liboctave/util/oct-sparse.h
--- a/liboctave/util/oct-sparse.h
+++ b/liboctave/util/oct-sparse.h
@@ -28,6 +28,9 @@
 
 #include "octave-config.h"
 
+#include "dSparse.h"
+#include "CSparse.h"
+
 #if defined (HAVE_SUITESPARSE_AMD_H)
 #  include <suitesparse/amd.h>
 #elif defined (HAVE_UFSPARSE_AMD_H)
@@ -88,6 +91,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)
@@ -182,6 +189,51 @@
 
   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<double>& 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<double> &a, cholmod_common *cc1);
+
+  extern void
+  spqr_error_handler (const cholmod_common *cc);
 }
 
 #endif
