From bfeb29421d4231965a46d0c770800565a156aad1 Mon Sep 17 00:00:00 2001
From: Max Kahnt <max.kahnt@fu-berlin.de>
Date: Fri, 18 Nov 2016 14:06:04 +0100
Subject: [PATCH] Fix: Cope with omitted diagonal blocks in outer GS loop.
---
dune/solvers/iterationsteps/blockgssteps.hh | 14 +++++++++++---
1 file changed, 11 insertions(+), 3 deletions(-)
diff --git a/dune/solvers/iterationsteps/blockgssteps.hh b/dune/solvers/iterationsteps/blockgssteps.hh
index 01d17c64..bd68b4cb 100644
--- a/dune/solvers/iterationsteps/blockgssteps.hh
+++ b/dune/solvers/iterationsteps/blockgssteps.hh
@@ -40,16 +40,24 @@ void linearStep(const M& m, V& x, const V& b, const BitVector* ignore,
// Note: move-capture requires C++14
auto blockStep = [&, localSolver = std::move(localSolver) ](size_t i) {
const auto& row_i = m[i];
+
// Compute residual
auto ri = b[i];
- for (auto cIt = row_i.begin(); cIt != row_i.end(); ++cIt)
- cIt->mmv(x[cIt.index()], ri);
+ using Block = typename M::block_type;
+ const Block* diag = nullptr;
+ for (auto cIt = row_i.begin(); cIt != row_i.end(); ++cIt) {
+ size_t j = cIt.index();
+ cIt->mmv(x[j], ri);
+ if (j == i)
+ diag = &*cIt;
+ }
std::bitset<V::block_type::dimension> ignore_i(0);
if (ignore != nullptr)
ignore_i = (*ignore)[i];
+
// Update iterate with correction
- x[i] += localSolver(row_i[i], std::move(ri), ignore_i);
+ x[i] += localSolver(diag ? *diag : Block(0.0), std::move(ri), ignore_i);
};
if (direction != Direction::BACKWARD)
--
GitLab