From 1b8cc9af43374e1adf8cd7a2c18d94dddb6080a6 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 27 Mar 2015 10:55:00 +0100 Subject: [PATCH] Slight numerical stability improvement in 2x2 svd --- Eigen/src/SVD/JacobiSVD.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index fcf01f518..6cef87f5e 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -425,12 +425,13 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, // If d!=0, then t/d cannot overflow because the magnitude of the // entries forming d are not too small compared to the ones forming t. RealScalar u = t / d; - rot1.s() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); - rot1.c() = rot1.s() * u; + RealScalar tmp = sqrt(RealScalar(1) + numext::abs2(u)); + rot1.s() = RealScalar(1) / tmp; + rot1.c() = u / tmp; } m.applyOnTheLeft(0,1,rot1); j_right->makeJacobi(m,0,1); - *j_left = rot1 * j_right->transpose(); + *j_left = rot1 * j_right->transpose(); } template @@ -680,6 +681,8 @@ JacobiSVD::compute(const MatrixType& matrix, unsig const RealScalar precision = RealScalar(2) * NumTraits::epsilon(); // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) + // FIXME What about considerering any denormal numbers as zero, using: + // const RealScalar considerAsZero = (std::numeric_limits::min)(); const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits::denorm_min(); // Scaling factor to reduce over/under-flows