Fix bug #791: infinite loop in JacobiSVD in the presence of NaN.

(grafted from d6236d3b26
)
This commit is contained in:
Gael Guennebaud 2014-09-10 11:54:20 +02:00
parent c67a7148c4
commit e0ab58d815
2 changed files with 14 additions and 5 deletions

View File

@ -861,7 +861,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
using std::max;
RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)),
abs(m_workMatrix.coeff(q,q))));
if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold)
// We compare both values to threshold instead of calling max to be robust to NaN (See bug 791)
if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold)
{
finished = false;

View File

@ -321,16 +321,23 @@ void jacobisvd_inf_nan()
VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
Scalar some_nan = zero<Scalar>() / zero<Scalar>();
VERIFY(some_nan != some_nan);
svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
VERIFY(nan != nan);
svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV);
MatrixType m = MatrixType::Zero(10,10);
m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
svd.compute(m, ComputeFullU | ComputeFullV);
m = MatrixType::Zero(10,10);
m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan;
svd.compute(m, ComputeFullU | ComputeFullV);
// regression test for bug 791
m.resize(3,3);
m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5,
0, -0.5, 0,
nan, 0, 0;
svd.compute(m, ComputeFullU | ComputeFullV);
}
@ -434,6 +441,7 @@ void test_jacobisvd()
// Test on inf/nan matrix
CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
CALL_SUBTEST_10( jacobisvd_inf_nan<MatrixXd>() );
}
CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));