diff options
-rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 10 | ||||
-rw-r--r-- | test/jacobisvd.cpp | 16 |
2 files changed, 21 insertions, 5 deletions
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 5f6139998..e5e682489 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -590,6 +590,9 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig // only worsening the precision of U and V as we accumulate more rotations const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); + // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) + const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits<RealScalar>::denorm_min(); + /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ if(!internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix) @@ -617,10 +620,11 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig { // if this 2x2 sub-matrix is not diagonal already... // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't - // keep us iterating forever. + // keep us iterating forever. Similarly, small denormal numbers are considered zero. using std::max; - if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) - > (max)(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision) + RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)), + internal::abs(m_workMatrix.coeff(q,q)))); + if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold) { finished = false; diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 14c965133..f44e66c94 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -244,6 +244,17 @@ void jacobisvd_inf_nan() svd.compute(m, ComputeFullU | ComputeFullV); } +// Regression test for bug 286: JacobiSVD loops indefinitely with some +// matrices containing denormal numbers. +void jacobisvd_bug286() +{ + Matrix2d M; + M << -7.90884e-313, -4.94e-324, + 0, 5.60844e-313; + JacobiSVD<Matrix2d> svd; + svd.compute(M); // just check we don't loop indefinitely +} + void jacobisvd_preallocate() { Vector3f v(3.f, 2.f, 1.f); @@ -280,8 +291,6 @@ void jacobisvd_preallocate() internal::set_is_malloc_allowed(false); svd2.compute(m, ComputeFullU|ComputeFullV); internal::set_is_malloc_allowed(true); - - } void test_jacobisvd() @@ -336,4 +345,7 @@ void test_jacobisvd() // Check that preallocation avoids subsequent mallocs CALL_SUBTEST_9( jacobisvd_preallocate() ); + + // Regression check for bug 286 + CALL_SUBTEST_2( jacobisvd_bug286() ); } |