From 77b658f8da62673f60e0a3c8b331e33d6d64d021 Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@googlemail.com>
Date: Fri, 25 Sep 2020 14:29:03 +0200
Subject: [PATCH] [test.la] fix undefined behavior
---
dune/xt/la/solver/istl/schurcomplement.hh | 19 ++++++++++++++-----
dune/xt/test/la/solver.tpl | 4 ++++
2 files changed, 18 insertions(+), 5 deletions(-)
diff --git a/dune/xt/la/solver/istl/schurcomplement.hh b/dune/xt/la/solver/istl/schurcomplement.hh
index 06c71fb8b..38806e292 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 7bca7c8ca..6cf370c3d 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.;
--
GitLab