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

2 files changed, 487 insertions, 0 deletions

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