diff options
| -rw-r--r-- | Eigen/src/SVD/JacobiSVD.h | 11 | ||||
| -rw-r--r-- | 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() ); } |