Fix LDLT with semi-definite complex matrices: owing to round-off errors, the diagonal was not real. Also exploit the fact that the diagonal is real in the rest of LDLT

2 files changed, 7 insertions, 7 deletions

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index ef81030f8..6881e1ca8 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -295,7 +295,7 @@ template<> struct ldlt_inplace<Lower>
 
       if(k>0)
       {
-        temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint();
+        temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint();
         mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
         if(rs>0)
           A21.noalias() -= A20 * temp.head(k);
@@ -305,10 +305,10 @@ template<> struct ldlt_inplace<Lower>
       // was smaller than the cutoff value. However, soince LDLT is not rank-revealing
       // we should only make sure we do not introduce INF or NaN values.
       // LAPACK also uses 0 as the cutoff value.
-      if((rs>0) && (abs(mat.coeffRef(k,k)) > RealScalar(0)))
-        A21 /= mat.coeffRef(k,k);
-
       RealScalar realAkk = numext::real(mat.coeffRef(k,k));
+      if((rs>0) && (abs(realAkk) > RealScalar(0)))
+        A21 /= realAkk;
+
       if (sign == PositiveSemiDef) {
         if (realAkk < 0) sign = Indefinite;
       } else if (sign == NegativeSemiDef) {
@@ -487,7 +487,7 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
     EIGEN_USING_STD_MATH(max);
     typedef typename LDLTType::MatrixType MatrixType;
     typedef typename LDLTType::RealScalar RealScalar;
-    const Diagonal<const MatrixType> vectorD = dec().vectorD();
+    const typename Diagonal<const MatrixType>::RealReturnType vectorD(dec().vectorD());
     // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon
     // as motivated by LAPACK's xGELSS:
     // RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest());
@@ -552,7 +552,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
   // L^* P
   res = matrixU() * res;
   // D(L^*P)
-  res = vectorD().asDiagonal() * res;
+  res = vectorD().real().asDiagonal() * res;
   // L(DL^*P)
   res = matrixL() * res;
   // P^T (LDL^*P)
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 9aa7da7bb..1d147bd7a 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -209,7 +209,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
       VERIFY_IS_APPROX(A * vecX, vecB);
     }
     
-    // check matrices with wide spectrum
+    // check matrices with a wide spectrum
     if(rows>=3)
     {
       RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8);