diff options
author | Gael Guennebaud <g.gael@free.fr> | 2012-06-04 13:24:41 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2012-06-04 13:24:41 +0200 |
commit | 945179b26c952f12ca6ea9ad293330cdde7554a0 (patch) | |
tree | 46248287e6deea2c321a62051ee3150bc0d0e069 | |
parent | 5f5a4d4546f785821612a53efafba3552ecb048e (diff) |
CholmodDecomposition now has explicit variants. These variants will allow to provide access to the underlying factors.
-rw-r--r-- | Eigen/CholmodSupport | 10 | ||||
-rw-r--r-- | Eigen/src/CholmodSupport/CholmodSupport.h | 296 | ||||
-rw-r--r-- | test/cholmod_support.cpp | 26 |
3 files changed, 265 insertions, 67 deletions
diff --git a/Eigen/CholmodSupport b/Eigen/CholmodSupport index f26aaf64e..0e1ff48e0 100644 --- a/Eigen/CholmodSupport +++ b/Eigen/CholmodSupport @@ -12,6 +12,16 @@ extern "C" { /** \ingroup Support_modules * \defgroup CholmodSupport_Module CholmodSupport module * + * This module provides an interface to the Cholmod library which is part of the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">suitesparse</a> package. + * It provides the two following main factorization classes: + * - class CholmodSupernodalLLT: a supernodal LLT Cholesky factorization. + * - class CholmodDecomposiiton: a general L(D)LT Cholesky factorization with automatic or explicit runtime selection of the underlying factorization method (supernodal or simplicial). + * + * For the sake of completeness, this module also propose the two following classes: + * - class CholmodSimplicialLLT + * - class CholmodSimplicialLDLT + * Note that these classes does not bring any particular advantage compared to the built-in + * SimplicialLLT and SimplicialLDLT factorization classes. * * \code * #include <Eigen/CholmodSupport> diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h index 5b0c34959..a06c429f0 100644 --- a/Eigen/src/CholmodSupport/CholmodSupport.h +++ b/Eigen/src/CholmodSupport/CholmodSupport.h @@ -160,24 +160,14 @@ enum CholmodMode { CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt }; + /** \ingroup CholmodSupport_Module - * \class CholmodDecomposition - * \brief A Cholesky factorization and solver based on Cholmod - * - * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization - * using the Cholmod library. The sparse matrix A must be selfajoint and positive definite. The vectors or matrices - * X and B can be either dense or sparse. - * - * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> - * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower - * or Upper. Default is Lower. - * - * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; isCompressed() or unisCompressed(). - * - * \sa \ref TutorialSparseDirectSolvers + * \class CholmodBase + * \brief The base class for the direct Cholesky factorization of Cholmod + * \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT */ -template<typename _MatrixType, int _UpLo = Lower> -class CholmodDecomposition +template<typename _MatrixType, int _UpLo, typename Derived> +class CholmodBase : internal::noncopyable { public: typedef _MatrixType MatrixType; @@ -189,21 +179,20 @@ class CholmodDecomposition public: - CholmodDecomposition() + CholmodBase() : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) { cholmod_start(&m_cholmod); - setMode(CholmodLDLt); } - CholmodDecomposition(const MatrixType& matrix) + CholmodBase(const MatrixType& matrix) : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) { cholmod_start(&m_cholmod); compute(matrix); } - ~CholmodDecomposition() + ~CholmodBase() { if(m_cholmodFactor) cholmod_free_factor(&m_cholmodFactor, &m_cholmod); @@ -213,31 +202,8 @@ class CholmodDecomposition inline Index cols() const { return m_cholmodFactor->n; } inline Index rows() const { return m_cholmodFactor->n; } - void setMode(CholmodMode mode) - { - switch(mode) - { - case CholmodAuto: - m_cholmod.final_asis = 1; - m_cholmod.supernodal = CHOLMOD_AUTO; - break; - case CholmodSimplicialLLt: - m_cholmod.final_asis = 0; - m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; - m_cholmod.final_ll = 1; - break; - case CholmodSupernodalLLt: - m_cholmod.final_asis = 1; - m_cholmod.supernodal = CHOLMOD_SUPERNODAL; - break; - case CholmodLDLt: - m_cholmod.final_asis = 1; - m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; - break; - default: - break; - } - } + Derived& derived() { return *static_cast<Derived*>(this); } + const Derived& derived() const { return *static_cast<const Derived*>(this); } /** \brief Reports whether previous computation was successful. * @@ -251,10 +217,11 @@ class CholmodDecomposition } /** Computes the sparse Cholesky decomposition of \a matrix */ - void compute(const MatrixType& matrix) + Derived& compute(const MatrixType& matrix) { analyzePattern(matrix); factorize(matrix); + return derived(); } /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -262,13 +229,13 @@ class CholmodDecomposition * \sa compute() */ template<typename Rhs> - inline const internal::solve_retval<CholmodDecomposition, Rhs> + inline const internal::solve_retval<CholmodBase, Rhs> solve(const MatrixBase<Rhs>& b) const { eigen_assert(m_isInitialized && "LLT is not initialized."); eigen_assert(rows()==b.rows() && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b"); - return internal::solve_retval<CholmodDecomposition, Rhs>(*this, b.derived()); + return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived()); } /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -276,13 +243,13 @@ class CholmodDecomposition * \sa compute() */ template<typename Rhs> - inline const internal::sparse_solve_retval<CholmodDecomposition, Rhs> + inline const internal::sparse_solve_retval<CholmodBase, Rhs> solve(const SparseMatrixBase<Rhs>& b) const { eigen_assert(m_isInitialized && "LLT is not initialized."); eigen_assert(rows()==b.rows() && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b"); - return internal::sparse_solve_retval<CholmodDecomposition, Rhs>(*this, b.derived()); + return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived()); } /** Performs a symbolic decomposition on the sparcity of \a matrix. @@ -372,7 +339,7 @@ class CholmodDecomposition template<typename Stream> void dumpMemory(Stream& s) {} - + protected: mutable cholmod_common m_cholmod; cholmod_factor* m_cholmodFactor; @@ -380,18 +347,225 @@ class CholmodDecomposition bool m_isInitialized; int m_factorizationIsOk; int m_analysisIsOk; +}; + +/** \ingroup CholmodSupport_Module + * \class CholmodSimplicialLLT + * \brief A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod + * + * This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization + * using the Cholmod library. + * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Thefore, it has little practical interest. + * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * X and B can be either dense or sparse. + * + * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> + * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower + * or Upper. Default is Lower. + * + * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. + * + * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLLT + */ +template<typename _MatrixType, int _UpLo = Lower> +class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> > +{ + typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base; + using Base::m_cholmod; - private: - CholmodDecomposition(CholmodDecomposition& ) {} + public: + + typedef _MatrixType MatrixType; + + CholmodSimplicialLLT() : Base() { init(); } + + CholmodSimplicialLLT(const MatrixType& matrix) : Base() + { + init(); + compute(matrix); + } + + ~CholmodSimplicialLLT() {} + protected: + void init() + { + m_cholmod.final_asis = 0; + m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; + m_cholmod.final_ll = 1; + } +}; + + +/** \ingroup CholmodSupport_Module + * \class CholmodSimplicialLDLT + * \brief A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod + * + * This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization + * using the Cholmod library. + * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Thefore, it has little practical interest. + * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * X and B can be either dense or sparse. + * + * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> + * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower + * or Upper. Default is Lower. + * + * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. + * + * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLDLT + */ +template<typename _MatrixType, int _UpLo = Lower> +class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> > +{ + typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base; + using Base::m_cholmod; + + public: + + typedef _MatrixType MatrixType; + + CholmodSimplicialLDLT() : Base() { init(); } + + CholmodSimplicialLDLT(const MatrixType& matrix) : Base() + { + init(); + compute(matrix); + } + + ~CholmodSimplicialLDLT() {} + protected: + void init() + { + m_cholmod.final_asis = 1; + m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; + } +}; + +/** \ingroup CholmodSupport_Module + * \class CholmodSupernodalLLT + * \brief A supernodal Cholesky (LLT) factorization and solver based on Cholmod + * + * This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization + * using the Cholmod library. + * This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM. + * The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * X and B can be either dense or sparse. + * + * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> + * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower + * or Upper. Default is Lower. + * + * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. + * + * \sa \ref TutorialSparseDirectSolvers + */ +template<typename _MatrixType, int _UpLo = Lower> +class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> > +{ + typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base; + using Base::m_cholmod; + + public: + + typedef _MatrixType MatrixType; + + CholmodSupernodalLLT() : Base() { init(); } + + CholmodSupernodalLLT(const MatrixType& matrix) : Base() + { + init(); + compute(matrix); + } + + ~CholmodSupernodalLLT() {} + protected: + void init() + { + m_cholmod.final_asis = 1; + m_cholmod.supernodal = CHOLMOD_SUPERNODAL; + } +}; + +/** \ingroup CholmodSupport_Module + * \class CholmodDecomposition + * \brief A general Cholesky factorization and solver based on Cholmod + * + * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization + * using the Cholmod library. The sparse matrix A must be selfajoint and positive definite. The vectors or matrices + * X and B can be either dense or sparse. + * + * This variant permits to change the underlying Cholesky method at runtime. + * On the other hand, it does not provide access to the result of the factorization. + * The default is to let Cholmod automatically choose between a simplicial and supernodal factorization. + * + * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> + * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower + * or Upper. Default is Lower. + * + * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed. + * + * \sa \ref TutorialSparseDirectSolvers + */ +template<typename _MatrixType, int _UpLo = Lower> +class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> > +{ + typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base; + using Base::m_cholmod; + + public: + + typedef _MatrixType MatrixType; + + CholmodDecomposition() : Base() { init(); } + + CholmodDecomposition(const MatrixType& matrix) : Base() + { + init(); + compute(matrix); + } + + ~CholmodDecomposition() {} + + void setMode(CholmodMode mode) + { + switch(mode) + { + case CholmodAuto: + m_cholmod.final_asis = 1; + m_cholmod.supernodal = CHOLMOD_AUTO; + break; + case CholmodSimplicialLLt: + m_cholmod.final_asis = 0; + m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; + m_cholmod.final_ll = 1; + break; + case CholmodSupernodalLLt: + m_cholmod.final_asis = 1; + m_cholmod.supernodal = CHOLMOD_SUPERNODAL; + break; + case CholmodLDLt: + m_cholmod.final_asis = 1; + m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; + break; + default: + break; + } + } + protected: + void init() + { + m_cholmod.final_asis = 1; + m_cholmod.supernodal = CHOLMOD_AUTO; + } }; namespace internal { -template<typename _MatrixType, int _UpLo, typename Rhs> -struct solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs> - : solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs> +template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs> +struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> + : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> { - typedef CholmodDecomposition<_MatrixType,_UpLo> Dec; + typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec; EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) template<typename Dest> void evalTo(Dest& dst) const @@ -400,11 +574,11 @@ struct solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs> } }; -template<typename _MatrixType, int _UpLo, typename Rhs> -struct sparse_solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs> - : sparse_solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs> +template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs> +struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> + : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> { - typedef CholmodDecomposition<_MatrixType,_UpLo> Dec; + typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec; EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) template<typename Dest> void evalTo(Dest& dst) const diff --git a/test/cholmod_support.cpp b/test/cholmod_support.cpp index 2dd0acc32..1ebba2438 100644 --- a/test/cholmod_support.cpp +++ b/test/cholmod_support.cpp @@ -28,13 +28,27 @@ template<typename T> void test_cholmod_T() { - CholmodDecomposition<SparseMatrix<T>, Lower> chol_colmajor_lower; chol_colmajor_lower.setMode(CholmodSupernodalLLt); - CholmodDecomposition<SparseMatrix<T>, Upper> chol_colmajor_upper; chol_colmajor_upper.setMode(CholmodSupernodalLLt); - CholmodDecomposition<SparseMatrix<T>, Lower> llt_colmajor_lower; llt_colmajor_lower.setMode(CholmodSimplicialLLt); - CholmodDecomposition<SparseMatrix<T>, Upper> llt_colmajor_upper; llt_colmajor_upper.setMode(CholmodSimplicialLLt); - CholmodDecomposition<SparseMatrix<T>, Lower> ldlt_colmajor_lower; ldlt_colmajor_lower.setMode(CholmodLDLt); - CholmodDecomposition<SparseMatrix<T>, Upper> ldlt_colmajor_upper; ldlt_colmajor_upper.setMode(CholmodLDLt); + CholmodDecomposition<SparseMatrix<T>, Lower> g_chol_colmajor_lower; g_chol_colmajor_lower.setMode(CholmodSupernodalLLt); + CholmodDecomposition<SparseMatrix<T>, Upper> g_chol_colmajor_upper; g_chol_colmajor_upper.setMode(CholmodSupernodalLLt); + CholmodDecomposition<SparseMatrix<T>, Lower> g_llt_colmajor_lower; g_llt_colmajor_lower.setMode(CholmodSimplicialLLt); + CholmodDecomposition<SparseMatrix<T>, Upper> g_llt_colmajor_upper; g_llt_colmajor_upper.setMode(CholmodSimplicialLLt); + CholmodDecomposition<SparseMatrix<T>, Lower> g_ldlt_colmajor_lower; g_ldlt_colmajor_lower.setMode(CholmodLDLt); + CholmodDecomposition<SparseMatrix<T>, Upper> g_ldlt_colmajor_upper; g_ldlt_colmajor_upper.setMode(CholmodLDLt); + + CholmodSupernodalLLT<SparseMatrix<T>, Lower> chol_colmajor_lower; + CholmodSupernodalLLT<SparseMatrix<T>, Upper> chol_colmajor_upper; + CholmodSimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower; + CholmodSimplicialLLT<SparseMatrix<T>, Upper> llt_colmajor_upper; + CholmodSimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower; + CholmodSimplicialLDLT<SparseMatrix<T>, Upper> ldlt_colmajor_upper; + check_sparse_spd_solving(g_chol_colmajor_lower); + check_sparse_spd_solving(g_chol_colmajor_upper); + check_sparse_spd_solving(g_llt_colmajor_lower); + check_sparse_spd_solving(g_llt_colmajor_upper); + check_sparse_spd_solving(g_ldlt_colmajor_lower); + check_sparse_spd_solving(g_ldlt_colmajor_upper); + check_sparse_spd_solving(chol_colmajor_lower); check_sparse_spd_solving(chol_colmajor_upper); check_sparse_spd_solving(llt_colmajor_lower); |