diff options
author | Gael Guennebaud <g.gael@free.fr> | 2018-04-03 16:49:35 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2018-04-03 16:49:35 +0200 |
commit | c5b56f1fb27fb5b85eefef6b93dd71f4edb400db (patch) | |
tree | 76ee1f3bd6e379cc9029163c3a62df3d22b495f8 | |
parent | 8d0ffe36552aeeb5f46d9c652edc45c68e536cdd (diff) |
bug #1528: better use numeric_limits::min() instead of 1/highest() that with underflow.
-rw-r--r-- | Eigen/src/Cholesky/LDLT.h | 7 |
1 files changed, 4 insertions, 3 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 5be58377b..2dfeac333 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -569,13 +569,14 @@ void LDLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) cons // more precisely, use pseudo-inverse of D (see bug 241) using std::abs; const typename Diagonal<const MatrixType>::RealReturnType vecD(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: + // In some previous versions, tolerance was set to the max of 1/highest (or rather numeric_limits::min()) + // and the maximal diagonal entry * epsilon as motivated by LAPACK's xGELSS: // RealScalar tolerance = numext::maxi(vecD.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 leads to numerical issues in some cases. // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. - RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest(); + // Using numeric_limits::min() gives us more robustness to denormals. + RealScalar tolerance = (std::numeric_limits<RealScalar>::min)(); for (Index i = 0; i < vecD.size(); ++i) { |