diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 88bc0688e..f08776bc6 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -350,7 +350,7 @@ template struct svd_precondition_2x2_block_to_be_real { typedef JacobiSVD SVD; - static void run(typename SVD::WorkMatrixType&, SVD&, Index, Index) {} + static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, const typename MatrixType::RealScalar&) { return true; } }; template @@ -359,7 +359,7 @@ struct svd_precondition_2x2_block_to_be_real typedef JacobiSVD SVD; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; - static void run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q) + static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, const typename MatrixType::RealScalar& precision) { using std::sqrt; Scalar z; @@ -368,7 +368,7 @@ struct svd_precondition_2x2_block_to_be_real if(n==0) { - // make sure firt column is zero (deflation) + // make sure first column is zero work_matrix.coeffRef(p,p) = work_matrix.coeffRef(q,p) = Scalar(0); if(work_matrix.coeff(p,q)!=Scalar(0)) { @@ -404,6 +404,12 @@ struct svd_precondition_2x2_block_to_be_real if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); } } + + const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits::denorm_min(); + RealScalar threshold = numext::maxi(considerAsZero, + precision * numext::maxi(abs(work_matrix.coeff(p,p)), abs(work_matrix.coeff(q,q)))); + // return true if we still have some work to do + return numext::abs(work_matrix(p,q)) > threshold || numext::abs(work_matrix(q,p)) > threshold; } }; @@ -735,18 +741,20 @@ JacobiSVD::compute(const MatrixType& matrix, unsig if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) { finished = false; - // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal - internal::svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q); - JacobiRotation j_left, j_right; - internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); + // the complex to real operation returns true is the updated 2x2 block is not already diagonal + if(internal::svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q, precision)) + { + JacobiRotation j_left, j_right; + internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); - // accumulate resulting Jacobi rotations - m_workMatrix.applyOnTheLeft(p,q,j_left); - if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); + // accumulate resulting Jacobi rotations + m_workMatrix.applyOnTheLeft(p,q,j_left); + if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); - m_workMatrix.applyOnTheRight(p,q,j_right); - if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); + m_workMatrix.applyOnTheRight(p,q,j_right); + if(computeV()) m_matrixV.applyOnTheRight(p,q,j_right); + } } } }