Fix: Cope with omitted diagonal blocks in outer GS loop.

The former implementation had two shortcomings:

1. It assumed the diagonal block to exist. This is not necessarily the case, e.g., in the case of a BCRSMatrix where a non-retrievable block within the dimensions of the matrix is equivalent to it being zero. Nevertheless we should be able to cope with this case, in particular because semi-definite matrices should be covered by this approach due to the inherent regularization.
2. Access to the diagonal element was performed via the random access operator which is relatively costly.

The new implementation tackles both problems. In case the diagonal block does not exist, a temporary block is constructed and passed to the local solver instead.

On a side note: A more generic approach to creating and initializing the temporary block might be needed for more complex nested matrix schemes.

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);
 
     std::bitset<V::block_type::dimension> ignore_i(0);
     if (ignore != nullptr)
       ignore_i = (*ignore)[i];
+
+    // Compute residual
+    auto ri = b[i];
+    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;
+    }
+
     // 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)
@@ -172,9 +180,10 @@ namespace LocalSolverFromLinearSolver {
  */
 template <class LinearSolver>
 auto truncateSymmetrically(LinearSolver&& linearSolver) {
-  return [linearSolver = std::move(linearSolver) ](
-      const auto& m, const auto& b, const auto& ignore) {
-    using Return = typename std::result_of<LinearSolver(decltype(m), decltype(b))>::type;
+  return [linearSolver = std::move(linearSolver)](const auto& m, const auto& b,
+                                                  const auto& ignore) {
+    using Return =
+        typename std::result_of<LinearSolver(decltype(m), decltype(b))>::type;
     if (ignore.all())
       return Return(0);
 
@@ -270,7 +279,8 @@ auto cg(size_t maxIter = LinearSolvers::defaultCgMaxIter,
     auto cgSolver = std::bind(
         LinearSolvers::cg<std::decay_t<decltype(m)>, std::decay_t<decltype(b)>>,
         _1, _2, maxIter, tol);
-    auto cgTruncatedSolver = LocalSolverFromLinearSolver::truncateSymmetrically(cgSolver);
+    auto cgTruncatedSolver =
+        LocalSolverFromLinearSolver::truncateSymmetrically(cgSolver);
     auto cgTruncatedRegularizedSolver =
         LocalSolverRegularizer::diagRegularize(r, cgTruncatedSolver);
     return cgTruncatedRegularizedSolver(m, b, ignore);