Fix conjugate-gradient for very small rhs

This commit is contained in:
Gael Guennebaud 2018-09-13 23:53:28 +02:00
parent 7f3b17e403
commit 1141bcf794

View File

@ -50,7 +50,8 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
tol_error = 0;
return;
}
RealScalar threshold = tol*tol*rhsNorm2;
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
RealScalar threshold = numext::maxi(tol*tol*rhsNorm2,considerAsZero);
RealScalar residualNorm2 = residual.squaredNorm();
if (residualNorm2 < threshold)
{
@ -58,7 +59,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
tol_error = sqrt(residualNorm2 / rhsNorm2);
return;
}
VectorType p(n);
p = precond.solve(residual); // initial search direction