# HG changeset patch
# User simon.hau79@gmail.com
# Date 1570380260 -7200
#      Sun Oct 06 18:44:20 2019 +0200
# Node ID c91ea5884eddaf7ec1879688e5774379840d8fde
# Parent  b470734a025ea552044418630fb25a2d92856032
test patch for replacement of cxsparse by spqr

diff -r b470734a025e -r c91ea5884edd libinterp/dldfcn/qr.cc
--- a/libinterp/dldfcn/qr.cc	Sat Oct 05 14:04:36 2019 -0700
+++ b/libinterp/dldfcn/qr.cc	Sun Oct 06 18:44:20 2019 +0200
@@ -328,49 +328,116 @@ orthogonal basis of @code{span (A)}.
         is_cmplx = true;
     }
 
-  if (arg.issparse ())
-    {
-      if (nargout > 2)
-        error ("qr: Permutation output is not supported for sparse input");
-
-      if (is_cmplx)
-        {
-          octave::math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ());
-
-          if (have_b)
-            {
-              retval = ovl (q.C (args(1).complex_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 ());
-            }
-          else if (nargout > 1)
-            retval = ovl (q.Q (), q.R (economy));
-          else
-            retval = ovl (q.R (economy));
-        }
-      else
-        {
-          octave::math::sparse_qr<SparseMatrix> q (arg.sparse_matrix_value ());
-
-          if (have_b)
-            {
-              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 ());
-            }
-          else if (nargout > 1)
-            retval = ovl (q.Q (), q.R (economy));
-          else
-            retval = ovl (q.R (economy));
-        }
-    }
+    if (arg.issparse ())
+      {
+	if (nargout > 3)
+	  error ("qr: too many output arguments");
+	if (is_cmplx)
+	  {
+	    if (have_b)
+	      {
+		octave_idx_type info;
+		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
+	      {
+		octave::math::sparse_qr<SparseComplexMatrix> q (
+		    arg.sparse_complex_matrix_value ());
+		if (nargout > 2)
+		  retval = ovl (q.Q (), q.R (economy), q.E ());
+		else if (nargout > 1)
+		  retval = ovl (q.Q (), q.R (economy));
+		else
+		  retval = ovl (q.R (economy));
+	      }
+	  }
+	else
+	  {
+	    if (have_b)
+	      {
+		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
+	      {
+		octave::math::sparse_qr<SparseMatrix> q (
+		    arg.sparse_matrix_value ());
+		if (nargout > 2)
+		  retval = ovl (q.Q (), q.R (economy), q.E ());
+		else if (nargout > 1)
+		  retval = ovl (q.Q (), q.R (economy));
+		else
+		  retval = ovl (q.R (economy));
+	      }
+	  }
+      }
   else
     {
       if (arg.is_single_type ())
@@ -1243,9 +1310,12 @@ If @var{orient} is @qcode{"row"}, full f
       if (argq.is_single_type () || argr.is_single_type ()
           || argx.is_single_type ())
         {
-          FloatComplexMatrix Q = argq.float_complex_matrix_value ();
-          FloatComplexMatrix R = argr.float_complex_matrix_value ();
-          FloatComplexMatrix x = argx.float_complex_matrix_value ();
+          FloatComplexMatrix Q =
+            argq.float_complex_matrix_value ();
+          FloatComplexMatrix R =
+            argr.float_complex_matrix_value ();
+          FloatComplexMatrix x =
+            argx.float_complex_matrix_value ();
 
           octave::math::qr<FloatComplexMatrix> fact (Q, R);
 
@@ -1432,8 +1502,10 @@ If @var{orient} is @qcode{"row"}, full f
       // complex case
       if (argq.is_single_type () || argr.is_single_type ())
         {
-          FloatComplexMatrix Q = argq.float_complex_matrix_value ();
-          FloatComplexMatrix R = argr.float_complex_matrix_value ();
+          FloatComplexMatrix Q =
+            argq.float_complex_matrix_value ();
+          FloatComplexMatrix R =
+            argr.float_complex_matrix_value ();
 
           octave::math::qr<FloatComplexMatrix> fact (Q, R);
 
diff -r b470734a025e -r c91ea5884edd liboctave/numeric/sparse-qr.cc
--- a/liboctave/numeric/sparse-qr.cc	Sat Oct 05 14:04:36 2019 -0700
+++ b/liboctave/numeric/sparse-qr.cc	Sun Oct 06 18:44:20 2019 +0200
@@ -1,5 +1,4 @@
 /*
-
 Copyright (C) 2016-2019 John W. Eaton
 Copyright (C) 2005-2019 David Bateman
 
@@ -24,7 +23,7 @@ along with Octave; see the file COPYING.
 #if defined (HAVE_CONFIG_H)
 #  include "config.h"
 #endif
-
+#include <iostream>
 #include "CMatrix.h"
 #include "CSparse.h"
 #include "MArray.h"
@@ -38,6 +37,7 @@ along with Octave; see the file COPYING.
 #include "quit.h"
 #include "sparse-qr.h"
 
+
 namespace octave
 {
   namespace math
@@ -53,7 +53,11 @@ namespace octave
     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
@@ -67,7 +71,10 @@ namespace octave
     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
@@ -80,9 +87,8 @@ namespace octave
     class sparse_qr<SPARSE_T>::sparse_qr_rep
     {
     public:
-
-      sparse_qr_rep (const SPARSE_T& a, int order);
-
+       sparse_qr_rep(const SPARSE_T& a, int order);
+      
       // No copying!
 
       sparse_qr_rep (const sparse_qr_rep&) = delete;
@@ -93,7 +99,9 @@ namespace octave
 
       bool ok (void) const
       {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+	return (H && Htau && HPinv && R_ && E_ && cc);
+#elif defined (HAVE_CXSPARSE)
         return (N && S);
 #else
         return false;
@@ -105,6 +113,8 @@ namespace octave
       ColumnVector Pinv (void) const;
 
       ColumnVector P (void) const;
+      
+      ColumnVector E (void) const;
 
       SPARSE_T R (bool econ) const;
 
@@ -114,13 +124,21 @@ namespace octave
       typename SPARSE_T::dense_matrix_type
       Q (void) const;
 
-      refcount<octave_idx_type> count;
-
+      refcount<int> 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;
+#else
+#endif
 
       template <typename RHS_T, typename RET_T>
       RET_T
@@ -129,8 +147,20 @@ namespace octave
       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
+#elif defined (HAVE_CXSPARSE)
+#else
+#endif
     };
-
+    
     template <typename SPARSE_T>
     ColumnVector
     sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const
@@ -155,7 +185,13 @@ namespace octave
     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);
 
@@ -170,21 +206,73 @@ namespace octave
 
 #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;
+#elif defined (HAVE_CXSPARSE)
+    return ColumnVector ();
+#else
+      return ColumnVector ();
+#endif
+    }
+    
     // Specializations.
 
     // Real-valued matrices.
 
-    template <>
-    sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep
-    (const SparseMatrix& a, int order)
-      : count (1), nrows (a.rows ()), ncols (a.columns ())
-#if defined (HAVE_CXSPARSE)
-      , S (nullptr), N (nullptr)
-#endif
-    {
-#if defined (HAVE_CXSPARSE)
-
+#if defined (HAVE_SPQR)
+       /*
+       * 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 ()),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)
+   template<>
+   sparse_qr<SparseMatrix>::sparse_qr_rep::sparse_qr_rep (
+    const SparseMatrix &a, int order) :
+    count (1), nrows (a.rows ()), ncols (a.columns ())
+    , S (nullptr), N (nullptr)
+   {
       CXSPARSE_DNAME () A;
 
       A.nzmax = a.nnz ();
@@ -208,6 +296,8 @@ namespace octave
         (*current_liboctave_error_handler)
           ("sparse_qr: sparse matrix QR factorization filled");
 
+      count = 1;
+
 #else
 
       octave_unused_parameter (order);
@@ -216,22 +306,32 @@ namespace octave
         ("sparse_qr: support for CXSparse was unavailable or disabled when liboctave was built");
 
 #endif
-    }
-
-    template <>
+      }
+
+    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);
+#else
 #endif
     }
 
-    template <>
+    template<>
     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?
@@ -266,11 +366,29 @@ namespace octave
 #endif
     }
 
-    template <>
+    template<>
     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?
@@ -306,13 +424,45 @@ namespace octave
       return SparseMatrix ();
 
 #endif
-    }
-
-    template <>
+      }
+
+    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) ("matrix dimension mismatch");
+      else if (b_nc <= 0 || b_nr <= 0)
+	(*current_liboctave_error_handler)
+	("matrix dimension with negative or zero size");
+      else
+	{
+	  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 ();
@@ -372,17 +522,45 @@ namespace octave
       return Matrix ();
 
 #endif
-    }
-
-    template <>
+      }
+
+    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");
@@ -422,7 +600,7 @@ namespace octave
                 }
 
               for (octave_idx_type i = 0; i < nr; i++)
-                vec[i+idx] = buf[i];
+                ret_vec[i+idx] = buf[i];
 
               bvec[j] = 0.0;
             }
@@ -437,15 +615,66 @@ namespace octave
 #endif
     }
 
-    template <>
-    template <>
+
+    template<>
+    template<>
     Matrix
-    sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>, Matrix>
-    (const MArray<double>& b, octave_idx_type& info) const
+    sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<double>,
+	Matrix> (const MArray<double> &b, octave_idx_type &info) const
     {
       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");
+      else
+	{
+	  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;
@@ -500,15 +729,70 @@ namespace octave
 #endif
     }
 
-    template <>
-    template <>
+    template<>
+    template<>
     Matrix
-    sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>, Matrix>
-    (const MArray<double>& b, octave_idx_type& info) const
+    sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<double>,
+	Matrix> (const MArray<double> &b, octave_idx_type &info) const
     {
-      info = -1;
-
-#if defined (HAVE_CXSPARSE)
+	  info = -1;
+#if defined (HAVE_SPQR)
+	  octave_idx_type nr = ncols;
+	  octave_idx_type nc = nrows;
+
+	  cholmod_common Common1;
+	  cholmod_common *cc1 = &Common1;
+	  CHOLMOD_NAME(start) (cc1);
+	  octave_idx_type b_nc = b.cols ();
+	  octave_idx_type b_nr = b.rows ();
+	  const double *bvec = b.data ();
+	  Matrix x (nc, b_nc); // X = Q * (R' \ E' * b)
+	  double *vec = x.fortran_vec ();
+	  cholmod_dense *tmp_X;
+	  cholmod_dense *tmp_B =
+	      static_cast<cholmod_dense*> (CHOLMOD_NAME(allocate_dense) (
+		  nc, b_nc, nc, CHOLMOD_REAL, cc1));
+	  OCTAVE_LOCAL_BUFFER(double, buf, nc);
+	  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;
+	  for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc;
+	      i++, idx += nc, bidx += b_nr)
+	    {
+	      octave_quit ();
+
+	      for (volatile octave_idx_type j = nr; j < nc; j++)
+		buf[j] = 0.;
+
+	      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+		CXSPARSE_DNAME (_ipvec) (E_, bvec + bidx, buf, nr);
+		CXSPARSE_DNAME (_utsolve) (&R2, buf);
+		END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	      for (volatile octave_idx_type j = 0; j < nc; j++)
+		{
+		  static_cast<double*> (tmp_B->x)[idx + j] = buf[j];
+		}
+	    }
+	  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	  tmp_X = SuiteSparseQR_qmult<double> (SPQR_QX, H, Htau, HPinv, tmp_B,
+					       cc1);
+	  END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	  spqr_error_handler (cc1);
+	  for (volatile octave_idx_type i = 0; i < nc * b_nc; i++)
+	    {
+	      vec[i] = static_cast<double*> (tmp_X->x)[i];
+	    }
+	  CHOLMOD_NAME(free_dense) (&tmp_X, cc1);
+	  CHOLMOD_NAME(free_dense) (&tmp_B, cc1);
+	  CHOLMOD_NAME(finish) (cc1);
+	  info = 0;
+	  return x;
+#elif defined (HAVE_CXSPARSE)
 
       // These are swapped because the original matrix was transposed in
       // sparse_qr<SparseMatrix>::solve.
@@ -568,15 +852,109 @@ namespace octave
 #endif
     }
 
-    template <>
-    template <>
-    SparseMatrix
-    sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix, SparseMatrix>
-    (const SparseMatrix& b, octave_idx_type& info) const
-    {
-      info = -1;
-
-#if defined (HAVE_CXSPARSE)
+    template<>
+      template<>
+	SparseMatrix
+	sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<SparseMatrix,
+	    SparseMatrix> (const SparseMatrix &b, octave_idx_type &info) const
+	{
+	  info = -1;
+
+#if defined (HAVE_SPQR)
+	  octave_idx_type nr = nrows;
+	  octave_idx_type nc = ncols;
+	  octave_idx_type b_nr = b.rows ();
+	  octave_idx_type b_nc = b.cols ();
+
+	  SparseMatrix x (nc, b_nc, b.nnz ()); //X = E*(R\(Q'*B))
+	  x.xcidx (0) = 0;
+	  volatile octave_idx_type x_nz = b.nnz ();
+	  volatile octave_idx_type ii = 0;
+	  OCTAVE_LOCAL_BUFFER(double, buf, nc);
+	  OCTAVE_LOCAL_BUFFER(double, Xx, nc);
+	  cholmod_common Common1;
+	  cholmod_common *cc1 = &Common1;
+	  CHOLMOD_NAME(start) (cc1);
+	  if (nc > nr)
+	    (*current_liboctave_error_handler) (
+		"nr columns > nr rows is not allowed in tall solve");
+	  if (nrows < 0 || ncols < 0 || nrows != b_nr)
+	    (*current_liboctave_error_handler) ("matrix dimension mismatch");
+	  if (nrows == 0 || ncols == 0 || b_nc == 0)
+	    (*current_liboctave_error_handler) (
+		"matrix dimension with size zero");
+	  else
+	    {
+	      cholmod_sparse *QTB; //Q' * B
+	      const cholmod_sparse tmp_B = ros2rcs (b);
+	      //convert R_ into CXSARSE matrix
+	      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;
+	      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+		QTB = SuiteSparseQR_qmult<double> (SPQR_QTX, H, Htau, HPinv,
+						   const_cast<cholmod_sparse*>(&tmp_B), cc1);
+		END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	      suitesparse_integer col_start_idx;
+	      suitesparse_integer col_end_idx;
+	      for (volatile octave_idx_type i = 0; i < b_nc; i++)
+		{
+		  //buf = QTB(:,i)
+		  col_start_idx = static_cast<suitesparse_integer*> (QTB->p)[i];
+		  col_end_idx =
+		      static_cast<suitesparse_integer*> (QTB->p)[i + 1];
+		  for (volatile octave_idx_type j = 0; j < nc; j++)
+		    {
+		      buf[j] = 0.0;
+		      Xx[j] = 0.0;
+		    }
+		  suitesparse_integer k;
+		  for (volatile suitesparse_integer j = col_start_idx;
+		      j < col_end_idx; j++)
+		    {
+		      k = static_cast<suitesparse_integer*> (QTB->i)[j];
+		      if (k < nc)
+			buf[k] = static_cast<double*> (QTB->x)[j];
+		    }
+		  BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+		    CXSPARSE_DNAME (_usolve) (&R2, buf);
+		    CXSPARSE_DNAME (_ipvec) (E_, buf, Xx, nc);
+		    END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+		  for (volatile octave_idx_type j = 0; j < nc; j++)
+		    {
+		      double tmp = Xx[j];
+
+		      if (tmp != 0.0)
+			{
+			  if (ii == x_nz)
+			    {
+			      // Resize the sparse matrix
+			      octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+			      sz = (sz > 10 ? sz : 10) + x_nz;
+			      x.change_capacity (sz);
+			      x_nz = sz;
+			    }
+
+			  x.xdata (ii) = tmp;
+			  x.xridx (ii++) = j;
+			}
+		    }
+		  x.xcidx (i + 1) = ii;
+		}
+	      CHOLMOD_NAME(free_sparse) (&QTB, cc1);
+	      CHOLMOD_NAME(finish) (cc1);
+	      info = 0;
+	      return x;
+	    }
+
+
+
+#elif defined (HAVE_CXSPARSE)
 
       octave_idx_type nr = nrows;
       octave_idx_type nc = ncols;
@@ -664,8 +1042,107 @@ namespace octave
     (const SparseMatrix& b, octave_idx_type& info) const
     {
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+
+      octave_idx_type nr = ncols;
+      octave_idx_type nc = nrows;
+      octave_idx_type b_nc = b.cols ();
+
+      SparseMatrix x (nc, b_nc, b.nnz ()); // X = Q * (R' \ E_' * b)
+      x.xcidx (0) = 0;
+
+      octave_idx_type col_start_idx;
+      octave_idx_type col_end_idx;
+      octave_idx_type nz_counter;
+
+      //sparse column vector for storing intermediate result
+      cholmod_sparse *S;
+      cholmod_sparse * QS; //Q * S
+
+      cholmod_common Common1;
+      cholmod_common *cc1 = &Common1;
+      CHOLMOD_NAME(start) (cc1);
+
+      volatile octave_idx_type x_nz = b.nnz ();
+      volatile octave_idx_type ii = 0;
+      //FIXME free b_column and buf?
+      //b_column contains column of matrix b
+      OCTAVE_LOCAL_BUFFER (double,b_column,nr);
+      //buf stores intermediate result
+      OCTAVE_LOCAL_BUFFER (double,buf,nr);
+      //convert R_ to a CXSPARSE matrix
+      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;
+      //calculate x(:,i)
+      for (volatile octave_idx_type i = 0; i < b_nc;i++)
+	 {
+	   octave_quit ();
+	   for (octave_idx_type j = 0; j < nr; j++) {
+	     buf[j] = 0.;
+	     b_column[j] = 0.;
+	   }
+	   //b_column = b(:,i)
+	   col_start_idx = b.cidx(i);
+	   col_end_idx = b.cidx(i+1);
+	   for (octave_idx_type j = col_start_idx; j < col_end_idx; j++) {
+	       b_column[b.ridx(j)] = b.data(j);
+	   }
+	   //buf = R'\(E_'* b(:,i))
+	   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	   CXSPARSE_DNAME (_ipvec) (E_,b_column,buf,nr);
+	   CXSPARSE_DNAME (_utsolve) (&R2,buf);
+	   END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	   //copy buf into S
+	   nz_counter = 0;
+	   for(octave_idx_type j = 0;j<nr;j++) if(buf[j] != 0.0) nz_counter++;
+	   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	   S = static_cast<cholmod_sparse*>(CHOLMOD_NAME(allocate_sparse)(nc,1,nz_counter,0,1,0,CHOLMOD_REAL,cc1));
+	   END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	   static_cast<octave_idx_type*>(S -> p)[1] = nz_counter;
+	   nz_counter = 0;
+	   for(volatile octave_idx_type j = 0;j<nr;j++) if(buf[j] != 0.0) {
+	       static_cast<octave_idx_type*>(S -> i)[nz_counter] = j;
+	       static_cast<double*>(S -> x)[nz_counter] = buf[j];
+	       nz_counter++;
+	   };
+	   //calculate Q * S
+	   BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	   QS = SuiteSparseQR_qmult<double> (SPQR_QX, H, Htau, HPinv,S,
+				      cc1);
+	   END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+	   //fill x(:,i). Resize x if necessary.
+	   col_start_idx = static_cast<octave_idx_type*>(QS -> p)[0];
+	   col_end_idx = static_cast<octave_idx_type*>(QS -> p)[1];
+	   for (volatile octave_idx_type k = col_start_idx; k < col_end_idx; k++)
+	    {
+	      octave_idx_type j = static_cast<octave_idx_type*>(QS -> i)[k];
+	      double tmp = static_cast<double*>(QS -> x)[k];
+	      if (ii == x_nz)
+		{
+		  // Resize the sparse matrix
+		  octave_idx_type sz = x_nz * (b_nc - i) / b_nc;
+		  sz = (sz > 10 ? sz : 10) + x_nz;
+		  x.change_capacity (sz);
+		  x_nz = sz;
+		}
+	      x.xdata (ii) = tmp;
+	      x.xridx (ii++) = j;
+	    }
+	    x.xcidx (i+1) = ii;
+	    CHOLMOD_NAME(free_sparse)(&QS,cc1);
+	    CHOLMOD_NAME(free_sparse)(&S,cc1);
+	 }
+	 CHOLMOD_NAME(finish)(cc1);
+	 info = 0;
+	 return x;
+
+#elif defined (HAVE_CXSPARSE)
 
       // These are swapped because the original matrix was transposed in
       // sparse_qr<SparseMatrix>::solve.
@@ -758,10 +1235,8 @@ namespace octave
     sparse_qr<SparseMatrix>::sparse_qr_rep::tall_solve<MArray<Complex>, ComplexMatrix>
     (const MArray<Complex>& b, octave_idx_type& info) const
     {
+#if defined (HAVE_CXSPARSE)
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
-
       octave_idx_type nr = nrows;
       octave_idx_type nc = ncols;
 
@@ -849,10 +1324,8 @@ namespace octave
     sparse_qr<SparseMatrix>::sparse_qr_rep::wide_solve<MArray<Complex>, ComplexMatrix>
     (const MArray<Complex>& b, octave_idx_type& info) const
     {
+#if defined (HAVE_CXSPARSE)
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
-
       // These are swapped because the original matrix was transposed in
       // sparse_qr<SparseMatrix>::solve.
 
@@ -942,17 +1415,38 @@ namespace octave
     }
 
     // Complex-valued matrices.
-
+#if defined (HAVE_SPQR)
+    template <>
+    sparse_qr<SparseComplexMatrix>::sparse_qr_rep::sparse_qr_rep
+    (const SparseComplexMatrix& a, int order)
+    : count (1), nrows (a.rows ()), ncols (a.columns ()),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 = 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)
     template <>
     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_CXSPARSE)
-
       CXSPARSE_ZNAME () A;
 
       A.nzmax = a.nnz ();
@@ -977,6 +1471,8 @@ namespace octave
         (*current_liboctave_error_handler)
           ("sparse_qr: sparse matrix QR factorization filled");
 
+      count = 1;
+
 #else
 
       octave_unused_parameter (order);
@@ -990,7 +1486,14 @@ namespace octave
     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
@@ -1038,7 +1541,25 @@ namespace octave
     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?
 
@@ -1080,7 +1601,40 @@ namespace octave
     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");
+      else if (b_nc <= 0 || b_nr <= 0)
+	(*current_liboctave_error_handler) (
+	  "matrix dimension with negative or zero size");
+      else
+	{
+	  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,7 +1697,34 @@ namespace octave
     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);
@@ -1448,10 +2029,8 @@ namespace octave
                                                               ComplexMatrix>
       (const MArray<double>& b, octave_idx_type& info) const
     {
+#if defined (HAVE_CXSPARSE)
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
-
       octave_idx_type nr = nrows;
       octave_idx_type nc = ncols;
 
@@ -1515,10 +2094,8 @@ namespace octave
                                                               ComplexMatrix>
       (const MArray<double>& b, octave_idx_type& info) const
     {
+#if defined (HAVE_CXSPARSE)
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
-
       // These are swapped because the original matrix was transposed in
       // sparse_qr<SparseComplexMatrix>::solve.
 
@@ -1590,10 +2167,8 @@ namespace octave
                                                               SparseComplexMatrix>
       (const SparseMatrix& b, octave_idx_type& info) const
     {
+#if defined (HAVE_CXSPARSE)
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
-
       octave_idx_type nr = nrows;
       octave_idx_type nc = ncols;
 
@@ -1686,10 +2261,8 @@ namespace octave
                                                               SparseComplexMatrix>
       (const SparseMatrix& b, octave_idx_type& info) const
     {
+#if defined (HAVE_CXSPARSE)
       info = -1;
-
-#if defined (HAVE_CXSPARSE)
-
       // These are swapped because the original matrix was transposed in
       // sparse_qr<SparseComplexMatrix>::solve.
 
@@ -2186,6 +2759,13 @@ namespace octave
     {
       return rep->P ();
     }
+    
+    template <typename SPARSE_T>
+    ColumnVector
+    sparse_qr<SPARSE_T>::E (void) const
+    {
+      return rep->E();
+    }
 
     template <typename SPARSE_T>
     SPARSE_T
@@ -2226,7 +2806,11 @@ namespace octave
     cxsparse_defaults<SparseMatrix>
     {
     public:
-      enum { order = 3 };
+#if defined (HAVE_SPQR)
+      enum { order = SPQR_ORDERING_DEFAULT };
+#elif defined (HAVE_CXSPARSE)
+      enum {order = 3}
+#endif
     };
 
     template <>
@@ -2234,15 +2818,42 @@ namespace octave
     cxsparse_defaults<SparseComplexMatrix>
     {
     public:
-      enum { order = 2 };
+#if defined (HAVE_SPQR)
+      enum { order = SPQR_ORDERING_DEFAULT };
+#elif defined (HAVE_CXSPARSE)
+      enum {order = 2}
+#endif
     };
 
-    template <typename SPARSE_T>
-    template <typename RHS_T, typename RET_T>
+    template<typename SPARSE_T>
+    template<typename RHS_T, typename RET_T>
     RET_T
-    sparse_qr<SPARSE_T>::solve (const SPARSE_T& a, const RHS_T& b,
-                                octave_idx_type& info)
-    {
+    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 ();
@@ -2275,8 +2886,302 @@ namespace octave
 
           return q.ok () ? q.wide_solve<RHS_T, RET_T> (b, info) : RET_T ();
         }
+#else
+
+#endif
     }
 
+#if defined (HAVE_SPQR)
+    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);
+#endif
+
+    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;
+    }
+
+
     template <typename SPARSE_T>
     template <typename RHS_T, typename RET_T>
     RET_T
diff -r b470734a025e -r c91ea5884edd liboctave/numeric/sparse-qr.h
--- a/liboctave/numeric/sparse-qr.h	Sat Oct 05 14:04:36 2019 -0700
+++ b/liboctave/numeric/sparse-qr.h	Sun Oct 06 18:44:20 2019 +0200
@@ -23,6 +23,9 @@ along with Octave; see the file COPYING.
 
 #if ! defined (octave_sparse_qr_h)
 #define octave_sparse_qr_h 1
+#if ! defined (HAVE_SPQR)
+#define HAVE_SPQR 1
+#endif
 
 #include "octave-config.h"
 
@@ -51,9 +54,13 @@ namespace octave
     public:
 
       sparse_qr (void);
-
+#if (HAVE_SPQR)
+      //order = 7 selects SPQR default ordering
+      sparse_qr (const SPARSE_T& a, int order = 7);
+#elif(HAVE_CXSPARSE)
       sparse_qr (const SPARSE_T& a, int order = 0);
-
+#else
+#endif
       sparse_qr (const sparse_qr& a);
 
       ~sparse_qr (void);
@@ -61,6 +68,8 @@ namespace octave
       sparse_qr& operator = (const sparse_qr& a);
 
       bool ok (void) const;
+      
+      ColumnVector E (void) const; 
 
       SPARSE_T V (void) const;
 
@@ -83,6 +92,11 @@ namespace octave
 
     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 -r b470734a025e -r c91ea5884edd liboctave/util/oct-sparse.cc
--- a/liboctave/util/oct-sparse.cc	Sat Oct 05 14:04:36 2019 -0700
+++ b/liboctave/util/oct-sparse.cc	Sun Oct 06 18:44:20 2019 +0200
@@ -26,6 +26,8 @@ along with Octave; see the file COPYING.
 
 #include "lo-error.h"
 #include "oct-sparse.h"
+#include "dSparse.h"
+#include "CSparse.h"
 
 #if (defined (HAVE_AMD) || defined (HAVE_CCOLAMD)               \
      || defined (HAVE_CHOLMOD) || defined (HAVE_COLAMD)         \
@@ -75,6 +77,180 @@ namespace octave
 
     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)
+      {
+	switch (cc->status)
+	  {
+	  case CHOLMOD_OUT_OF_MEMORY:
+	    (*current_liboctave_error_handler) (
+		"sparse_qr: sparse matrix QR factorization failed - out of memory");
+	    break;
+	  case CHOLMOD_TOO_LARGE:
+	    (*current_liboctave_error_handler) (
+		"sparse_qr: sparse matrix QR factorization failed - integer overflow occurred");
+	    break;
+	  default:
+	    (*current_liboctave_error_handler) (
+		"sparse_qr: sparse matrix QR factorization failed");
+	    break;
+	  }
+      }
+    //FIXME: Free memory?
+    //FIXME: Can cc-status > 0 (CHOLMOD_NOT_POSDEF, CHOLMOD_DSMALL) occur?
+  }
 }
 
 #endif
diff -r b470734a025e -r c91ea5884edd liboctave/util/oct-sparse.h
--- a/liboctave/util/oct-sparse.h	Sat Oct 05 14:04:36 2019 -0700
+++ b/liboctave/util/oct-sparse.h	Sun Oct 06 18:44:20 2019 +0200
@@ -22,7 +22,15 @@ along with Octave; see the file COPYING.
 
 #if ! defined (octave_oct_sparse_h)
 #define octave_oct_sparse_h 1
+#if ! defined (HAVE_SUITESPARSE_SPQR_H)
+#define HAVE_SUITESPARSE_SPQR_H 1
+#endif
+#if ! defined (HAVE_SPQR)
+#define HAVE_SPQR 1
+#endif
 
+#include "dSparse.h"
+#include "CSparse.h"
 #include "octave-config.h"
 
 #if defined (HAVE_SUITESPARSE_AMD_H)
@@ -85,6 +93,9 @@ along with Octave; see the file COPYING.
 #  include <umfpack.h>
 #endif
 
+#if defined (HAVE_SUITESPARSE_SPQR_H)
+# include <suitesparse/SuiteSparseQR.hpp>
+#endif
 // Cope with new SuiteSparse versions
 
 #if defined (SUITESPARSE_VERSION)
@@ -179,6 +190,59 @@ namespace octave
 
   extern const octave_idx_type*
   to_octave_idx_type_ptr (const suitesparse_integer *i);
+  /**
+  * Converts a real sparse octave matrix to a real
+  * sparse cholmod matrix. Returns a "shallow" copy of a.
+  */
+  extern const cholmod_sparse
+  ros2rcs (const SparseMatrix& a);
+  /**
+  * Converts a real sparse cholmod matrix to a real
+  * sparse octave matrix. Returns a "shallow" copy of a.
+  */
+  extern SparseMatrix
+  rcs2ros (const cholmod_sparse *y);
+  /**
+  * Converts a real dense octave matrix a to a real
+  * dense cholmod matrix. Returns a "shallow" copy of a.
+  */
+  extern const cholmod_dense
+  rod2rcd (const MArray<double>& a);
+  /**
+  * Converts a complex dense octave matrix to a complex
+  * dense cholmod matrix. Returns a "shallow" copy of a.
+  */
+  extern const cholmod_dense
+  cod2ccd (const ComplexMatrix &a);
+  /**
+  * Converts a complex sparse octave matrix to a complex
+  * sparse cholmod matrix. Returns a "shallow" copy of a.
+  */
+  extern const cholmod_sparse
+  cos2ccs (const SparseComplexMatrix &a);
+  /**
+  * Converts a complex sparse cholmod matrix to a complex
+  * sparse octave matrix. Returns a "deep" copy of a.
+  */
+  extern SparseComplexMatrix
+  ccs2cos (const cholmod_sparse *a);
+  /**
+  * Converts a real sparse octave matrix to a 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);
+  /**
+  * Converts a real dense octave matrix to a 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
