GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-13 20:26:03 +08:00
Code Issues Packages Projects Releases Wiki Activity

Make the ordering method of SimplicialL[D]LT user configurable.

Browse Source
This commit is contained in:
Gael Guennebaud 2014-07-20 14:22:58 +02:00
parent 8e19027130
commit d4cc1bdc7f
2 changed files with 46 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
Eigen/src/SparseCholesky/SimplicialCholesky.h
View File

@ -37,6 +37,7 @@ class SimplicialCholeskyBase : internal::noncopyable
{ {
public: public:
typedef typename internal::traits<Derived>::MatrixType MatrixType; typedef typename internal::traits<Derived>::MatrixType MatrixType;
typedef typename internal::traits<Derived>::OrderingType OrderingType;
enum { UpLo = internal::traits<Derived>::UpLo }; enum { UpLo = internal::traits<Derived>::UpLo };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
@ -240,15 +241,16 @@ class SimplicialCholeskyBase : internal::noncopyable
RealScalar m_shiftScale; RealScalar m_shiftScale;
}; };
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT; template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT; template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky; template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky;
namespace internal { namespace internal {
template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> > template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef _Ordering OrderingType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -259,9 +261,10 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixTyp
static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
}; };
template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> > template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef _Ordering OrderingType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -272,9 +275,10 @@ template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixTyp
static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
}; };
template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> > template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef _Ordering OrderingType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
}; };
@ -294,11 +298,12 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_Matr
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \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 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower. * or Upper. Default is Lower.
* \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
* *
* \sa class SimplicialLDLT * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo, typename _Ordering>
class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> > class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
@ -382,11 +387,12 @@ public:
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \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 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower. * or Upper. Default is Lower.
* \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
* *
* \sa class SimplicialLLT * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo, typename _Ordering>
class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> > class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
@ -467,8 +473,8 @@ public:
* *
* \sa class SimplicialLDLT, class SimplicialLLT * \sa class SimplicialLDLT, class SimplicialLLT
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo, typename _Ordering>
class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> > class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
@ -612,15 +618,13 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixTy
{ {
eigen_assert(a.rows()==a.cols()); eigen_assert(a.rows()==a.cols());
const Index size = a.rows(); const Index size = a.rows();
// TODO allows to configure the permutation
// Note that amd compute the inverse permutation // Note that amd compute the inverse permutation
{ {
CholMatrixType C; CholMatrixType C;
C = a.template selfadjointView<UpLo>(); C = a.template selfadjointView<UpLo>();
// remove diagonal entries:
// seems not to be needed OrderingType ordering;
// C.prune(keep_diag()); ordering(C,m_Pinv);
internal::minimum_degree_ordering(C, m_Pinv);
} }
if(m_Pinv.size()>0) if(m_Pinv.size()>0)

41
test/simplicial_cholesky.cpp
View File

@ -11,26 +11,31 @@
template<typename T> void test_simplicial_cholesky_T() template<typename T> void test_simplicial_cholesky_T()
{ {
SimplicialCholesky<SparseMatrix<T>, Lower> chol_colmajor_lower; SimplicialCholesky<SparseMatrix<T>, Lower> chol_colmajor_lower_amd;
SimplicialCholesky<SparseMatrix<T>, Upper> chol_colmajor_upper; SimplicialCholesky<SparseMatrix<T>, Upper> chol_colmajor_upper_amd;
SimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower; SimplicialLLT<SparseMatrix<T>, Lower> llt_colmajor_lower_amd;
SimplicialLDLT<SparseMatrix<T>, Upper> llt_colmajor_upper; SimplicialLLT<SparseMatrix<T>, Upper> llt_colmajor_upper_amd;
SimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower; SimplicialLDLT<SparseMatrix<T>, Lower> ldlt_colmajor_lower_amd;
SimplicialLDLT<SparseMatrix<T>, Upper> ldlt_colmajor_upper; SimplicialLDLT<SparseMatrix<T>, Upper> ldlt_colmajor_upper_amd;
SimplicialLDLT<SparseMatrix<T>, Lower, NaturalOrdering<int> > ldlt_colmajor_lower_nat;
SimplicialLDLT<SparseMatrix<T>, Upper, NaturalOrdering<int> > ldlt_colmajor_upper_nat;
check_sparse_spd_solving(chol_colmajor_lower); check_sparse_spd_solving(chol_colmajor_lower_amd);
check_sparse_spd_solving(chol_colmajor_upper); check_sparse_spd_solving(chol_colmajor_upper_amd);
check_sparse_spd_solving(llt_colmajor_lower); check_sparse_spd_solving(llt_colmajor_lower_amd);
check_sparse_spd_solving(llt_colmajor_upper); check_sparse_spd_solving(llt_colmajor_upper_amd);
check_sparse_spd_solving(ldlt_colmajor_lower); check_sparse_spd_solving(ldlt_colmajor_lower_amd);
check_sparse_spd_solving(ldlt_colmajor_upper); check_sparse_spd_solving(ldlt_colmajor_upper_amd);
check_sparse_spd_determinant(chol_colmajor_lower); check_sparse_spd_determinant(chol_colmajor_lower_amd);
check_sparse_spd_determinant(chol_colmajor_upper); check_sparse_spd_determinant(chol_colmajor_upper_amd);
check_sparse_spd_determinant(llt_colmajor_lower); check_sparse_spd_determinant(llt_colmajor_lower_amd);
check_sparse_spd_determinant(llt_colmajor_upper); check_sparse_spd_determinant(llt_colmajor_upper_amd);
check_sparse_spd_determinant(ldlt_colmajor_lower); check_sparse_spd_determinant(ldlt_colmajor_lower_amd);
check_sparse_spd_determinant(ldlt_colmajor_upper); check_sparse_spd_determinant(ldlt_colmajor_upper_amd);
check_sparse_spd_solving(ldlt_colmajor_lower_nat);
check_sparse_spd_solving(ldlt_colmajor_upper_nat);
} }
void test_simplicial_cholesky() void test_simplicial_cholesky()