Sparse module: refactoring of the cholesky factorization,

now the backends are well separated from the default impl, etc.

4 files changed, 241 insertions, 60 deletions

diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index 0a57badc7..ed031b264 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/Eigen/src/Sparse/CholmodSupport.h
@@ -99,25 +99,109 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
 }
 
 template<typename MatrixType>
-void SparseCholesky<MatrixType>::computeUsingCholmod(const MatrixType& a)
+class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
 {
-  cholmod_common c;
-  cholmod_start(&c);
+  protected:
+    typedef SparseCholesky<MatrixType> Base;
+    using Base::Scalar;
+    using Base::RealScalar;
+    using Base::MatrixLIsDirty;
+    using Base::SupernodalFactorIsDirty;
+    using Base::m_flags;
+    using Base::m_matrix;
+    using Base::m_status;
+
+  public:
+
+    SparseCholesky(const MatrixType& matrix, int flags = 0)
+      : Base(matrix, flags), m_cholmodFactor(0)
+    {
+      cholmod_start(&m_cholmod);
+      compute(matrix);
+    }
+
+    ~SparseCholesky()
+    {
+      if (m_cholmodFactor)
+        cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
+      cholmod_finish(&m_cholmod);
+    }
+
+    inline const typename Base::CholMatrixType& matrixL(void) const;
+
+    template<typename Derived>
+    void solveInPlace(MatrixBase<Derived> &b) const;
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    mutable cholmod_common m_cholmod;
+    cholmod_factor* m_cholmodFactor;
+};
+
+template<typename MatrixType>
+void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
+{
+  if (m_cholmodFactor)
+  {
+    cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
+    m_cholmodFactor = 0;
+  }
+
   cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
-  if (!(m_flags&CholPartial))
+  if (m_flags&IncompleteFactorization)
+  {
+    m_cholmod.nmethods = 1;
+    m_cholmod.method [0].ordering = CHOLMOD_NATURAL;
+    m_cholmod.postorder = 0;
+  }
+  else
+  {
+    m_cholmod.nmethods = 1;
+    m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
+    m_cholmod.postorder = 0;
+  }
+  m_cholmod.final_ll = 1;
+  m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
+  cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
+
+  m_status = (m_status & ~SupernodalFactorIsDirty) | MatrixLIsDirty;
+}
+
+template<typename MatrixType>
+inline const typename SparseCholesky<MatrixType>::CholMatrixType&
+SparseCholesky<MatrixType,Cholmod>::matrixL() const
+{
+  if (m_status & MatrixLIsDirty)
+  {
+    ei_assert(!(m_status & SupernodalFactorIsDirty));
+
+    cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
+    const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*cmRes);
+    free(cmRes);
+
+    m_status = (m_status & ~MatrixLIsDirty);
+  }
+  return m_matrix;
+}
+
+template<typename MatrixType>
+template<typename Derived>
+void SparseCholesky<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows());
+
+  if (m_status & MatrixLIsDirty)
+  {
+//     ei_assert(!(m_status & SupernodalFactorIsDirty));
+//     taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
+    matrixL();
+  }
+//   else
   {
-    c.nmethods = 1;
-    c.method [0].ordering = CHOLMOD_NATURAL;
-    c.postorder = 0;
+    Base::solveInPlace(b);
   }
-  c.final_ll = 1;
-  cholmod_factor *L = cholmod_analyze(&A, &c);
-  cholmod_factorize(&A, L, &c);
-  cholmod_sparse* cmRes = cholmod_factor_to_sparse(L, &c);
-  m_matrix = CholMatrixType::Map(*cmRes);
-  free(cmRes);
-  cholmod_free_factor(&L, &c);
-  cholmod_finish(&c);
 }
 
 #endif // EIGEN_CHOLMODSUPPORT_H
diff --git a/Eigen/src/Sparse/SparseCholesky.h b/Eigen/src/Sparse/SparseCholesky.h
index 5af84ebed..d8d9f6ab9 100644
--- a/Eigen/src/Sparse/SparseCholesky.h
+++ b/Eigen/src/Sparse/SparseCholesky.h
@@ -25,12 +25,25 @@
 #ifndef EIGEN_SPARSECHOLESKY_H
 #define EIGEN_SPARSECHOLESKY_H
 
