# HG changeset patch
# User address@hidden
# Date 1563570038 -7200
#      Fri Jul 19 23:00:38 2019 +0200
# Node ID 254f96a055e185f57e109fdc90466d1322033584
# Parent  596312d4f25d5dc500fcc482811e01cf35ea4b19
test patch for testing sparse-qr methods Q and R with spqr instead of cxsparse

diff -r 596312d4f25d -r 254f96a055e1 libinterp/dldfcn/qr.cc
--- a/libinterp/dldfcn/qr.cc	Fri Jul 19 10:02:20 2019 -0400
+++ b/libinterp/dldfcn/qr.cc	Fri Jul 19 23:00:38 2019 +0200
@@ -330,9 +330,8 @@ orthogonal basis of @code{span (A)}.
 
   if (arg.issparse ())
     {
-      if (nargout > 2)
-        error ("qr: Permutation output is not supported for sparse input");
-
+      if (nargout > 3)
+        error ("qr: too many output arguments");
       if (is_cmplx)
         {
           octave::math::sparse_qr<SparseComplexMatrix> q (arg.sparse_complex_matrix_value ());
@@ -355,7 +354,6 @@ orthogonal basis of @code{span (A)}.
       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));
@@ -365,10 +363,12 @@ orthogonal basis of @code{span (A)}.
                          "x%" OCTAVE_IDX_TYPE_FORMAT,
                          arg.rows (), arg.columns ());
             }
+          else 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));
+              retval = ovl (q.Q (), q.R (economy));
+          else 
+            retval = ovl (q.R (economy)); 
         }
     }
   else
diff -r 596312d4f25d -r 254f96a055e1 liboctave/numeric/sparse-qr.cc
--- a/liboctave/numeric/sparse-qr.cc	Fri Jul 19 10:02:20 2019 -0400
+++ b/liboctave/numeric/sparse-qr.cc	Fri Jul 19 23:00:38 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"
@@ -37,6 +36,8 @@ along with Octave; see the file COPYING.
 #include "oct-sparse.h"
 #include "quit.h"
 #include "sparse-qr.h"
+#include "suitesparse/SuiteSparseQR.hpp"
+
 
 namespace octave
 {
@@ -53,7 +54,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
@@ -80,9 +85,8 @@ namespace octave
     class sparse_qr<SPARSE_T>::sparse_qr_rep
     {
     public:
-
       sparse_qr_rep (const SPARSE_T& a, int order);
-
+      
       // No copying!
 
       sparse_qr_rep (const sparse_qr_rep&) = delete;
@@ -105,6 +109,8 @@ namespace octave
       ColumnVector Pinv (void) const;
 
       ColumnVector P (void) const;
+      
+      ColumnVector E (void) const;
 
       SPARSE_T R (bool econ) const;
 
@@ -118,10 +124,21 @@ namespace octave
 
       octave_idx_type nrows;
       octave_idx_type ncols;
-
+#if defined (HAVE_SPQR)
+      cholmod_common* cc; 
+      cholmod_sparse* Rfactor; 
+      suitesparse_integer* Evector;
+      cholmod_sparse *Hfactor; 
+      suitesparse_integer *HPinv; 
+      cholmod_dense *Htau;
       typename cxsparse_types<SPARSE_T>::symbolic_type *S;
       typename cxsparse_types<SPARSE_T>::numeric_type *N;
-
+#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
       tall_solve (const RHS_T& b, octave_idx_type& info) const;
@@ -130,7 +147,7 @@ namespace octave
       RET_T
       wide_solve (const RHS_T& b, octave_idx_type& info) const;
     };
-
+    
     template <typename SPARSE_T>
     ColumnVector
     sparse_qr<SPARSE_T>::sparse_qr_rep::Pinv (void) const
@@ -155,7 +172,12 @@ namespace octave
     ColumnVector
     sparse_qr<SPARSE_T>::sparse_qr_rep::P (void) const
     {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+    ColumnVector ret(nrows);
+    for (octave_idx_type i = 0; i < nrows; i++)
+    	ret.xelem(HPinv[i]) = i;
+    return ret;
+#elif defined (HAVE_CXSPARSE)
 
       ColumnVector ret (N->L->m);
 
@@ -170,7 +192,26 @@ 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 (int i = 0; i < ncols; i++)
+    	ret(i) = static_cast<octave_idx_type> (Evector[i] + 1);
+    return ret;
+#elif defined (HAVE_CXSPARSE)
+    return ColumnVector ();
+#else
+      return ColumnVector ();
+#endif
+    }
+    
+    
+    
+    
     // Specializations.
 
     // Real-valued matrices.
@@ -179,11 +220,47 @@ namespace octave
     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)
+    
+#elif (HAVE_CXSPARSE)
       , S (nullptr), N (nullptr)
