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

In simplicial cholesky: avoid deep copy of the input matrix is this later can be used readily

Browse Source
This commit is contained in:
Gael Guennebaud 2014-12-08 17:56:33 +01:00
parent a910a7466e
commit 41a20994cc
1 changed files with 53 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -17,6 +17,27 @@ enum SimplicialCholeskyMode {
SimplicialCholeskyLDLT SimplicialCholeskyLDLT
}; };
namespace internal {
template<typename CholMatrixType, typename InputMatrixType>
struct simplicial_cholesky_grab_input {
typedef CholMatrixType const * ConstCholMatrixPtr;
static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
{
tmp = input;
pmat = &tmp;
}
};
template<typename MatrixType>
struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
typedef MatrixType const * ConstMatrixPtr;
static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/)
{
pmat = &input;
}
};
} // end namespace internal
/** \ingroup SparseCholesky_Module /** \ingroup SparseCholesky_Module
* \brief A direct sparse Cholesky factorizations * \brief A direct sparse Cholesky factorizations
* *
@ -46,6 +67,7 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived>
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
typedef CholMatrixType const * ConstCholMatrixPtr;
typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Scalar,Dynamic,1> VectorType;
public: public:
@ -169,10 +191,11 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived>
{ {
eigen_assert(matrix.rows()==matrix.cols()); eigen_assert(matrix.rows()==matrix.cols());
Index size = matrix.cols(); Index size = matrix.cols();
CholMatrixType ap(size,size); CholMatrixType tmp(size,size);
ordering(matrix, ap); ConstCholMatrixPtr pmat;
analyzePattern_preordered(ap, DoLDLT); ordering(matrix, pmat, tmp);
factorize_preordered<DoLDLT>(ap); analyzePattern_preordered(*pmat, DoLDLT);
factorize_preordered<DoLDLT>(*pmat);
} }
template<bool DoLDLT> template<bool DoLDLT>
@ -180,9 +203,21 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived>
{ {
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 tmp(size,size);
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); ConstCholMatrixPtr pmat;
factorize_preordered<DoLDLT>(ap);
if(m_P.size()==0 && (UpLo&Upper)==Upper)
{
// If there is no ordering, try to directly use the input matrix without any copy
internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
}
else
{
tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
pmat = &tmp;
}
factorize_preordered<DoLDLT>(*pmat);
} }
template<bool DoLDLT> template<bool DoLDLT>
@ -192,13 +227,14 @@ class SimplicialCholeskyBase : public SparseSolverBase<Derived>
{ {
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 tmp(size,size);
ordering(a, ap); ConstCholMatrixPtr pmat;
analyzePattern_preordered(ap,doLDLT); ordering(a, pmat, tmp);
analyzePattern_preordered(*pmat,doLDLT);
} }
void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
void ordering(const MatrixType& a, CholMatrixType& ap); void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
/** keeps off-diagonal entries; drops diagonal entries */ /** keeps off-diagonal entries; drops diagonal entries */
struct keep_diag { struct keep_diag {
@ -603,10 +639,11 @@ public:
}; };
template<typename Derived> template<typename Derived>
void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap) void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
{ {
eigen_assert(a.rows()==a.cols()); eigen_assert(a.rows()==a.cols());
const Index size = a.rows(); const Index size = a.rows();
pmat = &ap;
// Note that ordering methods compute the inverse permutation // Note that ordering methods compute the inverse permutation
if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value) if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
{ {
@ -628,13 +665,14 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixTy
{ {
m_Pinv.resize(0); m_Pinv.resize(0);
m_P.resize(0); m_P.resize(0);
if(UpLo==Lower) if(UpLo==Lower || MatrixType::IsRowMajor)
{ {
// we have to transpose the lower part to to the upper one
ap.resize(size,size); ap.resize(size,size);
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
} }
else else
ap = a.template triangularView<Upper>(); internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
} }
} }