From 654eab3bd6dc842c668d8979eb36602cef929165 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 19 Nov 2013 11:53:48 +0100
Subject: [PATCH] Add scaling in JacobiSVD to avoid overflows

---
 Eigen/src/SVD/JacobiSVD.h | 11 +++++++++--
 test/jacobisvd.cpp        | 13 +++++++++++--
 2 files changed, 20 insertions(+), 4 deletions(-)

diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index bf1213656..eee31ca97 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -411,8 +411,8 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
 
 template<typename MatrixType, typename RealScalar, typename Index>
 void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
-                            JacobiRotation<RealScalar> *j_left,
-                            JacobiRotation<RealScalar> *j_right)
+                         JacobiRotation<RealScalar> *j_left,
+                         JacobiRotation<RealScalar> *j_right)
 {
   using std::sqrt;
   Matrix<RealScalar,2,2> m;
@@ -831,6 +831,11 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
     if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
     if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
   }
+  
+  // Scaling factor to reducover/under-flows
+  RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff();
+  if(scale==RealScalar(0)) scale = RealScalar(1);
+  m_workMatrix /= scale;
 
   /*** step 2. The main Jacobi SVD iteration. ***/
 
@@ -879,6 +884,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
     m_singularValues.coeffRef(i) = a;
     if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
   }
+  
+  m_singularValues *= scale;
 
   /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
 
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index a0d8562c1..d441a6eca 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -285,7 +285,7 @@ void jacobisvd_inf_nan()
 
 // Regression test for bug 286: JacobiSVD loops indefinitely with some
 // matrices containing denormal numbers.
-void jacobisvd_bug286()
+void jacobisvd_underoverflow()
 {
 #if defined __INTEL_COMPILER
 // shut up warning #239: floating point underflow
@@ -300,6 +300,15 @@ void jacobisvd_bug286()
 #endif
   JacobiSVD<Matrix2d> svd;
   svd.compute(M); // just check we don't loop indefinitely
+  
+  // Check for overflow:
+  Matrix3d M3;
+  M3 << 4.4331978442502944e+307, -5.8585363752028680e+307,  6.4527017443412964e+307,
+        3.7841695601406358e+307,  2.4331702789740617e+306, -3.5235707140272905e+307,
+       -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307;
+
+  JacobiSVD<Matrix3d> svd3;
+  svd3.compute(M3); // just check we don't loop indefinitely
 }
 
 void jacobisvd_preallocate()
@@ -398,5 +407,5 @@ void test_jacobisvd()
   CALL_SUBTEST_9( jacobisvd_preallocate() );
 
   // Regression check for bug 286
-  CALL_SUBTEST_2( jacobisvd_bug286() );
+  CALL_SUBTEST_2( jacobisvd_underoverflow() );
 }