diff --git a/dune/istl/istlexception.hh b/dune/istl/istlexception.hh
index d72f6feb770b6c56455bafc15024c1f4d830a89c..4a1da2d3f4473f6cdc9c8ff1d487f3cd4592992f 100644
--- a/dune/istl/istlexception.hh
+++ b/dune/istl/istlexception.hh
@@ -33,6 +33,15 @@ namespace Dune {
: public BCRSMatrixError
{};
+ //! Thrown when a solver aborts due to some problem.
+ /**
+ * Problems that may cause the solver to abort include a NaN detected during
+ * the convergence check (which may be caused by invalid input data), or
+ * breakdown conditions (which can happen e.g. in BiCGSTABSolver or
+ * RestartedGMResSolver).
+ */
+ class SolverAbort : public ISTLError {};
+
/** @} end documentation */
} // end namespace
diff --git a/dune/istl/solver.hh b/dune/istl/solver.hh
index 8898d9adfc26acc7958d6eb208a5ba3b2a502267..d010fe09b679c85d0607e7eaaba713aadb6f87be 100644
--- a/dune/istl/solver.hh
+++ b/dune/istl/solver.hh
@@ -95,6 +95,9 @@ namespace Dune
\param x The left hand side to store the result in.
\param b The right hand side
\param res Object to store the statistics about applying the operator.
+
+ \throw SolverAbort When the solver detects a problem and cannot
+ continue
*/
virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0;
@@ -107,6 +110,9 @@ namespace Dune
\param b The right hand side
\param reduction The minimum defect reduction to achieve.
\param res Object to store the statistics about applying the operator.
+
+ \throw SolverAbort When the solver detects a problem and cannot
+ continue
*/
virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0;
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index 7f9e1c1b787735d18d738eab582e678e12d15e3b..320354f0187e9e6e8c928caf403cac7636d97c6c 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -19,6 +19,7 @@
#include "solver.hh"
#include "preconditioner.hh"
#include <dune/common/deprecated.hh>
+#include <dune/common/exceptions.hh>
#include <dune/common/timer.hh>
#include <dune/common/ftraits.hh>
#include <dune/common/typetraits.hh>
@@ -410,9 +411,17 @@ namespace Dune {
\brief Apply inverse operator.
\copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
+
+ \note Currently, the CGSolver aborts when a NaN or infinite defect is
+ detected. However, -ffinite-math-only (implied by -ffast-math)
+ can inhibit a result from becoming NaN that really should be NaN.
+ E.g. numeric_limits<double>::quiet_NaN()*0.0==0.0 with gcc-5.3
+ -ffast-math.
*/
virtual void apply (X& x, X& b, InverseOperatorResult& res)
{
+ using std::isfinite;
+
res.clear(); // clear solver statistics
Timer watch; // start a timer
_prec.pre(x,b); // prepare preconditioner
@@ -422,6 +431,16 @@ namespace Dune {
X q(x); // a temporary vector
real_type def0 = _sp.norm(b); // compute norm
+
+ if (!isfinite(def0)) // check for inf or NaN
+ {
+ if (_verbose>0)
+ std::cout << "=== CGSolver: abort due to infinite or NaN initial defect"
+ << std::endl;
+ DUNE_THROW(SolverAbort, "CGSolver: initial defect=" << def0
+ << " is infinite or NaN");
+ }
+
if (def0<1E-30) // convergence check
{
res.converged = true;
@@ -472,6 +491,15 @@ namespace Dune {
this->printOutput(std::cout,real_type(i),defnew,def);
def = defnew; // update norm
+ if (!isfinite(def)) // check for inf or NaN
+ {
+ if (_verbose>0)
+ std::cout << "=== CGSolver: abort due to infinite or NaN defect"
+ << std::endl;
+ DUNE_THROW(SolverAbort,
+ "CGSolver: defect=" << def << " is infinite or NaN");
+ }
+
if (def<def0*_reduction || def<1E-30) // convergence check
{
res.converged = true;
@@ -513,6 +541,12 @@ namespace Dune {
\brief Apply inverse operator with given reduction factor.
\copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
+
+ \note Currently, the CGSolver aborts when a NaN or infinite defect is
+ detected. However, -ffinite-math-only (implied by -ffast-math)
+ can inhibit a result from becoming NaN that really should be NaN.
+ E.g. numeric_limits<double>::quiet_NaN()*0.0==0.0 with gcc-5.3
+ -ffast-math.
*/
virtual void apply (X& x, X& b, double reduction,
InverseOperatorResult& res)
@@ -583,6 +617,8 @@ namespace Dune {
\brief Apply inverse operator.
\copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
+
+ \note Currently, the BiCGSTABSolver aborts when it detects a breakdown.
*/
virtual void apply (X& x, X& b, InverseOperatorResult& res)
{
@@ -661,11 +697,11 @@ namespace Dune {
// look if breakdown occurred
if (abs(rho) <= EPSILON)
- DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho "
+ DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - rho "
<< rho << " <= EPSILON " << EPSILON
<< " after " << it << " iterations");
if (abs(omega) <= EPSILON)
- DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega "
+ DUNE_THROW(SolverAbort,"breakdown in BiCGSTAB - omega "
<< omega << " <= EPSILON " << EPSILON
<< " after " << it << " iterations");
@@ -691,7 +727,7 @@ namespace Dune {
h = _sp.dot(rt,v);
if (abs(h) < EPSILON)
- DUNE_THROW(ISTLError,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
+ DUNE_THROW(SolverAbort,"abs(h) < EPSILON in BiCGSTAB - abs(h) "
<< abs(h) << " < EPSILON " << EPSILON
<< " after " << it << " iterations");
@@ -785,6 +821,8 @@ namespace Dune {
\brief Apply inverse operator with given reduction factor.
\copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
+
+ \note Currently, the BiCGSTABSolver aborts when it detects a breakdown.
*/
virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res)
{
@@ -1191,7 +1229,14 @@ namespace Dune {
"P and S must have the same category!");
}
- //! \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
+ /*!
+ \brief Apply inverse operator.
+
+ \copydoc InverseOperator::apply(X&,Y&,InverseOperatorResult&)
+
+ \note Currently, the RestartedGMResSolver aborts when it detects a
+ breakdown.
+ */
virtual void apply (X& x, Y& b, InverseOperatorResult& res)
{
apply(x,b,_reduction,res);
@@ -1201,6 +1246,9 @@ namespace Dune {
\brief Apply inverse operator.
\copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
+
+ \note Currently, the RestartedGMResSolver aborts when it detects a
+ breakdown.
*/
virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res)
{
@@ -1277,7 +1325,7 @@ namespace Dune {
}
H[i+1][i] = _sp.norm(w);
if(abs(H[i+1][i]) < EPSILON)
- DUNE_THROW(ISTLError,
+ DUNE_THROW(SolverAbort,
"breakdown in GMRes - |w| == 0.0 after " << j << " iterations");
// normalize new vector
diff --git a/dune/istl/test/CMakeLists.txt b/dune/istl/test/CMakeLists.txt
index 0a998ea23bd7bcd00f55bf6b1897237235b0fc59..33f7d295b0b6e46082c9634ca8d4e62068bbd38a 100644
--- a/dune/istl/test/CMakeLists.txt
+++ b/dune/istl/test/CMakeLists.txt
@@ -45,6 +45,9 @@ dune_add_test(SOURCES inverseoperator2prectest.cc)
dune_add_test(SOURCES scaledidmatrixtest.cc)
+dune_add_test(SOURCES solveraborttest.cc
+ SKIP_ON_77)
+
# Pardiso tests
dune_add_test(SOURCES test_pardiso.cc
SKIP_ON_77)
diff --git a/dune/istl/test/solveraborttest.cc b/dune/istl/test/solveraborttest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..57ae8f9eb4f5d39191fbd3ed83c59fd288c36620
--- /dev/null
+++ b/dune/istl/test/solveraborttest.cc
@@ -0,0 +1,113 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include "config.h"
+
+#include <iostream>
+#include <limits>
+#include <string>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/istl/operators.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solver.hh>
+#include <dune/istl/solvers.hh>
+
+template<class Solver, class Vector>
+void checkSolverAbort(int &status, const std::string &name,
+ Solver &solver, Vector &x, Vector &b)
+{
+ try {
+ Dune::InverseOperatorResult res;
+ solver.apply(x, b, res);
+ std::cout << "Error: " << name << "::apply() did not abort, which is "
+ << "unexpected.\n"
+ << "converged = " << std::boolalpha << res.converged << "\n"
+ << "iterations = " << res.iterations << "\n"
+ << "reduction = " << res.reduction << "\n"
+ << "conv_rate = " << res.conv_rate << "\n"
+ << "elapsed = " << res.elapsed << "\n"
+ << "solution = {" << x << "}" << std::endl;
+ status = 1; // FAIL
+ }
+ catch(const Dune::SolverAbort &e)
+ {
+ std::cout << name << "::apply() aborted, as expected" << std::endl
+ << "Abort message was: " << e << std::endl;
+ if(status == 77) status = 0; // PASS
+ }
+ catch(const std::exception &e)
+ {
+ std::cout << "Error: " << name << "::apply() aborted with an exception "
+ << "not derived from Dune::SolverAbort, which is unexpected.\n"
+ << "e.what(): " << e.what() << std::endl;
+ status = 1; // FAIL
+ }
+ catch(...)
+ {
+ std::cout << "Error: " << name << "::apply() aborted with an exception "
+ << "not derived from Dune::SolverAbort, which is unexpected.\n"
+ << "In addition, the exception is not derived from "
+ << "std::exception, so there is no further information, sorry."
+ << std::endl;
+ status = 1; // FAIL
+ }
+}
+
+int main()
+{
+
+ int status = 77;
+
+ // How verbose the solvers should be. Use 2 (maximum verbosity) by default,
+ // this will include all information in the logs, and for the casual user of
+ // the unit tests ctest will hide the output anyway.
+ int verbose = 2;
+
+ { // CGSolver
+ std::cout << "Checking CGSolver with an unsolvable system...\n"
+ << "Expecting SolverAbort with a NaN defect" << std::endl;
+
+ using Matrix = Dune::FieldMatrix<double, 2, 2>;
+ using Vector = Dune::FieldVector<double, 2>;
+
+ Matrix matrix = { { 1, 1 },
+ { 1, 1 } };
+ Vector b = { 1, 2 };
+ Vector x = { 0, 0 };
+
+ Dune::MatrixAdapter<Matrix, Vector, Vector> op(matrix);
+ Dune::Richardson<Vector, Vector> richardson;
+ Dune::CGSolver<Vector> solver(op, richardson, 1e-10, 5000, verbose);
+
+ checkSolverAbort(status, "CGSolver", solver, x, b);
+ }
+
+ { // BiCGSTABSolver
+ std::cout << "Checking BiCGSTABSolver with an unsolvable system...\n"
+ << "Expecting abs(h) < EPSILON" << std::endl;
+
+ using Matrix = Dune::FieldMatrix<double, 2, 2>;
+ using Vector = Dune::FieldVector<double, 2>;
+
+ Matrix matrix = { { 1, 1 },
+ { 1, 1 } };
+ Vector b = { 1, 2 };
+ Vector x = { 0, 0 };
+
+ Dune::MatrixAdapter<Matrix, Vector, Vector> op(matrix);
+ Dune::Richardson<Vector, Vector> richardson;
+ Dune::BiCGSTABSolver<Vector> solver(op, richardson, 1e-10, 5000, verbose);
+
+ checkSolverAbort(status, "BiCGSTABSolver", solver, x, b);
+ }
+
+ // TODO:
+ // - trigger "breakdown in BiCGSTAB - rho"
+ // - trigger "breakdown in BiCGSTAB - omega"
+ // - trigger "breakdown in GMRes"
+
+ return status;
+}