diff options
| author |  Gael Guennebaud <g.gael@free.fr> | 2014-07-08 10:04:27 +0200 |
|---|---|---|
| committer | Gael Guennebaud <g.gael@free.fr> | 2014-07-08 10:04:27 +0200 |
| commit | 0dfb73d46ada6a9749a24a946186c8a8c2472dc5 (patch) | |
| tree | b9c53d969232bfbcc70a87438fbffd07ab4953cf /Eigen/src/Cholesky | |
| parent | 7fa83e7374b38e9900e77cf2b29b54919f65d1ef (diff) | |
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
Diffstat (limited to 'Eigen/src/Cholesky')
| -rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 12 |
1 files changed, 6 insertions, 6 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) |