From e2aa58b6316ba869be039e49dce73f99647f4139 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 12 Jul 2016 17:21:03 +0200 Subject: [PATCH] Consider denormals as zero in makeJacobi and 2x2 SVD. This also fix serious issues with x387 for which values can be much smaller than the smallest denormal! --- Eigen/src/Jacobi/Jacobi.h | 5 +++-- Eigen/src/misc/RealSvd2x2.h | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 55de15e87..d25af8e90 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -85,7 +85,8 @@ bool JacobiRotation::makeJacobi(const RealScalar& x, const Scalar& y, co using std::sqrt; using std::abs; typedef typename NumTraits::Real RealScalar; - if(y == Scalar(0)) + RealScalar deno = RealScalar(2)*abs(y); + if(deno < (std::numeric_limits::min)()) { m_c = Scalar(1); m_s = Scalar(0); @@ -93,7 +94,7 @@ bool JacobiRotation::makeJacobi(const RealScalar& x, const Scalar& y, co } else { - RealScalar tau = (x-z)/(RealScalar(2)*abs(y)); + RealScalar tau = (x-z)/deno; RealScalar w = sqrt(numext::abs2(tau) + RealScalar(1)); RealScalar t; if(tau>RealScalar(0)) diff --git a/Eigen/src/misc/RealSvd2x2.h b/Eigen/src/misc/RealSvd2x2.h index dfaaa0b17..abb4d3c2f 100644 --- a/Eigen/src/misc/RealSvd2x2.h +++ b/Eigen/src/misc/RealSvd2x2.h @@ -28,7 +28,8 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation rot1; RealScalar t = m.coeff(0,0) + m.coeff(1,1); RealScalar d = m.coeff(1,0) - m.coeff(0,1); - if(d == RealScalar(0)) + + if(abs(d) < (std::numeric_limits::min)()) { rot1.s() = RealScalar(0); rot1.c() = RealScalar(1);