GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-08 22:21:49 +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 bbe9e22d60
commit 08b0c08e5e
2 changed files with 7 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -311,7 +311,7 @@ template<> struct ldlt_inplace<Lower>
if(k>0) 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(); mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
if(rs>0) if(rs>0)
A21.noalias() -= A20 * temp.head(k); A21.noalias() -= A20 * temp.head(k);
@ -321,10 +321,10 @@ template<> struct ldlt_inplace<Lower>
// was smaller than the cutoff value. However, soince LDLT is not rank-revealing // 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. // we should only make sure we do not introduce INF or NaN values.
// LAPACK also uses 0 as the cutoff value. // 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)); RealScalar realAkk = numext::real(mat.coeffRef(k,k));
if((rs>0) && (abs(realAkk) > RealScalar(0)))
A21 /= realAkk;
if (sign == PositiveSemiDef) { if (sign == PositiveSemiDef) {
if (realAkk < 0) sign = Indefinite; if (realAkk < 0) sign = Indefinite;
} else if (sign == NegativeSemiDef) { } else if (sign == NegativeSemiDef) {
@ -504,8 +504,7 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
typedef typename LDLTType::MatrixType MatrixType; typedef typename LDLTType::MatrixType MatrixType;
typedef typename LDLTType::Scalar Scalar; typedef typename LDLTType::Scalar Scalar;
typedef typename LDLTType::RealScalar RealScalar; 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 // 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: // as motivated by LAPACK's xGELSS:
// RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); // RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest());
@ -571,7 +570,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
// L^* P // L^* P
res = matrixU() * res; res = matrixU() * res;
// D(L^*P) // D(L^*P)
res = vectorD().asDiagonal() * res; res = vectorD().real().asDiagonal() * res;
// L(DL^*P) // L(DL^*P)
res = matrixL() * res; res = matrixL() * res;
// P^T (LDL^*P) // P^T (LDL^*P)

2
test/cholesky.cpp
View File

@ -213,7 +213,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
VERIFY_IS_APPROX(A * vecX, vecB); VERIFY_IS_APPROX(A * vecX, vecB);
} }
// check matrices with wide spectrum // check matrices with a wide spectrum
if(rows>=3) if(rows>=3)
{ {
RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8); RealScalar s = (std::min)(16,std::numeric_limits<RealScalar>::max_exponent10/8);