+#else
 #endif
     {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+      cholmod_common Common;
+      cc = &Common;
+      Evector = NULL; 
+      HPinv = NULL; 
+      cholmod_l_start(cc);   
+      cholmod_sparse tmp_A;
+      cholmod_sparse * A = &tmp_A;  // = cholmod_l_allocate_sparse(a.rows(),a.cols(),a.nnz(),0,1,0,CHOLMOD_REAL,cc); 
+      A -> ncol = a.cols();
+      A -> nrow = a.rows();
+      A -> itype = CHOLMOD_LONG;
+      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->p = NULL ;
+      A->i = NULL ;
+      A->x = NULL ;
+      A->z = NULL ;
+      //TODO: Check casts between SuiteSparse_long, suitesparse_integer and octave_idx_type.
+      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 ());
+      BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      SuiteSparseQR<double>(SPQR_ORDERING_DEFAULT,SPQR_DEFAULT_TOL,A -> nrow,A,&Rfactor,&Evector,&Hfactor,&HPinv,&Htau,cc);
+      //TODO: Check cc->status and call octave error handler if necessary.
+      END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+      
+#elif defined (HAVE_CXSPARSE)
 
       CXSPARSE_DNAME () A;
 
@@ -223,9 +300,15 @@ namespace octave
     template <>
     sparse_qr<SparseMatrix>::sparse_qr_rep::~sparse_qr_rep (void)
     {
-#if defined (HAVE_CXSPARSE)
+#if defined (HAVE_SPQR)
+     cholmod_l_free_sparse(&Rfactor,cc);
+     cholmod_l_free_sparse(&Hfactor,cc);
+     cholmod_l_free_dense(&Htau,cc);
+     cholmod_l_finish(cc); 
+#elif defined (HAVE_CXSPARSE)
       CXSPARSE_DNAME (_sfree) (S);
       CXSPARSE_DNAME (_nfree) (N);
+#else
 #endif
     }
 
@@ -272,7 +355,21 @@ namespace octave
     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> (Rfactor -> nrow); 
+      octave_idx_type nc = static_cast<octave_idx_type> (Rfactor -> ncol);
+      octave_idx_type nz = static_cast<octave_idx_type> (Rfactor -> nzmax);
+      SparseMatrix ret ((econ ? (nc > nr ? nr : nc) : nr), nc, nz);
+      for (octave_idx_type j = 0; j < nc+1; j++)
+        ret.xcidx (j) = (static_cast<octave_idx_type *>(Rfactor -> p))[j];
+
+      for (octave_idx_type j = 0; j < nz; j++)
+        {
+          ret.xridx (j) = (static_cast<octave_idx_type *>(Rfactor -> i))[j];
+          ret.xdata (j) = (static_cast<double *>(Rfactor -> x))[j];
+      }
+      return ret;
+#elif defined (HAVE_CXSPARSE)
 
       // Drop zeros from R and sort
       // FIXME: Is the double transpose to sort necessary?
@@ -379,8 +476,34 @@ namespace octave
     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; 
+     //TODO: Use cc instead of cc1.
+     cholmod_common * cc1 = &Common1; 
+     cholmod_l_start(cc1); 
+     cholmod_dense* cholmod_ret; 
+     cholmod_dense* X = static_cast<cholmod_dense*>(cholmod_l_allocate_dense(nrows,nrows,nrows,CHOLMOD_REAL,cc1)); 
+     for (int i = 0; i < nrows * nrows; i++) (static_cast<double *>(X -> x))[i] = 0.0; 
+     for (int i = 0; i < nrows; i++)(static_cast<double *>(X -> x))[i * nrows + i] = 1.0;
+     BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+     cholmod_ret = SuiteSparseQR_qmult<double>(SPQR_QX, Hfactor, Htau, HPinv, X, cc1); 
+     //TODO: Check cc1->status.
+     END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE;
+     double * cholmod_x = static_cast<double*>(cholmod_ret -> x); 
+     //TODO: Does cholmod store matrices in column major format? Avoid this copy operation.
+     for (int j = 0; j < nrows; j++) {
+         for (int i = 0;i < nrows; i++) {
+             ret(i,j) = cholmod_x[j * nrows + i];
+         }
+     }
+     cholmod_l_free_dense(&cholmod_ret,cc1);
+     cholmod_l_free_dense(&X,cc1);
+     cholmod_l_finish(cc1);
+     return ret;
+#elif defined (HAVE_CXSPARSE)
+        
       octave_idx_type nc = N->L->n;
       octave_idx_type nr = nrows;
       Matrix ret (nr, nr);
@@ -438,7 +561,6 @@ namespace octave
 
 #endif
     }
-
     template <>
     template <>
     Matrix
@@ -2190,6 +2312,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
diff -r 596312d4f25d -r 254f96a055e1 liboctave/numeric/sparse-qr.h
--- a/liboctave/numeric/sparse-qr.h	Fri Jul 19 10:02:20 2019 -0400
+++ b/liboctave/numeric/sparse-qr.h	Fri Jul 19 23:00:38 2019 +0200
@@ -23,6 +23,7 @@ along with Octave; see the file COPYING.
 
 #if ! defined (octave_sparse_qr_h)
 #define octave_sparse_qr_h 1
+#define HAVE_SPQR 1
 
 #include "octave-config.h"
 
@@ -61,6 +62,8 @@ namespace octave
       sparse_qr& operator = (const sparse_qr& a);
 
       bool ok (void) const;
+      
+      ColumnVector E (void) const; 
 
       SPARSE_T V (void) const;