mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Sparse module: refactoring of the cholesky factorization,

Browse Source 
now the backends are well separated from the default impl, etc.
This commit is contained in:
Gael Guennebaud 2008-10-05 20:19:47 +00:00
parent b8fc1edb2c
commit 22507fa645
4 changed files with 242 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

116
Eigen/src/Sparse/CholmodSupport.h
View File

@ -99,25 +99,109 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
}
template<typename MatrixType>
void SparseCholesky<MatrixType>::computeUsingCholmod(const MatrixType& a)
class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
{
cholmod_common c;
cholmod_start(&c);
cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
if (!(m_flags&CholPartial))
protected:
typedef SparseCholesky<MatrixType> Base;
using Base::Scalar;
using Base::RealScalar;
using Base::MatrixLIsDirty;
using Base::SupernodalFactorIsDirty;
using Base::m_flags;
using Base::m_matrix;
using Base::m_status;
public:
SparseCholesky(const MatrixType& matrix, int flags = 0)
: Base(matrix, flags), m_cholmodFactor(0)
{
cholmod_start(&m_cholmod);
compute(matrix);
}
~SparseCholesky()
{
if (m_cholmodFactor)
cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
cholmod_finish(&m_cholmod);
}
inline const typename Base::CholMatrixType& matrixL(void) const;
template<typename Derived>
void solveInPlace(MatrixBase<Derived> &b) const;
void compute(const MatrixType& matrix);
protected:
mutable cholmod_common m_cholmod;
cholmod_factor* m_cholmodFactor;
};
template<typename MatrixType>
void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
{
if (m_cholmodFactor)
{
c.nmethods = 1;
c.method [0].ordering = CHOLMOD_NATURAL;
c.postorder = 0;
cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
m_cholmodFactor = 0;
}
cholmod_sparse A = const_cast<MatrixType&>(a).asCholmodMatrix();
if (m_flags&IncompleteFactorization)
{
m_cholmod.nmethods = 1;
m_cholmod.method [0].ordering = CHOLMOD_NATURAL;
m_cholmod.postorder = 0;
}
else
{
m_cholmod.nmethods = 1;
m_cholmod.method[0].ordering = CHOLMOD_NATURAL;
m_cholmod.postorder = 0;
}
m_cholmod.final_ll = 1;
m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
m_status = (m_status & ~SupernodalFactorIsDirty) | MatrixLIsDirty;
}
template<typename MatrixType>
inline const typename SparseCholesky<MatrixType>::CholMatrixType&
SparseCholesky<MatrixType,Cholmod>::matrixL() const
{
if (m_status & MatrixLIsDirty)
{
ei_assert(!(m_status & SupernodalFactorIsDirty));
cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*cmRes);
free(cmRes);
m_status = (m_status & ~MatrixLIsDirty);
}
return m_matrix;
}
template<typename MatrixType>
template<typename Derived>
void SparseCholesky<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
if (m_status & MatrixLIsDirty)
{
// ei_assert(!(m_status & SupernodalFactorIsDirty));
// taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
matrixL();
}
// else
{
Base::solveInPlace(b);
}
c.final_ll = 1;
cholmod_factor *L = cholmod_analyze(&A, &c);
cholmod_factorize(&A, L, &c);
cholmod_sparse* cmRes = cholmod_factor_to_sparse(L, &c);
m_matrix = CholMatrixType::Map(*cmRes);
free(cmRes);
cholmod_free_factor(&L, &c);
cholmod_finish(&c);
}
#endif // EIGEN_CHOLMODSUPPORT_H

67
Eigen/src/Sparse/SparseCholesky.h
View File

@ -25,12 +25,25 @@
#ifndef EIGEN_SPARSECHOLESKY_H
#define EIGEN_SPARSECHOLESKY_H
enum SparseBackend {
DefaultBackend,
Taucs,
Cholmod,
SuperLU
};
enum {
CholFull = 0x0, // full is the default
CholPartial = 0x1,
CompleteFactorization = 0x0, // full is the default
IncompleteFactorization = 0x1,
MemoryEfficient = 0x2,
SupernodalMultifrontal = 0x4,
SupernodalLeftLooking = 0x8,
/*
CholUseEigen = 0x0, // Eigen's impl is the default
CholUseTaucs = 0x2,
CholUseCholmod = 0x4,
CholUseCholmod = 0x4*/
};
/** \ingroup Sparse_Module
@ -43,23 +56,22 @@ enum {
*
* \sa class Cholesky, class CholeskyWithoutSquareRoot
*/
template<typename MatrixType> class SparseCholesky
template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
{
private:
protected:
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
typedef SparseMatrix<Scalar,Lower> CholMatrixType;
enum {
PacketSize = ei_packet_traits<Scalar>::size,
AlignmentMask = int(PacketSize)-1
SupernodalFactorIsDirty = 0x10000,
MatrixLIsDirty = 0x20000,
};
public:
SparseCholesky(const MatrixType& matrix, int flags = 0)
: m_matrix(matrix.rows(), matrix.cols()), m_flags(flags)
: m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
{
compute(matrix);
}
@ -69,17 +81,11 @@ template<typename MatrixType> class SparseCholesky
/** \returns true if the matrix is positive definite */
inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
// TODO impl the solver
// template<typename Derived>
// typename Derived::Eval solve(const MatrixBase<Derived> &b) const;
template<typename Derived>
void solveInPlace(MatrixBase<Derived> &b) const;
void compute(const MatrixType& matrix);
protected:
void computeUsingEigen(const MatrixType& matrix);
void computeUsingTaucs(const MatrixType& matrix);
void computeUsingCholmod(const MatrixType& matrix);
protected:
/** \internal
* Used to compute and store L
@ -87,24 +93,14 @@ template<typename MatrixType> class SparseCholesky
*/
CholMatrixType m_matrix;
int m_flags;
mutable int m_status;
bool m_isPositiveDefinite;
};
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
*/
template<typename MatrixType>
void SparseCholesky<MatrixType>::compute(const MatrixType& a)
{
if (m_flags&CholUseTaucs)
computeUsingTaucs(a);
else if (m_flags&CholUseCholmod)
computeUsingCholmod(a);
else
computeUsingEigen(a);
}
template<typename MatrixType>
void SparseCholesky<MatrixType>::computeUsingEigen(const MatrixType& a)
template<typename MatrixType, int Backend>
void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
{
assert(a.rows()==a.cols());
const int size = a.rows();
@ -173,4 +169,15 @@ void SparseCholesky<MatrixType>::computeUsingEigen(const MatrixType& a)
m_matrix.endFill();
}
template<typename MatrixType, int Backend>
template<typename Derived>
void SparseCholesky<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
m_matrix.solveTriangularInPlace(b);
m_matrix.adjoint().solveTriangularInPlace(b);
}
#endif // EIGEN_BASICSPARSECHOLESKY_H