+enum SparseBackend {
+  DefaultBackend,
+  Taucs,
+  Cholmod,
+  SuperLU
+};
+
 enum {
-  CholFull        = 0x0,  // full is the default
-  CholPartial     = 0x1,
+  CompleteFactorization        = 0x0,  // full is the default
+  IncompleteFactorization      = 0x1,
+  MemoryEfficient              = 0x2,
+  SupernodalMultifrontal       = 0x4,
+  SupernodalLeftLooking        = 0x8,
+
+
+/*
   CholUseEigen    = 0x0,  // Eigen's impl is the default
   CholUseTaucs    = 0x2,
-  CholUseCholmod  = 0x4,
+  CholUseCholmod  = 0x4*/
 };
 
 /** \ingroup Sparse_Module
@@ -43,23 +56,22 @@ enum {
   *
   * \sa class Cholesky, class CholeskyWithoutSquareRoot
   */
-template<typename MatrixType> class SparseCholesky
+template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
 {
-  private:
+  protected:
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
-    typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
     typedef SparseMatrix<Scalar,Lower> CholMatrixType;
 
     enum {
-      PacketSize = ei_packet_traits<Scalar>::size,
-      AlignmentMask = int(PacketSize)-1
+      SupernodalFactorIsDirty      = 0x10000,
+      MatrixLIsDirty               = 0x20000,
     };
 
   public:
 
     SparseCholesky(const MatrixType& matrix, int flags = 0)
-      : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags)
+      : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
     {
       compute(matrix);
     }
@@ -69,42 +81,26 @@ template<typename MatrixType> class SparseCholesky
     /** \returns true if the matrix is positive definite */
     inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
 
-    // TODO impl the solver
-//     template<typename Derived>
-//     typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
+    template<typename Derived>
+    void solveInPlace(MatrixBase<Derived> &b) const;
 
     void compute(const MatrixType& matrix);
 
   protected:
-    void computeUsingEigen(const MatrixType& matrix);
-    void computeUsingTaucs(const MatrixType& matrix);
-    void computeUsingCholmod(const MatrixType& matrix);
-
-  protected:
     /** \internal
       * Used to compute and store L
       * The strict upper part is not used and even not initialized.
       */
     CholMatrixType m_matrix;
     int m_flags;
+    mutable int m_status;
     bool m_isPositiveDefinite;
 };
 
 /** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
   */
-template<typename MatrixType>
-void SparseCholesky<MatrixType>::compute(const MatrixType& a)
-{
-  if (m_flags&CholUseTaucs)
-    computeUsingTaucs(a);
-  else if (m_flags&CholUseCholmod)
-    computeUsingCholmod(a);
-  else
-    computeUsingEigen(a);
-}
-
-template<typename MatrixType>
-void SparseCholesky<MatrixType>::computeUsingEigen(const MatrixType& a)
+template<typename MatrixType, int Backend>
+void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
 {
   assert(a.rows()==a.cols());
   const int size = a.rows();
@@ -173,4 +169,15 @@ void SparseCholesky<MatrixType>::computeUsingEigen(const MatrixType& a)
   m_matrix.endFill();
 }
 
+template<typename MatrixType, int Backend>
+template<typename Derived>
+void SparseCholesky<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows());
+
+  m_matrix.solveTriangularInPlace(b);
+  m_matrix.adjoint().solveTriangularInPlace(b);
+}
+
 #endif // EIGEN_BASICSPARSECHOLESKY_H
diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h
index 55d2892ec..a0a5b4b71 100644
--- a/Eigen/src/Sparse/TaucsSupport.h
+++ b/Eigen/src/Sparse/TaucsSupport.h
@@ -79,12 +79,112 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& tau
 }
 
 template<typename MatrixType>
