diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index ed60fc4b8..c52b7d1a6 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -311,7 +311,7 @@ template<> struct ldlt_inplace 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); @@ -321,10 +321,10 @@ template<> struct ldlt_inplace // 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) { @@ -504,8 +504,7 @@ struct solve_retval, Rhs> typedef typename LDLTType::MatrixType MatrixType; typedef typename LDLTType::Scalar Scalar; typedef typename LDLTType::RealScalar RealScalar; - const Diagonal vectorD = dec().vectorD(); - + const typename Diagonal::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::epsilon(),RealScalar(1) / NumTraits::highest()); @@ -571,7 +570,7 @@ MatrixType LDLT::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 65d8585e7..64bcbccc4 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -213,7 +213,7 @@ template 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::max_exponent10/8);