From 0dfb73d46ada6a9749a24a946186c8a8c2472dc5 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 8 Jul 2014 10:04:27 +0200
Subject: 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

---
 Eigen/src/Cholesky/LDLT.h | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

(limited to 'Eigen/src/Cholesky/LDLT.h')

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)
-- 
cgit v1.2.3