From ee9f3e34b0c54d40c9613cce74233d86a1a96182 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 23 Jan 2012 17:28:23 +0100 Subject: LLT: improve rankUpdate to support downdates, LDLT: add the missing info() function, improve unit testing of rankUpdate() --- test/cholesky.cpp | 70 +++++++++++++++++++++++++++---------------------------- 1 file changed, 35 insertions(+), 35 deletions(-) (limited to 'test/cholesky.cpp') diff --git a/test/cholesky.cpp b/test/cholesky.cpp index d9806e5c3..1a1b2eeb5 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -41,6 +41,38 @@ static int nb_temporaries; VERIFY( (#XPR) && nb_temporaries==N ); \ } +template class CholType> void test_chol_update(const MatrixType& symm) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef Matrix VectorType; + + MatrixType symmLo = symm.template triangularView(); + MatrixType symmUp = symm.template triangularView(); + MatrixType symmCpy = symm; + + CholType chollo(symmLo); + CholType cholup(symmUp); + + for (int k=0; k<10; ++k) + { + VectorType vec = VectorType::Random(symm.rows()); + RealScalar sigma = internal::random(); + symmCpy += sigma * vec * vec.adjoint(); + + // we are doing some downdates, so it might be the case that the matrix is not SPD anymore + CholType chol(symmCpy); + if(chol.info()!=Success) + break; + + chollo.rankUpdate(vec, sigma); + VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); + + cholup.rankUpdate(vec, sigma); + VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); + } +} + template void cholesky(const MatrixType& m) { typedef typename MatrixType::Index Index; @@ -155,41 +187,9 @@ template void cholesky(const MatrixType& m) m2.noalias() -= symmLo.template selfadjointView().llt().solve(matB); VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView().llt().solve(matB)); - // LLT update/downdate - { - MatrixType symmLo = symm.template triangularView(); - MatrixType symmUp = symm.template triangularView(); - - VectorType vec = VectorType::Random(rows); - - MatrixType symmCpy = symm + vec * vec.adjoint(); - - LLT chollo(symmLo); - chollo.rankUpdate(vec); - VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); - - LLT cholup(symmUp); - cholup.rankUpdate(vec); - VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); - } - - // LDLT update/downdate - { - MatrixType symmLo = symm.template triangularView(); - MatrixType symmUp = symm.template triangularView(); - - VectorType vec = VectorType::Random(rows); - - MatrixType symmCpy = symm + vec * vec.adjoint(); - - LDLT chollo(symmLo); - chollo.rankUpdate(vec); - VERIFY_IS_APPROX(symmCpy, chollo.reconstructedMatrix()); - - LDLT cholup(symmUp); - cholup.rankUpdate(vec); - VERIFY_IS_APPROX(symmCpy, cholup.reconstructedMatrix()); - } + // update/downdate + CALL_SUBTEST(( test_chol_update(symm) )); + CALL_SUBTEST(( test_chol_update(symm) )); } template void cholesky_cplx(const MatrixType& m) -- cgit v1.2.3