From 5f42104e0a0d11d58a40f666bec38111c9d2fca1 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Tue, 23 Feb 2010 18:10:31 -0500
Subject: [PATCH] really fix the LDLt, at the expense of letting
 isPositiveDefinite() always return true (it was fundamentally broken anyway,
 especially as in 2.0 we don't even pivot at all). also, fix compilation

---
 Eigen/src/Cholesky/LDLT.h | 14 ++++----------
 test/cholesky.cpp         |  2 ++
 2 files changed, 6 insertions(+), 10 deletions(-)

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 8a466e1c7..cdc2da0fb 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -96,8 +96,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
   assert(a.rows()==a.cols());
   const int size = a.rows();
   m_matrix.resize(size, size);
-  m_isPositiveDefinite = true;
-  const RealScalar eps = precision<Scalar>() * a.cwise().abs().maxCoeff();
+  m_isPositiveDefinite = true; // always true. This decomposition is not rank-revealing anyway.
 
   if (size<=1)
   {
@@ -121,12 +120,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
     RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
     m_matrix.coeffRef(j,j) = tmp;
 
-    if (tmp < eps)
-    {
-      m_isPositiveDefinite = false;
-      return;
-    }
-
     int endSize = size-j-1;
     if (endSize>0)
     {
@@ -136,7 +129,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
       m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
                                    - _temporary.end(endSize).transpose();
 
-      m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
+      if(tmp != RealScalar(0))
+        m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
     }
   }
 }
@@ -192,7 +186,7 @@ template<typename Derived>
 inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::ldlt() const
 {
-  return LDLT<PlainObject>(derived());
+  return LDLT<PlainMatrixType>(derived());
 }
 
 #endif // EIGEN_LDLT_H
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 2e3353d21..4f26baa91 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -100,6 +100,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     VERIFY_IS_APPROX(symm * matX, matB);
   }
 
+#if 0 // cholesky is not rank-revealing anyway
   // test isPositiveDefinite on non definite matrix
   if (rows>4)
   {
@@ -109,6 +110,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     LDLT<SquareMatrixType> cholnosqrt(symm);
     VERIFY(!cholnosqrt.isPositiveDefinite());
   }
+#endif
 }
 
 void test_cholesky()