From ef3e690a0c410e7309cb7a0ff60154943375db03 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 17 Jan 2011 11:09:03 +0100
Subject: [PATCH] return the index of the first non positive diagonal entry
 (more useful than simply true or false)

---
 Eigen/src/Cholesky/LLT.h | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index ceb532055..53a0543fe 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -184,11 +184,12 @@ template<int UpLo> struct llt_inplace;
 template<> struct llt_inplace<Lower>
 {
   template<typename MatrixType>
-  static bool unblocked(MatrixType& mat)
+  static typename MatrixType::Index unblocked(MatrixType& mat)
   {
+    typedef typename MatrixType::Index Index;
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
-    typedef typename MatrixType::Index Index;
+    
     eigen_assert(mat.rows()==mat.cols());
     const Index size = mat.rows();
     for(Index k = 0; k < size; ++k)
@@ -202,16 +203,16 @@ template<> struct llt_inplace<Lower>
       RealScalar x = real(mat.coeff(k,k));
       if (k>0) x -= A10.squaredNorm();
       if (x<=RealScalar(0))
-        return false;
+        return k;
       mat.coeffRef(k,k) = x = sqrt(x);
       if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
       if (rs>0) A21 *= RealScalar(1)/x;
     }
-    return true;
+    return -1;
   }
 
   template<typename MatrixType>
-  static bool blocked(MatrixType& m)
+  static typename MatrixType::Index blocked(MatrixType& m)
   {
     typedef typename MatrixType::Index Index;
     eigen_assert(m.rows()==m.cols());
@@ -235,24 +236,25 @@ template<> struct llt_inplace<Lower>
       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);
       Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs);
 
-      if(!unblocked(A11)) return false;
+      Index ret;
+      if((ret=unblocked(A11))>=0) return k+ret;
       if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21);
       if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck
     }
-    return true;
+    return -1;
   }
 };
 
 template<> struct llt_inplace<Upper>
 {
   template<typename MatrixType>
-  static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
+  static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat)
   {
     Transpose<MatrixType> matt(mat);
     return llt_inplace<Lower>::unblocked(matt);
   }
   template<typename MatrixType>
-  static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
+  static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat)
   {
     Transpose<MatrixType> matt(mat);
     return llt_inplace<Lower>::blocked(matt);
@@ -266,7 +268,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
   inline static MatrixL getL(const MatrixType& m) { return m; }
   inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
   static bool inplace_decomposition(MatrixType& m)
-  { return llt_inplace<Lower>::blocked(m); }
+  { return llt_inplace<Lower>::blocked(m)==-1; }
 };
 
 template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
@@ -276,7 +278,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
   inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
   inline static MatrixU getU(const MatrixType& m) { return m; }
   static bool inplace_decomposition(MatrixType& m)
-  { return llt_inplace<Upper>::blocked(m); }
+  { return llt_inplace<Upper>::blocked(m)==-1; }
 };
 
 } // end namespace internal