From ee9f3e34b0c54d40c9613cce74233d86a1a96182 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
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<typename MatrixType,template <typename,int> class CholType> void test_chol_update(const MatrixType& symm)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+
+  MatrixType symmLo = symm.template triangularView<Lower>();
+  MatrixType symmUp = symm.template triangularView<Upper>();
+  MatrixType symmCpy = symm;
+
+  CholType<MatrixType,Lower> chollo(symmLo);
+  CholType<MatrixType,Upper> cholup(symmUp);
+
+  for (int k=0; k<10; ++k)
+  {
+    VectorType vec = VectorType::Random(symm.rows());
+    RealScalar sigma = internal::random<RealScalar>();
+    symmCpy += sigma * vec * vec.adjoint();
+
+    // we are doing some downdates, so it might be the case that the matrix is not SPD anymore
+    CholType<MatrixType,Lower> 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<typename MatrixType> void cholesky(const MatrixType& m)
 {
   typedef typename MatrixType::Index Index;
@@ -155,41 +187,9 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
   m2.noalias() -= symmLo.template selfadjointView<Lower>().llt().solve(matB);
   VERIFY_IS_APPROX(m2, m1 - symmLo.template selfadjointView<Lower>().llt().solve(matB));
 
-  // LLT update/downdate
-  {
-    MatrixType symmLo = symm.template triangularView<Lower>();
-    MatrixType symmUp = symm.template triangularView<Upper>();
-    
-    VectorType vec = VectorType::Random(rows);
-
-    MatrixType symmCpy = symm + vec * vec.adjoint();
-
-    LLT<MatrixType,Lower> chollo(symmLo);
-    chollo.rankUpdate(vec);
-    VERIFY_IS_APPROX(symmCpy,  chollo.reconstructedMatrix());
-
-    LLT<MatrixType,Upper> cholup(symmUp);
-    cholup.rankUpdate(vec);
-    VERIFY_IS_APPROX(symmCpy,  cholup.reconstructedMatrix());
-  }
-  
-  // LDLT update/downdate
-  {
-    MatrixType symmLo = symm.template triangularView<Lower>();
-    MatrixType symmUp = symm.template triangularView<Upper>();
-    
-    VectorType vec = VectorType::Random(rows);
-
-    MatrixType symmCpy = symm + vec * vec.adjoint();
-
-    LDLT<MatrixType,Lower> chollo(symmLo);
-    chollo.rankUpdate(vec);
-    VERIFY_IS_APPROX(symmCpy,  chollo.reconstructedMatrix());
-
-    LDLT<MatrixType,Upper> cholup(symmUp);
-    cholup.rankUpdate(vec);
-    VERIFY_IS_APPROX(symmCpy,  cholup.reconstructedMatrix());
-  }
+  // update/downdate
+  CALL_SUBTEST(( test_chol_update<SquareMatrixType,LLT>(symm)  ));
+  CALL_SUBTEST(( test_chol_update<SquareMatrixType,LDLT>(symm) ));
 }
 
 template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
-- 
cgit v1.2.3