From 22507fa6452a1de4d3af9b4fe0300e21599838e0 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Sun, 5 Oct 2008 20:19:47 +0000 Subject: Sparse module: refactoring of the cholesky factorization, now the backends are well separated from the default impl, etc. --- Eigen/src/Sparse/CholmodSupport.h | 114 +++++++++++++++++++++++++++++++----- Eigen/src/Sparse/SparseCholesky.h | 67 +++++++++++---------- Eigen/src/Sparse/TaucsSupport.h | 110 ++++++++++++++++++++++++++++++++-- Eigen/src/Sparse/TriangularSolver.h | 10 ---- 4 files changed, 241 insertions(+), 60 deletions(-) (limited to 'Eigen') 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 SparseMatrix::Map(cholmod_sparse& cm) } template -void SparseCholesky::computeUsingCholmod(const MatrixType& a) +class SparseCholesky : public SparseCholesky { - cholmod_common c; - cholmod_start(&c); + protected: + typedef SparseCholesky 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 + void solveInPlace(MatrixBase &b) const; + + void compute(const MatrixType& matrix); + + protected: + mutable cholmod_common m_cholmod; + cholmod_factor* m_cholmodFactor; +}; + +template +void SparseCholesky::compute(const MatrixType& a) +{ + if (m_cholmodFactor) + { + cholmod_free_factor(&m_cholmodFactor, &m_cholmod); + m_cholmodFactor = 0; + } + cholmod_sparse A = const_cast(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 +inline const typename SparseCholesky::CholMatrixType& +SparseCholesky::matrixL() const +{ + if (m_status & MatrixLIsDirty) + { + ei_assert(!(m_status & SupernodalFactorIsDirty)); + + cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod); + const_cast(m_matrix) = Base::CholMatrixType::Map(*cmRes); + free(cmRes); + + m_status = (m_status & ~MatrixLIsDirty); + } + return m_matrix; +} + +template +template +void SparseCholesky::solveInPlace(MatrixBase &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 class SparseCholesky +template class SparseCholesky { - private: + protected: typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; - typedef Matrix VectorType; typedef SparseMatrix CholMatrixType; enum { - PacketSize = ei_packet_traits::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,17 +81,11 @@ template 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::Eval solve(const MatrixBase &b) const; + template + void solveInPlace(MatrixBase &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 @@ -87,24 +93,14 @@ template class SparseCholesky */ 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 -void SparseCholesky::compute(const MatrixType& a) -{ - if (m_flags&CholUseTaucs) - computeUsingTaucs(a); - else if (m_flags&CholUseCholmod) - computeUsingCholmod(a); - else - computeUsingEigen(a); -} - -template -void SparseCholesky::computeUsingEigen(const MatrixType& a) +template +void SparseCholesky::compute(const MatrixType& a) { assert(a.rows()==a.cols()); const int size = a.rows(); @@ -173,4 +169,15 @@ void SparseCholesky::computeUsingEigen(const MatrixType& a) m_matrix.endFill(); } +template +template +void SparseCholesky::solveInPlace(MatrixBase &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 SparseMatrix::Map(taucs_ccs_matrix& tau } template -void SparseCholesky::computeUsingTaucs(const MatrixType& a) +class SparseCholesky : public SparseCholesky { - taucs_ccs_matrix taucsMatA = const_cast(a).asTaucsMatrix(); - taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0); - m_matrix = CholMatrixType::Map(*taucsRes); - free(taucsRes); + protected: + typedef SparseCholesky 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 + void solveInPlace(MatrixBase &b) const; + + void compute(const MatrixType& matrix); + + protected: + void* m_taucsSupernodalFactor; +}; + +template +void SparseCholesky::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(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(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 +inline const typename SparseCholesky::CholMatrixType& +SparseCholesky::matrixL() const +{ + if (m_status & MatrixLIsDirty) + { + ei_assert(!(m_status & SupernodalFactorIsDirty)); + + taucs_ccs_matrix* taucsL = taucs_supernodal_factor_to_ccs(m_taucsSupernodalFactor); + const_cast(m_matrix) = Base::CholMatrixType::Map(*taucsL); + free(taucsL); + m_status = (m_status & ~MatrixLIsDirty); + } + return m_matrix; +} + +template +template +void SparseCholesky::solveInPlace(MatrixBase &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 -// struct ei_sparse_trisolve_selector; - // forward substitution, row-major template struct ei_solve_triangular_selector -- cgit v1.2.3