diff --git a/dune/solvers/iterationsteps/blockgssteps.hh b/dune/solvers/iterationsteps/blockgssteps.hh
index 2bbdb0ff2f16ef24b60e386325ff906e9cb3ba20..7ed5a385c6cab2cd83b3c18d0472320ab4cbe19c 100644
--- a/dune/solvers/iterationsteps/blockgssteps.hh
+++ b/dune/solvers/iterationsteps/blockgssteps.hh
@@ -7,6 +7,7 @@
 #include <dune/common/parametertree.hh>
 
 #include <dune/matrix-vector/ldlt.hh>
+#include <dune/matrix-vector/algorithm.hh>
 
 #include <dune/solvers/common/defaultbitvector.hh>
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
@@ -45,12 +46,11 @@ void linearStep(const M& m, V& x, const V& b, const BitVector* ignore,
     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);
+    MatrixVector::sparseRangeFor(row_i, [&](const auto& mij, auto j) {
+      mij.mmv(x[j], ri);
       if (j == i)
-        diag = &*cIt;
-    }
+        diag = &mij;
+      });
 
     using Ignore = std::bitset<V::block_type::dimension>;
 
@@ -150,12 +150,12 @@ VBlock gs(const MBlock& m, const VBlock& b, double tol = defaultGsTol) {
     if (std::abs(mii) <= tol)
       continue;
     x[i] = b[i];
-    const auto& end = mi.end();
-    for (auto it = mi.begin(); it != end; ++it) {
-      auto j = it.index();
+
+    MatrixVector::sparseRangeFor(mi, [&](const auto& mij, auto j) {
       if (j != i)
-        x[i] -= (*it) * x[j];
-    }
+        x[i] -= mij * x[j];
+    });
+
     x[i] /= mii;
   }
   return x;