GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-01 16:24:28 +08:00
Code Issues Packages Projects Releases Wiki Activity

For consistency, Simplicial* now factorizes P A P^-1 (instead of P^-1 A P).

Browse Source 
Document how is applied the permutation in Simplicial* .
This commit is contained in:
Gael Guennebaud 2012-06-07 16:24:46 +02:00
parent 63c6ab3e42
commit 1e5e66b642
1 changed files with 20 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -77,6 +77,9 @@ enum SimplicialCholeskyMode {
* selfadjoint and positive definite. The factorization allows for solving A.X = B where * selfadjoint and positive definite. The factorization allows for solving A.X = B where
* X and B can be either dense or sparse. * X and B can be either dense or sparse.
* *
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
* such that the factorized matrix is P A P^-1.
*
* \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.
@ -208,7 +211,7 @@ class SimplicialCholeskyBase : internal::noncopyable
return; return;
if(m_P.size()>0) if(m_P.size()>0)
dest = m_Pinv * b; dest = m_P * b;
else else
dest = b; dest = b;
@ -222,7 +225,7 @@ class SimplicialCholeskyBase : internal::noncopyable
derived().matrixU().solveInPlace(dest); derived().matrixU().solveInPlace(dest);
if(m_P.size()>0) if(m_P.size()>0)
dest = m_P * dest; dest = m_Pinv * dest;
} }
/** \internal */ /** \internal */
@ -268,7 +271,7 @@ class SimplicialCholeskyBase : internal::noncopyable
eigen_assert(a.rows()==a.cols()); eigen_assert(a.rows()==a.cols());
int size = a.cols(); int size = a.cols();
CholMatrixType ap(size,size); CholMatrixType ap(size,size);
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv); ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
factorize_preordered<DoLDLT>(ap); factorize_preordered<DoLDLT>(ap);
} }
@ -359,6 +362,9 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_Matr
* selfadjoint and positive definite. The factorization allows for solving A.X = B where * selfadjoint and positive definite. The factorization allows for solving A.X = B where
* X and B can be either dense or sparse. * X and B can be either dense or sparse.
* *
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
* such that the factorized matrix is P A P^-1.
*
* \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.
@ -444,6 +450,9 @@ public:
* selfadjoint and positive definite. The factorization allows for solving A.X = B where * selfadjoint and positive definite. The factorization allows for solving A.X = B where
* X and B can be either dense or sparse. * X and B can be either dense or sparse.
* *
* In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
* such that the factorized matrix is P A P^-1.
*
* \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.
@ -628,7 +637,7 @@ public:
return; return;
if(Base::m_P.size()>0) if(Base::m_P.size()>0)
dest = Base::m_Pinv * b; dest = Base::m_P * b;
else else
dest = b; dest = b;
@ -652,7 +661,7 @@ public:
} }
if(Base::m_P.size()>0) if(Base::m_P.size()>0)
dest = Base::m_P * dest; dest = Base::m_Pinv * dest;
} }
Scalar determinant() const Scalar determinant() const
@ -678,22 +687,23 @@ 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 // TODO allows to configure the 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: // remove diagonal entries:
// seems not to be needed // seems not to be needed
// C.prune(keep_diag()); // C.prune(keep_diag());
internal::minimum_degree_ordering(C, m_P); internal::minimum_degree_ordering(C, m_Pinv);
} }
if(m_P.size()>0) if(m_Pinv.size()>0)
m_Pinv = m_P.inverse(); m_P = m_Pinv.inverse();
else else
m_Pinv.resize(0); m_P.resize(0);
ap.resize(size,size); ap.resize(size,size);
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv); ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
} }
template<typename Derived> template<typename Derived>