devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Big API change in Cholesky module:

Browse Source 
* rename Cholesky to LLT
 * rename CholeskyWithoutSquareRoot to LDLT
 * rename MatrixBase::cholesky() to llt()
 * rename MatrixBase::choleskyNoSqrt() to ldlt()
 * make {LLT,LDLT}::solve() API consistent with other modules

Note that we are going to keep a source compatibility untill the next beta release.
E.g., the "old" Cholesky* classes, etc are still available for some time.
To be clear, Eigen beta2 should be (hopefully) source compatible with beta1,
and so beta2 will contain all the deprecated API of beta1. Those features marked
as deprecated will be removed in beta3 (or in the final 2.0 if there is no beta 3 !).

Also includes various updated in sparse Cholesky.
This commit is contained in:
Gael Guennebaud
2008-10-13 15:53:27 +00:00
parent e2bd8623f8
commit 765219aa51
23 changed files with 608 additions and 183 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
Eigen/src/Sparse/AmbiVector.h
View File

@@ -28,7 +28,7 @@
/** \internal
* Hybrid sparse/dense vector class designed for intensive read-write operations.
*
* See BasicSparseCholesky and SparseProduct for usage examples.
* See BasicSparseLLT and SparseProduct for usage examples.
*/
template<typename _Scalar> class AmbiVector
{
@@ -269,12 +269,14 @@ class AmbiVector<_Scalar>::Iterator
/** Default constructor
* \param vec the vector on which we iterate
* \param nonZeroReferenceValue reference value used to prune zero coefficients.
* In practice, the coefficient are compared to \a nonZeroReferenceValue * precision<Scalar>().
* \param epsilon the minimal value used to prune zero coefficients.
* In practice, all coefficients having a magnitude smaller than \a epsilon
* are skipped.
*/
Iterator(const AmbiVector& vec, RealScalar nonZeroReferenceValue = RealScalar(0.1)) : m_vector(vec)
Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*precision<RealScalar>())
: m_vector(vec)
{
m_epsilon = nonZeroReferenceValue * precision<Scalar>();
m_epsilon = epsilon;
m_isDense = m_vector.m_mode==IsDense;
if (m_isDense)
{

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

@@ -99,10 +99,10 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(cholmod_sparse& cm)
}
template<typename MatrixType>
class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
{
protected:
typedef SparseCholesky<MatrixType> Base;
typedef SparseLLT<MatrixType> Base;
using Base::Scalar;
using Base::RealScalar;
using Base::MatrixLIsDirty;
@@ -113,14 +113,20 @@ class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
public:
SparseCholesky(const MatrixType& matrix, int flags = 0)
SparseLLT(int flags = 0)
: Base(flags), m_cholmodFactor(0)
{
cholmod_start(&m_cholmod);
}
SparseLLT(const MatrixType& matrix, int flags = 0)
: Base(matrix, flags), m_cholmodFactor(0)
{
cholmod_start(&m_cholmod);
compute(matrix);
}
~SparseCholesky()
~SparseLLT()
{
if (m_cholmodFactor)
cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
@@ -140,7 +146,7 @@ class SparseCholesky<MatrixType,Cholmod> : public SparseCholesky<MatrixType>
};
template<typename MatrixType>
void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a)
{
if (m_cholmodFactor)
{
@@ -169,8 +175,8 @@ void SparseCholesky<MatrixType,Cholmod>::compute(const MatrixType& a)
}
template<typename MatrixType>
inline const typename SparseCholesky<MatrixType>::CholMatrixType&
SparseCholesky<MatrixType,Cholmod>::matrixL() const
inline const typename SparseLLT<MatrixType>::CholMatrixType&
SparseLLT<MatrixType,Cholmod>::matrixL() const
{
if (m_status & MatrixLIsDirty)
{
@@ -187,7 +193,7 @@ SparseCholesky<MatrixType,Cholmod>::matrixL() const
template<typename MatrixType>
template<typename Derived>
void SparseCholesky<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());

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

@@ -38,25 +38,19 @@ enum {
MemoryEfficient = 0x2,
SupernodalMultifrontal = 0x4,
SupernodalLeftLooking = 0x8
/*
CholUseEigen = 0x0, // Eigen's impl is the default
CholUseTaucs = 0x2,
CholUseCholmod = 0x4*/
};
/** \ingroup Sparse_Module
*
* \class SparseCholesky
* \class SparseLLT
*
* \brief Standard Cholesky decomposition of a matrix and associated features
* \brief Standard LLT decomposition of a matrix and associated features
*
* \param MatrixType the type of the matrix of which we are computing the Cholesky decomposition
* \param MatrixType the type of the matrix of which we are computing the LLT decomposition
*
* \sa class Cholesky, class CholeskyWithoutSquareRoot
* \sa class LLT, class LDLT
*/
template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
template<typename MatrixType, int Backend = DefaultBackend> class SparseLLT
{
protected:
typedef typename MatrixType::Scalar Scalar;
@@ -70,66 +64,66 @@ template<typename MatrixType, int Backend = DefaultBackend> class SparseCholesky
public:
SparseCholesky(const MatrixType& matrix, int flags = 0)
SparseLLT(int flags = 0)
: m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
}
SparseLLT(const MatrixType& matrix, int flags = 0)
: m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0)
{
m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>();
compute(matrix);
}
inline const CholMatrixType& matrixL(void) const { return m_matrix; }
void setPrecision(RealScalar v) { m_precision = v; }
RealScalar precision() const { return m_precision; }
/** \returns true if the matrix is positive definite */
inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
void setFlags(int f) { m_flags = f; }
int flags() const { return m_flags; }
void compute(const MatrixType& matrix);
inline const CholMatrixType& matrixL(void) const { return m_matrix; }
template<typename Derived>
void solveInPlace(MatrixBase<Derived> &b) const;
void compute(const MatrixType& matrix);
/** \returns true if the factorization succeeded */
inline bool succeeded(void) const { return m_succeeded; }
protected:
/** \internal
* Used to compute and store L
* The strict upper part is not used and even not initialized.
*/
CholMatrixType m_matrix;
RealScalar m_precision;
int m_flags;
mutable int m_status;
bool m_isPositiveDefinite;
bool m_succeeded;
};
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
/** Computes / recomputes the LLT decomposition A = LL^* = U^*U of \a matrix
*/
template<typename MatrixType, int Backend>
void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
void SparseLLT<MatrixType,Backend>::compute(const MatrixType& a)
{
assert(a.rows()==a.cols());
const int size = a.rows();
m_matrix.resize(size, size);
const RealScalar eps = ei_sqrt(precision<Scalar>());
// const RealScalar eps = ei_sqrt(precision<Scalar>());
// allocate a temporary vector for accumulations
AmbiVector<Scalar> tempVector(size);
RealScalar density = a.nonZeros()/RealScalar(size*size);
// TODO estimate the number of nnz
m_matrix.startFill(a.nonZeros()*2);
for (int j = 0; j < size; ++j)
{
// std::cout << j << "\n";
Scalar x = ei_real(a.coeff(j,j));
int endSize = size-j-1;
// TODO estimate the number of non zero entries
// float ratioLhs = float(lhs.nonZeros())/float(lhs.rows()*lhs.cols());
// float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
// float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
// let's do a more accurate determination of the nnz ratio for the current column j of res
//float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
// FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
// float ratioColRes = ratioRes;
// if (ratioColRes>0.1)
// tempVector.init(IsSparse);
tempVector.init(IsDense);
// TODO better estimate the density !
tempVector.init(density>0.001? IsDense : IsSparse);
tempVector.setBounds(j+1,size);
tempVector.setZero();
// init with current matrix a
@@ -161,7 +155,7 @@ void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
RealScalar rx = ei_sqrt(ei_real(x));
m_matrix.fill(j,j) = rx;
Scalar y = Scalar(1)/rx;
for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
for (typename AmbiVector<Scalar>::Iterator it(tempVector, m_precision*rx); it; ++it)
{
m_matrix.fill(it.index(), j) = it.value() * y;
}
@@ -171,7 +165,7 @@ void SparseCholesky<MatrixType,Backend>::compute(const MatrixType& a)
template<typename MatrixType, int Backend>
template<typename Derived>
void SparseCholesky<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
void SparseLLT<MatrixType, Backend>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());

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

@@ -79,10 +79,10 @@ SparseMatrix<Scalar,Flags> SparseMatrix<Scalar,Flags>::Map(taucs_ccs_matrix& tau
}
template<typename MatrixType>
class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
{
protected:
typedef SparseCholesky<MatrixType> Base;
typedef SparseLLT<MatrixType> Base;
using Base::Scalar;
using Base::RealScalar;
using Base::MatrixLIsDirty;
@@ -93,13 +93,18 @@ class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
public:
SparseCholesky(const MatrixType& matrix, int flags = 0)
SparseLLT(int flags = 0)
: Base(flags), m_taucsSupernodalFactor(0)
{
}
SparseLLT(const MatrixType& matrix, int flags = 0)
: Base(matrix, flags), m_taucsSupernodalFactor(0)
{
compute(matrix);
}
~SparseCholesky()
~SparseLLT()
{
if (m_taucsSupernodalFactor)
taucs_supernodal_factor_free(m_taucsSupernodalFactor);
@@ -117,7 +122,7 @@ class SparseCholesky<MatrixType,Taucs> : public SparseCholesky<MatrixType>
};
template<typename MatrixType>
void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
{
if (m_taucsSupernodalFactor)
{
@@ -128,7 +133,7 @@ void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
if (m_flags & IncompleteFactorization)
{
taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, 0, 0);
taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0);
m_matrix = Base::CholMatrixType::Map(*taucsRes);
free(taucsRes);
m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
@@ -153,8 +158,8 @@ void SparseCholesky<MatrixType,Taucs>::compute(const MatrixType& a)
}
template<typename MatrixType>
inline const typename SparseCholesky<MatrixType>::CholMatrixType&
SparseCholesky<MatrixType,Taucs>::matrixL() const
inline const typename SparseLLT<MatrixType>::CholMatrixType&
SparseLLT<MatrixType,Taucs>::matrixL() const
{
if (m_status & MatrixLIsDirty)
{
@@ -170,7 +175,7 @@ SparseCholesky<MatrixType,Taucs>::matrixL() const
template<typename MatrixType>
template<typename Derived>
void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
{
const int size = m_matrix.rows();
ei_assert(size==b.rows());
@@ -179,9 +184,9 @@ void SparseCholesky<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) cons
{
// ei_assert(!(m_status & SupernodalFactorIsDirty));
// taucs_supernodal_solve_llt(m_taucsSupernodalFactor,double* b);
matrixL();
//matrixL();
}
// else
else
{
Base::solveInPlace(b);
}

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

@@ -102,7 +102,7 @@ struct ei_solve_triangular_selector<Lhs,Rhs,Lower,ColMajor|IsSparse>
typename Lhs::InnerIterator it(lhs, i);
if(!(Lhs::Flags & UnitDiagBit))
{
std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
// std::cerr << it.value() << " ; " << it.index() << " == " << i << "\n";
ei_assert(it.index()==i);
other.coeffRef(i,col) /= it.value();
}