-void SparseCholesky<MatrixType>::computeUsingTaucs(const MatrixType& a)
+class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
 {
-  taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
-  taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
-  m_matrix = CholMatrixType::Map(*taucsRes);
-  free(taucsRes);
+  protected:
+    typedef SparseCholesky<MatrixType> Base;
+    using Base::Scalar;
+    using Base::RealScalar;
+    using Base::MatrixLIsDirty;
+    using Base::SupernodalFactorIsDirty;
+    using Base::m_flags;
+    using Base::m_matrix;
+    using Base::m_status;
+
+  public:
+
+    SparseCholesky(const MatrixType& matrix, int flags = 0)
+      : Base(matrix, flags), m_taucsSupernodalFactor(0)
+    {
+      compute(matrix);
+    }
+
+    ~SparseCholesky()
+    {
+      if (m_taucsSupernodalFactor)
+        taucs_supernodal_factor_free(m_taucsSupernodalFactor);
+    }
+
+    inline const typename Base::CholMatrixType& matrixL(void) const;
+
+    template<typename Derived>
+    void solveInPlace(MatrixBase<Derived> &b) const;
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+    void* m_taucsSupernodalFactor;
+};
+
+template<typename MatrixType>
+void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
+{
+  if (m_taucsSupernodalFactor)
+  {
+    taucs_supernodal_factor_free(m_taucsSupernodalFactor);
+    m_taucsSupernodalFactor = 0;
+  }
+
+  if (m_flags & IncompleteFactorization)
+  {
+    taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
+    taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
+    m_matrix = Base::CholMatrixType::Map(*taucsRes);
+    free(taucsRes);
+    m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
+             | IncompleteFactorization
+             | SupernodalFactorIsDirty;
+  }
+  else
+  {
+    taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
+    if ( (m_flags & SupernodalLeftLooking)
+      || ((!(m_flags & SupernodalMultifrontal)) && (m_flags & MemoryEfficient)) )
+    {
+      m_taucsSupernodalFactor = taucs_ccs_factor_llt_ll(&taucsMatA);
+    }
+    else
+    {
+      // use the faster Multifrontal routine
+      m_taucsSupernodalFactor = taucs_ccs_factor_llt_ll(&taucsMatA);
+    }
+    m_status = (m_status & ~IncompleteFactorization) | CompleteFactorization | MatrixLIsDirty;
+  }
+}
+
+template<typename MatrixType>
+inline const typename SparseCholesky<MatrixType>::CholMatrixType&
+SparseCholesky<MatrixType,Taucs>::matrixL() const
+{
+  if (m_status & MatrixLIsDirty)
+  {
+    ei_assert(!(m_status & SupernodalFactorIsDirty));
+
+    taucs_ccs_matrix* taucsL = taucs_supernodal_factor_to_ccs(m_taucsSupernodalFactor);
+    const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*taucsL);
+    free(taucsL);
+    m_status = (m_status & ~MatrixLIsDirty);
+  }
+  return m_matrix;
+}
+
+template<typename MatrixType>
+template<typename Derived>
+void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
+{
+  const int size = m_matrix.rows();
+  ei_assert(size==b.rows());
+
+  if (m_status & MatrixLIsDirty)
+  {
+//     ei_assert(!(m_status & SupernodalFactorIsDirty));
+//     taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
+    matrixL();
+  }
+//   else
+  {
+    Base::solveInPlace(b);
+  }
 }
 
 #endif // EIGEN_TAUCSSUPPORT_H
diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h
index 5fb8396f2..3a1a33b3c 100644
--- a/Eigen/src/Sparse/TriangularSolver.h
+++ b/Eigen/src/Sparse/TriangularSolver.h
@@ -25,16 +25,6 @@
 #ifndef EIGEN_SPARSETRIANGULARSOLVER_H
 #define EIGEN_SPARSETRIANGULARSOLVER_H
 
-// template<typename Lhs, typename Rhs,
-//   int TriangularPart = (int(Lhs::Flags) & LowerTriangularBit)
-//                      ? Lower
-//                      : (int(Lhs::Flags) & UpperTriangularBit)
-//                      ? Upper
-//                      : -1,
-//   int StorageOrder = int(Lhs::Flags) & RowMajorBit ? RowMajor : ColMajor
-//   >
-// struct ei_sparse_trisolve_selector;
-
 // forward substitution, row-major
 template<typename Lhs, typename Rhs>
 struct ei_solve_triangular_selector<Lhs,Rhs,Lower,RowMajor|IsSparse>