CholmodDecomposition now has explicit variants. These variants will allow to provide access to the underlying factors.

This commit is contained in:
Gael Guennebaud 2012-06-04 13:24:41 +02:00
parent 5f5a4d4546
commit 945179b26c
3 changed files with 265 additions and 67 deletions

View File

@ -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>

View File

@ -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

View File

@ -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);