diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-07-02 23:04:46 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-07-02 23:04:46 +0200 |
commit | 998455a5707e1312d3058672d6d337a20158d86e (patch) | |
tree | 6d69322ac18ed90460599050dab87aa3e8a5d80e /Eigen/src/Cholesky/LDLT.h | |
parent | 0a8e4712d1f3cfe8830783f59c6c9987ea34701c (diff) |
LDLT is not rank-revealing, so we should not attempt to use the biggest diagonal elements as thresholds.
Diffstat (limited to 'Eigen/src/Cholesky/LDLT.h')
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 27 |
1 files changed, 13 insertions, 14 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 1d84438e2..67a99a0af 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -258,23 +258,13 @@ template<> struct ldlt_inplace<Lower> return true; } - RealScalar cutoff(0), biggest_in_corner; - for (Index k = 0; k < size; ++k) { // Find largest diagonal element Index index_of_biggest_in_corner; - biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); + mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); index_of_biggest_in_corner += k; - if(k == 0) - { - // The biggest overall is the point of reference to which further diagonals - // are compared; if any diagonal is negligible compared - // to the largest overall, the algorithm bails. - cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner); - } - transpositions.coeffRef(k) = index_of_biggest_in_corner; if(k != index_of_biggest_in_corner) { @@ -311,7 +301,11 @@ template<> struct ldlt_inplace<Lower> A21.noalias() -= A20 * temp.head(k); } - if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) + // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot + // 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)); @@ -495,8 +489,13 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs> typedef typename LDLTType::Scalar Scalar; typedef typename LDLTType::RealScalar RealScalar; const Diagonal<const MatrixType> vectorD = dec().vectorD(); - RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(), - RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS + // 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()); + // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest + // diagonal element is not well justified and to numerical issues in some cases. + // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. + RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest(); for (Index i = 0; i < vectorD.size(); ++i) { if(abs(vectorD(i)) > tolerance) dst.row(i) /= vectorD(i); |