Merge branch 'feature/blocksolver-missing-diagonal' into 'master'

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.*

See merge request !11

dune/solvers/iterationsteps/blockgssteps.hh

@@ -40,16 +40,23 @@ 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;
+    }
+
+    using Ignore = std::bitset<V::block_type::dimension>;
 
-    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 != nullptr ? *diag : Block(0.0), std::move(ri),
+                        ignore != nullptr ? (*ignore)[i] : Ignore(0));
   };
 
   if (direction != Direction::BACKWARD)