Failing convergence checks in solvers for good initial iterates
Recently I ran some user code that resulted in an attempt to solve a linear system using the BiCGSTABSolver
using dune-istl 2.6 with an initial iterate close to the solution. This resulted in a failure of the solver because the tolerance for convergence could not be reached.
After some digging I found the issue and noticed it still exists in the current master
version.
The line in question is:
solver.hh@473: _res.converged = (Simd::allTrue(def<_def0*_parent._reduction) || Simd::max(def)<1E-30); // convergence check
def
is the current residual, def0
is set to the initial residual.
If the initial residual is already close to the solution, def0*_parent._reduction < eps
with machine epsilon eps
and so the first condition cannot become true.
Likewise if for the type used eps > 1E-30
the second condition cannot become true either.
A possible solution for this is to set the magic number 1E-30
to a reasonable type-dependant value, whichever value that would be (std::numeric_limits<decltype(def)>::epsilon()
?).
Another way is to rethink the _reduction
parameter. Right now we only reduce the error by a fixed value, regardless of the quality of the initial error. An alternative would be to provide a target error size parameter and stop when this is reached (def < max_error).