110
Eigen/src/Sparse/TaucsSupport.h
View File

@ -79,12 +79,112 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& tau
}
template<typename MatrixType>
void SparseCholesky<MatrixType>::computeUsingTaucs(const MatrixType& a)
class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
{
taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
m_matrix = CholMatrixType::Map(*taucsRes);
free(taucsRes);
protected:
typedef SparseCholesky<MatrixType> Base;
using Base::Scalar;
using Base::RealScalar;
using Base::MatrixLIsDirty;
using Base::SupernodalFactorIsDirty;
using Base::m_flags;
using Base::m_matrix;
using Base::m_status;
public:
SparseCholesky(const MatrixType& matrix, int flags = 0)
: Base(matrix, flags), m_taucsSupernodalFactor(0)
{
compute(matrix);
}
~SparseCholesky()
{
if (m_taucsSupernodalFactor)
taucs_supernodal_factor_free(m_taucsSupernodalFactor);
}
inline const typename Base::CholMatrixType& matrixL(void) const;
template<typename Derived>
void solveInPlace(MatrixBase<Derived> &b) const;
void compute(const MatrixType& matrix);
protected:
void* m_taucsSupernodalFactor;
};
template<typename MatrixType>
void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
{
if (m_taucsSupernodalFactor)
{
taucs_supernodal_factor_free(m_taucsSupernodalFactor);
m_taucsSupernodalFactor = 0;
}
if (m_flags & IncompleteFactorization)
{
taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
m_matrix = Base::CholMatrixType::Map(*taucsRes);
free(taucsRes);
m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
| IncompleteFactorization
| SupernodalFactorIsDirty;
}
else
{
taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
if ( (m_flags & SupernodalLeftLooking)
|| ((!(m_flags & SupernodalMultifrontal)) && (m_flags & MemoryEfficient)) )
{
m_taucsSupernodalFactor = taucs_ccs_factor_llt_ll(&taucsMatA);
}
else
{
// use the faster Multifrontal routine
m_taucsSupernodalFactor = taucs_ccs_factor_llt_ll(&taucsMatA);
}
m_status = (m_status & ~IncompleteFactorization) | CompleteFactorization | MatrixLIsDirty;
}
}
template<typename MatrixType>
inline const typename SparseCholesky<MatrixType>::CholMatrixType&
SparseCholesky<MatrixType,Taucs>::matrixL() const
{
if (m_status & MatrixLIsDirty)
{
ei_assert(!(m_status & SupernodalFactorIsDirty));
taucs_ccs_matrix* taucsL = taucs_supernodal_factor_to_ccs(m_taucsSupernodalFactor);
const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*taucsL);
free(taucsL);
m_status = (m_status & ~MatrixLIsDirty);
}
return m_matrix;
}
template<typename MatrixType>
template<typename Derived>
void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
if (m_status & MatrixLIsDirty)
{
// ei_assert(!(m_status & SupernodalFactorIsDirty));
// taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
matrixL();
}
// else
{
Base::solveInPlace(b);
}
}
#endif // EIGEN_TAUCSSUPPORT_H

10
Eigen/src/Sparse/TriangularSolver.h
View File

@ -25,16 +25,6 @@
#ifndef EIGEN_SPARSETRIANGULARSOLVER_H
#define EIGEN_SPARSETRIANGULARSOLVER_H
// template<typename Lhs, typename Rhs,
// int TriangularPart = (int(Lhs::Flags) & LowerTriangularBit)
// ? Lower
// : (int(Lhs::Flags) & UpperTriangularBit)
// ? Upper
// : -1,
// int StorageOrder = int(Lhs::Flags) & RowMajorBit ? RowMajor : ColMajor
// >
// struct ei_sparse_trisolve_selector;
// forward substitution, row-major
template<typename Lhs, typename Rhs>
struct ei_solve_triangular_selector<Lhs,Rhs,Lower,RowMajor|IsSparse>