Fix bug #286: Infinite loop in JacobiSVD with denormals

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() );
 }