make jacobi SVD more robust after experimenting with very nasty matrices...

it turns out to be better to repeat the jacobi steps on a given (p,q) pair until it
is diagonal to machine precision, before going to the next (p,q) pair. it's also
an optimization as experiments show that in a majority of cases this allows to find out
that the (p,q) pair is already diagonal to machine precision.

1 files changed, 25 insertions, 16 deletions

diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h
index 18416e312..18c3db980 100644
--- a/Eigen/src/SVD/JacobiSquareSVD.h
+++ b/Eigen/src/SVD/JacobiSquareSVD.h
@@ -102,14 +102,16 @@ void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType&
   if(ComputeU) m_matrixU = MatrixUType::Identity(size,size);
   if(ComputeV) m_matrixV = MatrixUType::Identity(size,size);
   m_singularValues.resize(size);
-  while(true)
+
+sweep_again:
+  for(int p = 1; p < size; ++p)
   {
-    bool finished = true;
-    for(int p = 1; p < size; ++p)
+    for(int q = 0; q < p; ++q)
     {
-      for(int q = 0; q < p; ++q)
-      {
-        Scalar c, s;
+      Scalar c, s;
+      while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
+          > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*machine_epsilon<Scalar>())
+      {  
         if(work_matrix.makeJacobiForAtA(p,q,&c,&s))
         {
           work_matrix.applyJacobiOnTheRight(p,q,c,s);
@@ -117,9 +119,13 @@ void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType&
         }
         if(work_matrix.makeJacobiForAAt(p,q,&c,&s))
         {
+          Scalar x = ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(p,q));
+          Scalar y = ei_conj(work_matrix.coeff(q,p))*work_matrix.coeff(p,p) + ei_conj(work_matrix.coeff(q,q))*work_matrix.coeff(p,q);
+          Scalar z = ei_abs2(work_matrix.coeff(q,p)) + ei_abs2(work_matrix.coeff(q,q));
           work_matrix.applyJacobiOnTheLeft(p,q,c,s);
           if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s);
-          if(std::max(ei_abs(work_matrix.coeff(p,q)), ei_abs(work_matrix.coeff(q,p))) > std::max(ei_abs(work_matrix.coeff(q,q)), ei_abs(work_matrix.coeff(p,p)) ))
+          if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
+           > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q))) )
           {
             work_matrix.row(p).swap(work_matrix.row(q));
             if(ComputeU) m_matrixU.col(p).swap(m_matrixU.col(q));
@@ -127,15 +133,18 @@ void JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType&
         }
       }
     }
-    RealScalar biggest = work_matrix.diagonal().cwise().abs().maxCoeff();
-    for(int p = 0; p < size; ++p)
-    {
-      for(int q = 0; q < size; ++q)
-      {
-        if(p!=q && ei_abs(work_matrix.coeff(p,q)) > biggest * machine_epsilon<Scalar>()) finished = false;
-      }
-    }
-    if(finished) break;
+  }
+  
+  RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff();
+  RealScalar maxAllowedOffDiag = biggestOnDiag * machine_epsilon<Scalar>();
+  for(int p = 0; p < size; ++p)
+  {
+    for(int q = 0; q < p; ++q)
+      if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
+        goto sweep_again;
+    for(int q = p+1; q < size; ++q)
+      if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag)
+        goto sweep_again;
   }
 
   m_singularValues = work_matrix.diagonal().cwise().abs();