diff options
Diffstat (limited to 'Eigen/src')
| -rw-r--r-- | Eigen/src/Jacobi/Jacobi.h | 2 | ||||
| -rw-r--r-- | Eigen/src/SVD/JacobiSquareSVD.h | 41 |
2 files changed, 26 insertions, 17 deletions
diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 5866ac44f..48195c1ec 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -50,7 +50,7 @@ void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s template<typename Scalar> bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) { - if(ei_abs(y) < ei_abs(z-x) * 0.5 * machine_epsilon<Scalar>()) + if(ei_abs(y) <= ei_abs(z-x) * 0.5 * machine_epsilon<Scalar>()) { *c = Scalar(1); *s = Scalar(0); 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(); |