From c7303a876f8a051855ec9866a4aad1c694b2140d Mon Sep 17 00:00:00 2001 From: Hauke Heibel Date: Fri, 22 May 2009 15:58:20 +0200 Subject: [PATCH] Oops, here the actual LLT and LDLT patch. --- Eigen/src/Cholesky/LDLT.h | 42 ++++++++++++++++++++++++++++++++++----- Eigen/src/Cholesky/LLT.h | 25 +++++++++++++++++++++-- test/cholesky.cpp | 39 ++++++++++++++++++------------------ 3 files changed, 80 insertions(+), 26 deletions(-) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 5cd6f1bb1..b393afde8 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -62,16 +62,29 @@ template class LDLT typedef Matrix IntColVectorType; typedef Matrix 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 matrixL(void) const { return m_matrix; } + inline Part 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 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 vectorD(void) const { return m_matrix.diagonal(); } + inline Diagonal 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 bool solve(const MatrixBase &b, MatrixBase *result) const; @@ -110,6 +136,7 @@ template 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::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::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 bool LDLT ::solve(const MatrixBase &b, MatrixBase *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 template bool LDLT::solveInPlace(MatrixBase &bAndX) const { + ei_assert(m_isInitialized && "LDLT is not initialized."); const int size = m_matrix.rows(); ei_assert(size == bAndX.rows()); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 14ec51a68..b8ff3a3de 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -65,14 +65,27 @@ template 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 matrixL(void) const { return m_matrix; } + inline Part matrixL(void) const + { + ei_assert(m_isInitialized && "LLT is not initialized."); + return m_matrix; + } template bool solve(const MatrixBase &b, MatrixBase *result) const; @@ -88,6 +101,7 @@ template 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::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::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 template bool LLT::solve(const MatrixBase &b, MatrixBase *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 template bool LLT::solveInPlace(MatrixBase &bAndX) const { + ei_assert(m_isInitialized && "LLT is not initialized."); const int size = m_matrix.rows(); ei_assert(size==bAndX.rows()); matrixL().solveTriangularInPlace(bAndX); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index a49625e90..9bcf228d2 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -112,29 +112,25 @@ template void cholesky(const MatrixType& m) } -template -void doSomeRankPreservingOperations(Eigen::MatrixBase& m) +template 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(-1,1); - int i = Eigen::ei_random(0,m.rows()-1); // i is a random row number - int j; - do { - j = Eigen::ei_random(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(0,m.cols()-1); // i is a random column number - do { - j = Eigen::ei_random(0,m.cols()-1); - } while (i==j); // j is another one (must be different) - m.col(i) += d * m.col(j); - } + LLT llt; + VERIFY_RAISES_ASSERT(llt.matrixL()) + VERIFY_RAISES_ASSERT(llt.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(llt.solveInPlace(&tmp)) + + LDLT 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() ); + CALL_SUBTEST( cholesky_verify_assert() ); + CALL_SUBTEST( cholesky_verify_assert() ); + CALL_SUBTEST( cholesky_verify_assert() ); }