mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source
This commit is contained in:
Gael Guennebaud 2014-07-08 10:04:27 +02:00
parent 7fa83e7374
commit 0dfb73d46a
2 changed files with 7 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
Eigen/src/Cholesky/LDLT.h
View File

@ -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)

2
test/cholesky.cpp
View File

@ -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);