Add a sparse QR factorization and update the elimination tree in SparseLU

10 files changed, 618 insertions, 36 deletions

diff --git a/Eigen/SparseQR b/Eigen/SparseQR
new file mode 100644
index 000000000..d2382666a
--- /dev/null
+++ b/Eigen/SparseQR
@@ -0,0 +1,28 @@
+#ifndef EIGEN_SPARSEQR_MODULE_H
+#define EIGEN_SPARSEQR_MODULE_H
+
+#include "SparseCore"
+#include "OrderingMethods"
+#include "src/Core/util/DisableStupidWarnings.h"
+
+/** \defgroup SparseQR_Module SparseQR module
+  * 
+  * 
+  * This module provides a simplicial version of the left-looking Sparse QR decomposition. 
+  * The columns of the input matrix should be reordered to limit the fill-in during the 
+  * decomposition. Built-in methods (COLAMD, AMD) or external  methods (METIS) can be used to this end.
+  * See \link OrderingMethods_Module OrderingMethods_Module \endlink for the list 
+  * of built-in and external ordering methods.
+  * 
+  * \code
+  * #include <Eigen/SparseQR>
+  * \endcode
+  * 
+  * 
+  */
+
+#include "src/SparseQR/SparseQR.h"
+
+#include "src/Core/util/ReenableStupidWarnings.h"
+
+#endif
\ No newline at end of file
diff --git a/Eigen/src/MetisSupport/MetisSupport.h b/Eigen/src/MetisSupport/MetisSupport.h
index a762d96f6..3a723b384 100644
--- a/Eigen/src/MetisSupport/MetisSupport.h
+++ b/Eigen/src/MetisSupport/MetisSupport.h
@@ -122,10 +122,9 @@ public:
     //NOTE:  If Ap is the permuted matrix then perm and iperm vectors are defined as follows 
     // Row (column) i of Ap is the perm(i) row(column) of A, and row (column) i of A is the iperm(i) row(column) of Ap
     
-    // To be consistent with the use of the permutation in SparseLU module, we thus keep the iperm vector 
      matperm.resize(m);
      for (int j = 0; j < m; j++)
-       matperm.indices()(j) = iperm(j); 
+       matperm.indices()(iperm(j)) = j;
    
   }
   
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 2f515eb41..6003237fe 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -170,7 +170,7 @@ class SparseLU
      /** \brief Reports whether previous computation was successful.
       *
       * \returns \c Success if computation was succesful,
-      *          \c NumericalIssue if the PaStiX reports a problem
+      *          \c NumericalIssue if the LU factorization reports a problem, zero diagonal for instance
       *          \c InvalidInput if the input matrix is invalid
       *
       * \sa iparm()          
@@ -320,15 +320,14 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   }
     
   // Compute the column elimination tree of the permuted matrix 
-  /*if (m_etree.size() == 0)  */m_etree.resize(m_mat.cols());
-  
-  SparseLUBase<Scalar,Index>::LU_sp_coletree(m_mat, m_etree); 
+  IndexVector firstRowElt;
+  internal::coletree(m_mat, m_etree,firstRowElt); 
      
   // In symmetric mode, do not do postorder here
   if (!m_symmetricmode) {
     IndexVector post, iwork; 
     // Post order etree
-    SparseLUBase<Scalar,Index>::LU_TreePostorder(m_mat.cols(), m_etree, post); 
+    internal::treePostorder(m_mat.cols(), m_etree, post); 
       
    
     // Renumber etree in postorder 
@@ -445,13 +444,12 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   // Work on one 'panel' at a time. A panel is one of the following :
   //  (a) a relaxed supernode at the bottom of the etree, or
   //  (b) panel_size contiguous columns, <panel_size> defined by the user
-  int jcol,kcol; 
+  int jcol; 
   IndexVector panel_histo(n);
-  Index nextu, nextlu, jsupno, fsupc, new_next;
   Index pivrow; // Pivotal row number in the original row matrix
   int nseg1; // Number of segments in U-column above panel row jcol
   int nseg; // Number of segments in each U-column 
-  int irep, icol; 
+  int irep; 
   int i, k, jj; 
   for (jcol = 0; jcol < n; )
   {
diff --git a/Eigen/src/SparseLU/SparseLUBase.h b/Eigen/src/SparseLU/SparseLUBase.h
index 20338a646..893728b5e 100644
--- a/Eigen/src/SparseLU/SparseLUBase.h
+++ b/Eigen/src/SparseLU/SparseLUBase.h
@@ -22,10 +22,6 @@ struct SparseLUBase
   typedef LU_GlobalLU_t<IndexVector, ScalarVector> GlobalLU_t; 
   typedef SparseMatrix<Scalar,ColMajor,Index> MatrixType; 
   
-  static int etree_find (int i, IndexVector& pp); 
-  static int LU_sp_coletree(const MatrixType& mat, IndexVector& parent);
-  static void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum);
-  static void LU_TreePostorder(int n, IndexVector& parent, IndexVector& post);
   template <typename VectorType>
   static int expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions);
   static int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,  GlobalLU_t& glu); 
diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h
index d3bc36ea4..e373525ad 100644
--- a/Eigen/src/SparseLU/SparseLU_Coletree.h
+++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -30,9 +30,11 @@
  */
 #ifndef SPARSELU_COLETREE_H
 #define SPARSELU_COLETREE_H
+
+namespace internal {
 /** Find the root of the tree/set containing the vertex i : Use Path halving */ 
-template< typename Scalar,typename Index>
-int SparseLUBase<Scalar,Index>::etree_find (int i, IndexVector& pp)
+template<typename IndexVector>
+int etree_find (int i, IndexVector& pp)
 {
   int p = pp(i); // Parent 
   int gp = pp(p); // Grand parent 
@@ -50,45 +52,53 @@ int SparseLUBase<Scalar,Index>::etree_find (int i, IndexVector& pp)
   * NOTE : The matrix is supposed to be in column-major format. 
   * 
   */
-template <typename Scalar, typename Index>
-int SparseLUBase<Scalar,Index>::LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
+template <typename MatrixType, typename IndexVector>
+int coletree(const MatrixType& mat, IndexVector& parent, IndexVector& firstRowElt)
 {
-  int nc = mat.cols(); // Number of columns 
-  int nr = mat.rows(); // Number of rows 
-  
+  typedef typename MatrixType::Index Index;
+  Index nc = mat.cols(); // Number of columns 
+  Index m = mat.rows();
   IndexVector root(nc); // root of subtree of etree 
   root.setZero();
   IndexVector pp(nc); // disjoint sets 
   pp.setZero(); // Initialize disjoint sets 
-  IndexVector firstcol(nr); // First nonzero column in each row 
-  
+  parent.resize(mat.cols());
   //Compute first nonzero column in each row 
   int row,col; 
-  firstcol.setConstant(nc);  //for (row = 0; row < nr; firstcol(row++) = nc); 
+  firstRowElt.resize(m);
+  firstRowElt.setConstant(nc);
+  firstRowElt.segment(0, nc).setLinSpaced(nc, 0, nc-1);
+  bool found_diag;
   for (col = 0; col < nc; col++)
   {
     for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
-    { // Is it necessary to browse the whole matrix, the lower part should do the job ??
+    { 
       row = it.row();
-      firstcol(row) = (std::min)(firstcol(row), col);
+      firstRowElt(row) = (std::min)(firstRowElt(row), col);
     }
   }
   /* Compute etree by Liu's algorithm for symmetric matrices,
-          except use (firstcol[r],c) in place of an edge (r,c) of A.
+          except use (firstRowElt[r],c) in place of an edge (r,c) of A.
     Thus each row clique in A'*A is replaced by a star
     centered at its first vertex, which has the same fill. */
   int rset, cset, rroot; 
   for (col = 0; col < nc; col++) 
   {
+    found_diag = false;
     pp(col) = col; 
     cset = col; 
     root(cset) = col; 
     parent(col) = nc; 
-    for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
+    /* The diagonal element is treated here even if it does not exist in the matrix
+     * hence the loop is executed once more */ 
+    for (typename MatrixType::InnerIterator it(mat, col); it||!found_diag; ++it)
     { //  A sequence of interleaved find and union is performed 
-      row = firstcol(it.row());
+      int i = col;
+      if(it) i = it.index();
+      if (i == col) found_diag = true;
+      row = firstRowElt(i);
       if (row >= col) continue; 
-      rset = etree_find(row, pp); // Find the name of the set containing row
+      rset = internal::etree_find(row, pp); // Find the name of the set containing row
       rroot = root(rset);
       if (rroot != col) 
       {
@@ -106,8 +116,8 @@ int SparseLUBase<Scalar,Index>::LU_sp_coletree(const MatrixType& mat, IndexVecto
   * Depth-first search from vertex n.  No recursion.
   * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
 */
-template <typename Scalar, typename Index>
-void SparseLUBase<Scalar,Index>::LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum)
+template <typename IndexVector>
+void nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum)
 {
   int current = n, first, next;
   while (postnum != n) 
@@ -152,8 +162,8 @@ void SparseLUBase<Scalar,Index>::LU_nr_etdfs (int n, IndexVector& parent, IndexV
   * \param parent Input tree
   * \param post postordered tree
   */
-template <typename Scalar, typename Index>
-void SparseLUBase<Scalar,Index>::LU_TreePostorder(int n, IndexVector& parent, IndexVector& post)
+template <typename IndexVector>
+void treePostorder(int n, IndexVector& parent, IndexVector& post)
 {
   IndexVector first_kid, next_kid; // Linked list of children 
   int postnum; 
@@ -174,7 +184,8 @@ void SparseLUBase<Scalar,Index>::LU_TreePostorder(int n, IndexVector& parent, In
   
   // Depth-first search from dummy root vertex #n
   postnum = 0; 
-  LU_nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
+  internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
 }
 
+} // internal
 #endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
index 69e1d4da9..7f6b3bab2 100644
--- a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
+++ b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
@@ -44,7 +44,7 @@ void SparseLUBase<Scalar,Index>::LU_heap_relax_snode (const int n, IndexVector&
   
   // The etree may not be postordered, but its heap ordered  
   IndexVector post;
-  LU_TreePostorder(n, et, post); // Post order etree
+  internal::treePostorder(n, et, post); // Post order etree
   IndexVector inv_post(n+1); 
   int i;
   for (i = 0; i < n+1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()???
diff --git a/Eigen/src/SparseQR/CMakeLists.txt b/Eigen/src/SparseQR/CMakeLists.txt
new file mode 100644
index 000000000..15009e6b5
--- /dev/null
+++ b/Eigen/src/SparseQR/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_SparseQRSupport_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_SparseQRSupport_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseQRSupport/ COMPONENT Devel
+  )
diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h
new file mode 100644
index 000000000..f1e5509dd
--- /dev/null
+++ b/Eigen/src/SparseQR/SparseQR.h
@@ -0,0 +1,481 @@
+#ifndef EIGEN_SPARSE_QR_H
+#define EIGEN_SPARSE_QR_H
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
+// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+
+namespace Eigen {
+#include "../SparseLU/SparseLU_Coletree.h"
+
+template<typename MatrixType, typename OrderingType> class SparseQR;
+template<typename SparseQRType> struct SparseQRMatrixQReturnType;
+template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType;
+template<typename SparseQRType, typename Derived> struct SparseQR_QProduct;
+namespace internal {
+  template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> >
+  {
+    typedef typename SparseQRType::MatrixType ReturnType;
+  };
+  template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> >
+  {
+    typedef typename SparseQRType::MatrixType ReturnType;
+  };
+  template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> >
+  {
+    typedef typename Derived::PlainObject ReturnType;
+  };
+} // End namespace internal
+
+/**
+ * \ingroup SparseQRSupport_Module
+ * \class SparseQR
+ * \brief Sparse left-looking QR factorization
+ * 
+ * This class is used to perform a left-looking 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 :
+ * 
+ * P is the column permutation. Use colsPermutation() to get it.
+ * 
+ * Q is the orthogonal matrix represented as Householder reflectors. 
+ * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
+ * You can then apply it to a vector.
+ * 
+ * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
+ * 
+ * NOTE This is not a rank-revealing QR decomposition.
+ * 
+ * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
+ * \tparam _OrderingType The fill-reducing ordering method. See \link OrderingMethods_Module 
+ *  OrderingMethods_Module \endlink for the list of built-in and external ordering methods.
+ * 
+ * 
+ */
+template<typename _MatrixType, typename _OrderingType>
+class SparseQR
+{
+  public:
+    typedef _MatrixType MatrixType;
+    typedef _OrderingType OrderingType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename MatrixType::Index Index;
+    typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType;
+    typedef Matrix<Index, Dynamic, 1> IndexVector;
+    typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
+    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+  public:
+    SparseQR () : m_isInitialized(false),m_analysisIsok(false)
+    { }
+    
+    SparseQR(const MatrixType& mat) : m_isInitialized(false),m_analysisIsok(false)
+    {
+      compute(mat);
+    }
+    void compute(/*const*/ MatrixType& mat)
+    {
+      analyzePattern(mat);
+      factorize(mat);
+    }
+    void analyzePattern(const MatrixType& mat);
+    void factorize(/*const*/ MatrixType& mat);
+    
+    /** 
+     * Get the number of rows of the triangular matrix. 
+     */
+    inline Index rows() const { return m_R.rows(); }
+    
+    /** 
+     * Get the number of columns of the triangular matrix. 
+     */
+    inline Index cols() const { return m_R.cols();}
+    
+    /**
+     * This is the number of rows in the input matrix and the Q matrix
+     */
+    inline Index rowsQ() const {return m_pmat.rows(); }
+    
+    /// Get the sparse triangular matrix R. It is a sparse matrix
+    MatrixType matrixQR() const
+    {
+      return m_R;
+    }
+    /// Get an expression of the matrix Q
+    SparseQRMatrixQReturnType<SparseQR> matrixQ() const 
+    {
+      return SparseQRMatrixQReturnType<SparseQR>(*this);
+    }
+    
+    /// Get the permutation that was applied to columns of A
+    PermutationType colsPermutation() const
+    { 
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+      return m_perm_c;
+    }
+    template<typename Rhs, typename Dest>
+    bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const
+    {
+      eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
+      eigen_assert(this->rowsQ() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
+      Index rank = this->matrixQR().cols();
+      // Compute Q^T * b;
+      dest = this->matrixQ().transpose() * B;      
+      // Solve with the triangular matrix R
+      Dest y;
+      y = this->matrixQR().template triangularView<Upper>().solve(dest.derived().topRows(rank));
+      // Apply the column permutation
+      if (m_perm_c.size())
+        dest.topRows(rank) =  colsPermutation().inverse() * y;
+      else 
+        dest = y;
+      m_info = Success;
+      return true;
+    }
+    /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
+      *
+      * \sa compute()
+      */
+    template<typename Rhs>
+    inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 
+    {
+      eigen_assert(m_isInitialized && "The factorization should be called first, use compute()");
+      eigen_assert(this->rowsQ() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix");
+      return internal::solve_retval<SparseQR, Rhs>(*this, B.derived());
+    }
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful,
+      *          \c NumericalIssue if the QR factorization reports a problem
+      *          \c InvalidInput if the input matrix is invalid
+      *
+      * \sa iparm()          
+      */
+    ComputationInfo info() const
+    {
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+      return m_info;
+    }
+  protected:
+    bool m_isInitialized;
+    bool m_analysisIsok;
+    bool m_factorizationIsok;
+    mutable ComputationInfo m_info;
+    QRMatrixType m_pmat; // Temporary matrix
+    QRMatrixType m_R; // The triangular factor matrix
+    QRMatrixType m_Q; // The orthogonal reflectors
+    ScalarVector m_hcoeffs; // The Householder coefficients
+    PermutationType m_perm_c; // Column  permutation 
+    PermutationType m_perm_r; // Column permutation 
+    IndexVector m_etree; // Column elimination tree
+    IndexVector m_firstRowElt; // First element in each row
+    IndexVector m_found_diag_elem; // Existence of diagonal elements
+    template <typename, typename > friend struct SparseQR_QProduct;
+    
+};
+/**
+ * \brief Preprocessing step of a QR factorization 
+ * 
+ * In this step, the fill-reducing permutation is computed and applied to the columns of A
+ * and the column elimination tree is computed as well.
+ * \note In this step it is assumed that there is no empty row in the matrix \p mat
+ */
+template <typename MatrixType, typename OrderingType>
+void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
+{
+  // Compute the column fill reducing ordering
+  OrderingType ord; 
+  ord(mat, m_perm_c); 
+  Index n = mat.cols();
+  Index m = mat.rows();
+  // Permute the input matrix... only the column pointers are permuted
+  m_pmat = mat;
+  m_pmat.uncompress();
+  for (int i = 0; i < n; i++)
+  {
+    Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i;
+//     m_pmat.outerIndexPtr()[i] = mat.outerIndexPtr()[p]; 
+//     m_pmat.innerNonZeroPtr()[i] = mat.outerIndexPtr()[p+1] - mat.outerIndexPtr()[p]; 
+    m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; 
+    m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 
+  }
+  // Compute the column elimination tree of the permuted matrix
+  internal::coletree(m_pmat, m_etree, m_firstRowElt/*, m_found_diag_elem*/);
+  
+  m_R.resize(n, n);
+  m_Q.resize(m, m);
+  // Allocate space for nonzero elements : rough estimation
+  m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
+  m_Q.reserve(2*mat.nonZeros());
+  m_hcoeffs.resize(n);
+  m_analysisIsok = true;
+}
+
+/**
+ * \brief Perform the numerical QR factorization of the input matrix
+ * 
+ * \param mat The sparse column-major matrix
+ */
+template <typename MatrixType, typename OrderingType>
+void SparseQR<MatrixType,OrderingType>::factorize(MatrixType& mat)
+{
+  eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
+  Index m = mat.rows();
+  Index n = mat.cols();
+  IndexVector mark(m); mark.setConstant(-1);// Record the visited nodes
+  IndexVector Ridx(n),Qidx(m); // Store temporarily the row indexes for the current column of R and Q
+  Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
+  Index pcol;
+  ScalarVector tval(m); tval.setZero(); // Temporary vector
+  IndexVector iperm(m);
+  bool found_diag;
+  if (m_perm_c.size())
+    for(int i = 0; i < m; i++) iperm(m_perm_c.indices()(i)) = i;
+  else
+    iperm.setLinSpaced(m, 0, m-1);
+      
+  // Left looking QR factorization : Compute a column of R and Q at a time
+  for (Index col = 0; col < n; col++)
+  {
+    m_R.startVec(col);
+    m_Q.startVec(col);
+    mark(col) = col;
+    Qidx(0) = col; 
+    nzcolR = 0; nzcolQ = 1;
+    pcol = iperm(col);
+    found_diag = false;
+    // Find the nonzero locations of the column k of R, 
+    // i.e All the nodes (with indexes lower than k) reachable through the col etree rooted at node k
+    for (typename MatrixType::InnerIterator itp(mat, pcol); itp || !found_diag; ++itp)
+    {
+      Index curIdx = col;
+      if (itp) curIdx = itp.row();
+      if(curIdx == col) found_diag = true;
+      // Get the nonzeros indexes  of the current column of R
+      Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here 
+      if (st < 0 )
+      {
+        std::cerr << " Empty row found during Numerical factorization ... Abort \n";
+        m_info = NumericalIssue;
+        return;
+      }
+      // Traverse the etree 
+      Index bi = nzcolR;
+      for (; mark(st) != col; st = m_etree(st))
+      {
+        Ridx(nzcolR) = st; // Add this row to the list 
+        mark(st) = col; // Mark this row as visited
+        nzcolR++;
+      }
+      // Reverse the list to get the topological ordering
+      Index nt = nzcolR-bi;
+      for(int i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1));
+       
+      // Copy the current row value of mat
+      if (itp) tval(curIdx) = itp.value();
+      else tval(curIdx) = Scalar(0.);
+      
+      // Compute the pattern of Q(:,k)
+      if (curIdx > col && mark(curIdx) < col) 
+      {
+        Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q
+        mark(curIdx) = col; // And mark it as visited
+        nzcolQ++;
+      }
+    }
+    
+    //Browse all the indexes of R(:,col) in reverse order
+    for (Index i = nzcolR-1; i >= 0; i--)
+    {
+      Index curIdx = Ridx(i);
+      // Apply the <curIdx> householder vector  to tval
+      Scalar tdot(0.);
+      //First compute q'*tval
+      for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
+      {
+        tdot += internal::conj(itq.value()) * tval(itq.row());
+      }
+      tdot *= m_hcoeffs(curIdx);
+      // Then compute tval = tval - q*tau
+      for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
+      {
+        tval(itq.row()) -= itq.value() * tdot;
+      }
+      //With the topological ordering, updates for curIdx are fully done at this point
+      m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx);
+      tval(curIdx) = Scalar(0.);
+      
+      //Detect fill-in for the current column of Q
+      if(m_etree(curIdx) == col)
+      {
+        for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq)
+        {
+          Index iQ = itq.row();
+          if (mark(iQ) < col)
+          {
+            Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q
+            mark(iQ) = col; //And mark it as visited
+          }
+        }
+      }
+    } // End update current column of R
+    
+    // Record the current (unscaled) column of V.
+    for (Index itq = 0; itq < nzcolQ; ++itq)
+    {
+      Index iQ = Qidx(itq);      
+      m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ);
+      tval(iQ) = Scalar(0.);
+    }
+    // Compute the new Householder reflection
+    RealScalar sqrNorm =0.;
+    Scalar tau; RealScalar beta;
+    typename QRMatrixType::InnerIterator itq(m_Q, col);
+    Scalar c0 = (itq) ? itq.value() : Scalar(0.);
+    //First, the squared norm of Q((col+1):m, col)
+    if(itq) ++itq;
+    for (; itq; ++itq)
+    {
+      sqrNorm += internal::abs2(itq.value());
+    }
+    if(sqrNorm == RealScalar(0) && internal::imag(c0) == RealScalar(0))
+    {
+      tau = RealScalar(0);
+      beta = internal::real(c0);
+      typename QRMatrixType::InnerIterator it(m_Q,col);
+      it.valueRef() = 1; //FIXME A row permutation should be performed at this point
+    }
+    else
+    {
+      beta = std::sqrt(internal::abs2(c0) + sqrNorm);
+      if(internal::real(c0) >= RealScalar(0))
+        beta = -beta;
+      typename QRMatrixType::InnerIterator it(m_Q,col);
+      it.valueRef() = 1;
+      for (++it; it; ++it)
+      {
+        it.valueRef() /= (c0 - beta);
+      }
+      tau = internal::conj((beta-c0) / beta);
+        
+    }
+    m_hcoeffs(col) = tau;
+    m_R.insertBackByOuterInnerUnordered(col, col) = beta;
+  }
+  // Finalize the column pointers of the sparse matrices R and Q
+  m_R.finalize(); m_R.makeCompressed();
+  m_Q.finalize(); m_Q.makeCompressed();
+  m_isInitialized = true; 
+  m_factorizationIsok = true;
+  m_info = Success;
+  
+}
+
+namespace internal {
+  
+template<typename _MatrixType, typename OrderingType, typename Rhs>
+struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs>
+  : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs>
+{
+  typedef SparseQR<_MatrixType,OrderingType> Dec;
+  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dec()._solve(rhs(),dst);
+  }
+};
+
+} // end namespace internal
+
+template <typename SparseQRType, typename Derived>
+struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> >
+{
+  typedef typename SparseQRType::QRMatrixType MatrixType;
+  typedef typename SparseQRType::Scalar Scalar;
+  typedef typename SparseQRType::Index Index;
+  // Get the references 
+  SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 
+  m_qr(qr),m_other(other),m_transpose(transpose) {}
+  inline Index rows() const { return m_transpose ? m_qr.rowsQ() : m_qr.cols(); }
+  inline Index cols() const { return m_other.cols(); }
+  
+  // Assign to a vector
+  template<typename DesType>
+  void evalTo(DesType& res) const
+  {
+    Index m = m_qr.rowsQ();
+    Index n = m_qr.cols(); 
+    if (m_transpose)
+    {
+      eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
+      // Compute res = Q' * other :
+      res =  m_other;
+      for (Index k = 0; k < n; k++)
+      {
+        Scalar tau; 
+        // Or alternatively 
+        tau = m_qr.m_Q.col(k).tail(m-k).dot(res.tail(m-k)); 
+        tau = tau * m_qr.m_hcoeffs(k);
+        res -= tau * m_qr.m_Q.col(k);
+      }
+    }
+    else
+    {
+      eigen_assert(m_qr.m_Q.cols() == m_other.rows() && "Non conforming object sizes");
+      // Compute res = Q * other :
+      res = m_other;
+      for (Index k = n-1; k >=0; k--)
+      {
+        Scalar tau;
+        tau = m_qr.m_Q.col(k).tail(m-k).dot(res.tail(m-k));
+        tau = tau * m_qr.m_hcoeffs(k);
+        res -= tau * m_qr.m_Q.col(k);
+      }
+    }
+  }
+  
+  const SparseQRType& m_qr;
+  const Derived& m_other;
+  bool m_transpose;
+};
+
+template<typename SparseQRType>
+struct SparseQRMatrixQReturnType{
+  
+  SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {}
+  template<typename Derived>
+  SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other)
+  {
+    return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false);
+  }
+  SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const
+  {
+    return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
+  }
+  // To use for operations with the transpose of Q
+  SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
+  {
+    return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
+  }
+  const SparseQRType& m_qr;
+};
+
+template<typename SparseQRType>
+struct SparseQRMatrixQTransposeReturnType{
+  SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {}
+  template<typename Derived>
+  SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other)
+  {
+    return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true);
+  }
+  const SparseQRType& m_qr;
+};
+} // end namespace Eigen
+#endif
\ No newline at end of file
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 067f4a83a..f676ee9eb 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -219,6 +219,7 @@ ei_add_test(simplicial_cholesky)
 ei_add_test(conjugate_gradient)
 ei_add_test(bicgstab)
 ei_add_test(sparselu)
+ei_add_test(sparseqr)
 
 # ei_add_test(denseLM)
 
diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp
new file mode 100644
index 000000000..d34f7c48d
--- /dev/null
+++ b/test/sparseqr.cpp
@@ -0,0 +1,62 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+#include "sparse.h"
+#include <Eigen/SparseQR>
+
+
+template<typename MatrixType,typename DenseMat>
+int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300)
+{
+  eigen_assert(maxRows >= maxCols);
+  typedef typename MatrixType::Scalar Scalar;
+  int rows = internal::random<int>(1,maxRows);
+  int cols = internal::random<int>(1,rows);
+  double density = (std::max)(8./(rows*cols), 0.01);
+  
+  A.resize(rows,rows);
+  dA.resize(rows,rows);
+  initSparse<Scalar>(density, dA, A,ForceNonZeroDiag);
+  A.makeCompressed();
+  return rows;
+}
+
+template<typename Scalar> void test_sparseqr_scalar()
+{
+  typedef SparseMatrix<Scalar,ColMajor> MatrixType; 
+  MatrixType A;
+  Matrix<Scalar,Dynamic,Dynamic> dA;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+  DenseVector refX,x,b; 
+  SparseQR<MatrixType, AMDOrdering<int> > solver; 
+  generate_sparse_rectangular_problem(A,dA);
+  
+  int n = A.cols();
+  b = DenseVector::Random(n);
+  solver.compute(A);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse QR factorization failed\n";
+    exit(0);
+    return;
+  }
+  x = solver.solve(b);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse QR factorization failed\n";
+    exit(0);
+    return;
+  }  
+  //Compare with a dense QR solver
+  refX = dA.colPivHouseholderQr().solve(b);
+  VERIFY(x.isApprox(refX,test_precision<Scalar>()));
+}
+void test_sparseqr()
+{
+  CALL_SUBTEST_1(test_sparseqr_scalar<double>());
+  CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >());
+}
\ No newline at end of file