mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

Add scaling in JacobiSVD to avoid overflows

Browse Source
This commit is contained in:
Gael Guennebaud 2013-11-19 11:53:48 +01:00
parent 5d1291a4de
commit 654eab3bd6
2 changed files with 20 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/SVD/JacobiSVD.h
View File

@ -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 ***/

13
test/jacobisvd.cpp
View File

@ -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() );
}