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

Oops, here the actual LLT and LDLT patch.

Browse Source
This commit is contained in:
Hauke Heibel 2009-05-22 15:58:20 +02:00
parent 0523b64fe9
commit c7303a876f
3 changed files with 80 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
Eigen/src/Cholesky/LDLT.h
View File

@ -62,16 +62,29 @@ template<typename MatrixType> class LDLT
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef Matrix<int, 1, MatrixType::RowsAtCompileTime> IntRowVectorType;
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LDLT::compute(const MatrixType&).
*/
LDLT() : m_matrix(), m_p(), m_transpositions(), m_isInitialized(false) {}
LDLT(const MatrixType& matrix)
: m_matrix(matrix.rows(), matrix.cols()),
m_p(matrix.rows()),
m_transpositions(matrix.rows())
m_transpositions(matrix.rows()),
m_isInitialized(false)
{
compute(matrix);
}
/** \returns the lower triangular matrix L */
inline Part<MatrixType, UnitLowerTriangular> matrixL(void) const { return m_matrix; }
inline Part<MatrixType, UnitLowerTriangular> matrixL(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix;
}
/** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed,
* representing the P permutation i.e. the permutation of the rows. For its precise meaning,
@ -79,17 +92,30 @@ template<typename MatrixType> class LDLT
*/
inline const IntColVectorType& permutationP() const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_p;
}
/** \returns the coefficients of the diagonal matrix D */
inline Diagonal<MatrixType,0> vectorD(void) const { return m_matrix.diagonal(); }
inline Diagonal<MatrixType,0> vectorD(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix.diagonal();
}
/** \returns true if the matrix is positive (semidefinite) */
inline bool isPositive(void) const { return m_sign == 1; }
inline bool isPositive(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == 1;
}
/** \returns true if the matrix is negative (semidefinite) */
inline bool isNegative(void) const { return m_sign == -1; }
inline bool isNegative(void) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == -1;
}
template<typename RhsDerived, typename ResDerived>
bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
@ -110,6 +136,7 @@ template<typename MatrixType> class LDLT
IntColVectorType m_p;
IntColVectorType m_transpositions;
int m_sign;
bool m_isInitialized;
};
/** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
@ -126,6 +153,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
m_p.setZero();
m_transpositions.setZero();
m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1;
m_isInitialized = true;
return;
}
@ -205,6 +233,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
for(int k = size-1; k >= 0; --k) {
std::swap(m_p.coeffRef(k), m_p.coeffRef(m_transpositions.coeff(k)));
}
m_isInitialized = true;
}
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
@ -222,6 +252,7 @@ template<typename RhsDerived, typename ResDerived>
bool LDLT<MatrixType>
::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
const int size = m_matrix.rows();
ei_assert(size==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
*result = b;
@ -243,6 +274,7 @@ template<typename MatrixType>
template<typename Derived>
bool LDLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
{
ei_assert(m_isInitialized && "LDLT is not initialized.");
const int size = m_matrix.rows();
ei_assert(size == bAndX.rows());

25
Eigen/src/Cholesky/LLT.h
View File

@ -65,14 +65,27 @@ template<typename MatrixType> class LLT
public:
/**
* \brief Default Constructor.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via LLT::compute(const MatrixType&).
*/
LLT() : m_matrix(), m_isInitialized(false) {}
LLT(const MatrixType& matrix)
: m_matrix(matrix.rows(), matrix.cols())
: m_matrix(matrix.rows(), matrix.cols()),
m_isInitialized(false)
{
compute(matrix);
}
/** \returns the lower triangular matrix L */
inline Part<MatrixType, LowerTriangular> matrixL(void) const { return m_matrix; }
inline Part<MatrixType, LowerTriangular> matrixL(void) const
{
ei_assert(m_isInitialized && "LLT is not initialized.");
return m_matrix;
}
template<typename RhsDerived, typename ResDerived>
bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
@ -88,6 +101,7 @@ template<typename MatrixType> class LLT
* The strict upper part is not used and even not initialized.
*/
MatrixType m_matrix;
bool m_isInitialized;
};
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
@ -109,7 +123,10 @@ void LLT<MatrixType>::compute(const MatrixType& a)
x = ei_real(a.coeff(0,0));
m_matrix.coeffRef(0,0) = ei_sqrt(x);
if(size==1)
{
m_isInitialized = true;
return;
}
m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
for (int j = 1; j < size; ++j)
{
@ -130,6 +147,8 @@ void LLT<MatrixType>::compute(const MatrixType& a)
- m_matrix.col(j).end(endSize) ) / x;
}
}
m_isInitialized = true;
}
/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
@ -149,6 +168,7 @@ template<typename MatrixType>
template<typename RhsDerived, typename ResDerived>
bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
{
ei_assert(m_isInitialized && "LLT is not initialized.");
const int size = m_matrix.rows();
ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
return solveInPlace((*result) = b);
@ -169,6 +189,7 @@ template<typename MatrixType>
template<typename Derived>
bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
{
ei_assert(m_isInitialized && "LLT is not initialized.");
const int size = m_matrix.rows();
ei_assert(size==bAndX.rows());
matrixL().solveTriangularInPlace(bAndX);

39
test/cholesky.cpp
View File

@ -112,29 +112,25 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
}
template<typename Derived>
void doSomeRankPreservingOperations(Eigen::MatrixBase<Derived>& m)
template<typename MatrixType> void cholesky_verify_assert()
{
typedef typename Derived::RealScalar RealScalar;
for(int a = 0; a < 3*(m.rows()+m.cols()); a++)
{
RealScalar d = Eigen::ei_random<RealScalar>(-1,1);
int i = Eigen::ei_random<int>(0,m.rows()-1); // i is a random row number
int j;
do {
j = Eigen::ei_random<int>(0,m.rows()-1);
} while (i==j); // j is another one (must be different)
m.row(i) += d * m.row(j);
MatrixType tmp;
i = Eigen::ei_random<int>(0,m.cols()-1); // i is a random column number
do {
j = Eigen::ei_random<int>(0,m.cols()-1);
} while (i==j); // j is another one (must be different)
m.col(i) += d * m.col(j);
}
LLT<MatrixType> llt;
VERIFY_RAISES_ASSERT(llt.matrixL())
VERIFY_RAISES_ASSERT(llt.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp))
LDLT<MatrixType> ldlt;
VERIFY_RAISES_ASSERT(ldlt.matrixL())
VERIFY_RAISES_ASSERT(ldlt.permutationP())
VERIFY_RAISES_ASSERT(ldlt.vectorD())
VERIFY_RAISES_ASSERT(ldlt.isPositive())
VERIFY_RAISES_ASSERT(ldlt.isNegative())
VERIFY_RAISES_ASSERT(ldlt.solve(tmp,&tmp))
VERIFY_RAISES_ASSERT(ldlt.solveInPlace(&tmp))
}
void test_cholesky()
{
for(int i = 0; i < g_repeat; i++) {
@ -147,4 +143,9 @@ void test_cholesky()
CALL_SUBTEST( cholesky(MatrixXd(17,17)) );
CALL_SUBTEST( cholesky(MatrixXf(200,200)) );
}
CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() );
CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() );
CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() );
CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() );
}