diff --git a/dune/xt/la/solver/istl/schurcomplement.hh b/dune/xt/la/solver/istl/schurcomplement.hh index 06c71fb8b964ed5bf7cdd0276a55de37fa675a8e..38806e29222de75950879b444d37212259ea739b 100644 --- a/dune/xt/la/solver/istl/schurcomplement.hh +++ b/dune/xt/la/solver/istl/schurcomplement.hh @@ -15,8 +15,10 @@ #include <dune/istl/operators.hh> #include <dune/istl/solvers.hh> -#include <dune/xt/common/exceptions.hh> #include <dune/xt/common/configuration.hh> +#include <dune/xt/common/exceptions.hh> +#include <dune/xt/common/math.hh> + #include <dune/xt/la/container/istl.hh> #include <dune/xt/la/solver.hh> @@ -82,10 +84,17 @@ public: auto& B1x = m_vec_1_; B1_.mv(x, B1x); // calculate A^{-1} B1 x - auto& AinvB1x = m_vec_2_; - A_inv_.apply(B1x, AinvB1x, solver_opts_); - // apply B2^T - B2_.mtv(AinvB1x, y); + if (XT::Common::is_zero(B1x.two_norm())) { + // if B1x is zero, y will also be zero + // this special case is here to avoid calling the dune-istl CG solver with an + // initial defect of 0 which causes undefined behavior (division-by-zero) + y *= 0.; + } else { + auto& AinvB1x = m_vec_2_; + A_inv_.apply(B1x, AinvB1x, solver_opts_); + // apply B2^T + B2_.mtv(AinvB1x, y); + } // calculate Cx auto& Cx = n_vec_1_; C_.mv(x, Cx); diff --git a/dune/xt/test/la/solver.tpl b/dune/xt/test/la/solver.tpl index 7bca7c8cab1db73b500b91bb46b414a7e32308ac..6cf370c3dda5cedc9d014ebe533161dcc2b6b28f 100644 --- a/dune/xt/test/la/solver.tpl +++ b/dune/xt/test/la/solver.tpl @@ -64,6 +64,10 @@ struct SolverTest_{{T_NAME}} : public ::testing::Test options.report(out, " "); // dynamic tests + // Matrix is the identity matrix, solution and rhs are initially filled with ones + // Multiply by 2 to avoid calling the istl iterative solver with an + // initial defect of 0 which results in undefined behavior (division-by-zero) + solution *= 2.; solver.apply(rhs, solution, type); EXPECT_TRUE(XT::Common::FloatCmp::eq(solution, rhs)); solution *= 0.;