* remove EIGEN_DONT_INLINE that harm performance for small sizes

* normalize left Jacobi rotations to avoid having to swap rows
* set precision to 2*machine_epsilon instead of machine_epsilon, we lose 1 bit of precision
  but gain between 10% and 100% speed, plus reduce the risk that some day we hit a bad matrix
  where it's impossible to approach machine precision

1 files changed, 4 insertions, 11 deletions

diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h
index 18c3db980..82133f7be 100644
--- a/Eigen/src/SVD/JacobiSquareSVD.h
+++ b/Eigen/src/SVD/JacobiSquareSVD.h
@@ -102,6 +102,7 @@ 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);
+  const RealScalar precision = 2 * machine_epsilon<Scalar>();
 
 sweep_again:
   for(int p = 1; p < size; ++p)
@@ -110,7 +111,7 @@ sweep_again:
     {
       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>())
+          > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision)
       {  
         if(work_matrix.makeJacobiForAtA(p,q,&c,&s))
         {
@@ -119,24 +120,16 @@ sweep_again:
         }
         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));
+          ei_normalizeJacobi(&c, &s, work_matrix.coeff(p,p), work_matrix.coeff(q,p)),
           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(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));
-          }
         }
       }
     }
   }
   
   RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff();
-  RealScalar maxAllowedOffDiag = biggestOnDiag * machine_epsilon<Scalar>();
+  RealScalar maxAllowedOffDiag = biggestOnDiag * precision;
   for(int p = 0; p < size; ++p)
   {
     for(int q = 0; q < p; ++q)