Add more useful functions to SPQR interface

2 files changed, 51 insertions, 7 deletions

diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
index 35dba2e68..2647d22f0 100644
--- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
+++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
@@ -36,7 +36,7 @@ namespace Eigen {
  * \class SPQR
  * \brief Sparse QR factorization based on SuiteSparseQR library
  * 
- * This class is used to perform a multithreaded and multifrontal QR decomposition 
+ * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition 
  * of sparse matrices. The result is then used to solve linear leasts_square systems.
  * Clearly, a QR factorization is returned such that A*P = Q*R where :
  * 
@@ -61,6 +61,7 @@ class SPQR
     typedef typename _MatrixType::RealScalar RealScalar;
     typedef UF_long Index ; 
     typedef SparseMatrix<Scalar, _MatrixType::Flags, Index> MatrixType;
+    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
   public:
     SPQR() 
     : m_ordering(SPQR_ORDERING_DEFAULT),
@@ -115,8 +116,8 @@ class SPQR
         // Solves with the triangular matrix R
       Dest y; 
       y = this->matrixQR().template triangularView<Upper>().solve(dest.derived());
-      // Apply the column permutation //TODO Check the terminology behind the permutation
-      for (int j = 0; j < y.size(); j++)  dest(m_E[j]) = y(j); 
+      // Apply the column permutation 
+      dest = colsPermutation() * y;
       
       m_info = Success;
     }
@@ -133,13 +134,29 @@ class SPQR
       return SPQRMatrixQReturnType<SPQR>(*this);
     }
     /// Get the permutation that was applied to columns of A
-    Index *colsPermutation() { return m_E; }
-    
+    PermutationType colsPermutation() const
+    { 
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+      Index n = m_cR->ncol;
+      PermutationType colsPerm(n);
+      for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j];
+      return colsPerm; 
+      
+    }
+    /**
+     * Gets the rank of the matrix. 
+     * It should be equal to matrixQR().cols if the matrix is full-rank
+     */
+    Index rank() const
+    {
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+      return m_cc.SPQR_istat[4];
+    }
     /// Set the fill-reducing ordering method to be used
     void setOrdering(int ord) { m_ordering = ord;}
     /// Set the tolerance tol to treat columns with 2-norm < =tol as zero
-    void setTolerance(RealScalar tol) { m_tolerance = tol; }
-      
+    void setThreshold(RealScalar tol) { m_tolerance = tol; }
+    
     /// Return a pointer to SPQR workspace 
     cholmod_common *cc() const { return &m_cc; }
     cholmod_sparse * H() const { return m_H; }
diff --git a/cmake/FindSPQR.cmake b/cmake/FindSPQR.cmake
new file mode 100644
index 000000000..55c4d9f6e
--- /dev/null
+++ b/cmake/FindSPQR.cmake
@@ -0,0 +1,27 @@
+# SPQR lib usually requires linking to a blas and lapack library.
+# It is up to the user of this module to find a BLAS and link to it.
+
+# SPQR lib requires Cholmod, colamd and amd as well. 
+# FindCholmod.cmake can be used to find those packages before finding spqr
+
+if (SPQR_INCLUDES AND SPQR_LIBRARIES)
+  set(SPQR_FIND_QUIETLY TRUE)
+endif (SPQR_INCLUDES AND SPQR_LIBRARIES)
+
+find_path(SPQR_INCLUDES
+  NAMES
+  SuiteSparseQR.hpp
+  PATHS
+  $ENV{SPQRDIR}
+  ${INCLUDE_INSTALL_DIR}
+  PATH_SUFFIXES
+  suitesparse
+  ufsparse
+)
+
+find_library(SPQR_LIBRARIES spqr $ENV{SPQRDIR} ${LIB_INSTALL_DIR})
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(SPQR DEFAULT_MSG SPQR_INCLUDES SPQR_LIBRARIES)
+
+mark_as_advanced(SPQR_INCLUDES SPQR_LIBRARIES)
\ No newline at end of file