LDLT:

* fix bug thanks to Ben Goodrich: we were terminating at the wrong place, leaving some matrix coefficients with wrong values.
* don't use Higham's formula here: we're not trying to be rank-revealing.

1 files changed, 9 insertions, 12 deletions

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 4d3149d42..708b02375 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -202,11 +202,8 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
     {
       // The biggest overall is the point of reference to which further diagonals
       // are compared; if any diagonal is negligible compared
-      // to the largest overall, the algorithm bails.  This cutoff is suggested
-      // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
-      // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
-      // Algorithms" page 217, also by Higham.
-      cutoff = ei_abs(NumTraits<Scalar>::epsilon() * RealScalar(size) * biggest_in_corner);
+      // to the largest overall, the algorithm bails.
+      cutoff = ei_abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
 
       m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
     }
@@ -235,13 +232,6 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
                                                .dot(m_matrix.col(j).head(j)));
     m_matrix.coeffRef(j,j) = Djj;
 
-    // Finish early if the matrix is not full rank.
-    if(ei_abs(Djj) < cutoff)
-    {
-      for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
-      break;
-    }
-
     int endSize = size - j - 1;
     if (endSize > 0) {
       _temporary.tail(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
@@ -250,6 +240,13 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
       m_matrix.row(j).tail(endSize) = m_matrix.row(j).tail(endSize).conjugate()
                                    - _temporary.tail(endSize).transpose();
 
+      // Finish early if the matrix is not full rank.
+      if(ei_abs(Djj) < cutoff)
+      {
+        for(int i = j; i < size; i++) m_transpositions.coeffRef(i) = i;
+        break;
+      }
+
       m_matrix.col(j).tail(endSize) = m_matrix.row(j).tail(endSize) / Djj;
     }
   }