diff --git a/dune/matrix-vector/axpy.hh b/dune/matrix-vector/axpy.hh
index 792dc034b71d522634ca3aab203f1b54abb6e29e..3a5e13038de1a1d0d41f5f29e497f890a8b746d0 100644
--- a/dune/matrix-vector/axpy.hh
+++ b/dune/matrix-vector/axpy.hh
@@ -8,6 +8,7 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/scaledidmatrix.hh>
 
+#include "algorithm.hh"
 #include "matrixtraits.hh"
 #include "scalartraits.hh"
 
@@ -30,18 +31,13 @@ namespace MatrixVector {
         class ADummy = A,
         std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
     static void addProduct(A& a, const B& b, const C& c) {
-      typename B::ConstRowIterator bi = b.begin();
-      typename B::ConstRowIterator bEnd = b.end();
-      for (; bi != bEnd; ++bi) {
-        typename B::ConstColIterator bik = bi->begin();
-        typename B::ConstColIterator biEnd = bi->end();
-        for (; bik != biEnd; ++bik) {
-          typename C::ConstColIterator ckj = c[bik.index()].begin();
-          typename C::ConstColIterator ckEnd = c[bik.index()].end();
-          for (; ckj != ckEnd; ++ckj)
-            a[bi.index()][ckj.index()] += (*bik) * (*ckj);
-        }
-      }
+      rangeForEach(b, [&](auto&& bi, auto&& i) {
+        rangeForEach(bi, [&](auto&& bik, auto&& k) {
+          rangeForEach(c[k], [&](auto&& ckj, auto&& j) {
+            Dune::MatrixVector::addProduct(a[i][j], bik, ckj);
+          });
+        });
+      });
     }
   };
 
@@ -61,18 +57,13 @@ namespace MatrixVector {
         class ADummy = A,
         std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
     static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
-      typename B::ConstRowIterator bi = b.begin();
-      typename B::ConstRowIterator bEnd = b.end();
-      for (; bi != bEnd; ++bi) {
-        typename B::ConstColIterator bik = b[bi.index()].begin();
-        typename B::ConstColIterator biEnd = b[bi.index()].end();
-        for (; bik != biEnd; ++bik) {
-          typename C::ConstColIterator ckj = c[bik.index()].begin();
-          typename C::ConstColIterator ckEnd = c[bik.index()].end();
-          for (; ckj != ckEnd; ++ckj)
-            a[bi.index()][ckj.index()] += scalar * (*bik) * (*ckj);
-        }
-      }
+      rangeForEach(b, [&](auto&& bi, auto&& i) {
+        rangeForEach(bi, [&](auto&& bik, auto&& k) {
+          rangeForEach(c[k], [&](auto&& ckj, auto&& j) {
+            Dune::MatrixVector::addProduct(a[i][j], scalar, bik, ckj);
+          });
+        });
+      });
     }
   };
 
@@ -257,15 +248,11 @@ namespace MatrixVector {
         class ADummy = A,
         std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
     static void addProduct(A& a, const B& b, const C& c) {
-      typename C::ConstRowIterator ci = c.begin();
-      typename C::ConstRowIterator cEnd = c.end();
-      for (; ci != cEnd; ++ci) {
-        typename C::ConstColIterator cik = c[ci.index()].begin();
-        typename C::ConstColIterator ciEnd = c[ci.index()].end();
-        for (; cik != ciEnd; ++cik) {
-          a[ci.index()][cik.index()] += b * (*cik);
-        }
-      }
+      rangeForEach(c, [&](auto&& ci, auto && i) {
+        rangeForEach(ci, [&](auto&& cij, auto && j) {
+          Dune::MatrixVector::addProduct(a[i][j], b, cij);
+        });
+      });
     }
   };
 
@@ -286,15 +273,11 @@ namespace MatrixVector {
         class ADummy = A,
         std::enable_if_t<MatrixTraits<ADummy>::isMatrix, int> SFINAE_Dummy = 0>
     static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
-      typename C::ConstRowIterator ci = c.begin();
-      typename C::ConstRowIterator cEnd = c.end();
-      for (; ci != cEnd; ++ci) {
-        typename C::ConstColIterator cik = c[ci.index()].begin();
-        typename C::ConstColIterator ciEnd = c[ci.index()].end();
-        for (; cik != ciEnd; ++cik) {
-          a[ci.index()][cik.index()] += scalar * b * (*cik);
-        }
-      }
+      rangeForEach(c, [&](auto&& ci, auto&& i) {
+        rangeForEach(ci, [&](auto&& cij, auto&& j) {
+          Dune::MatrixVector::addProduct(a[i][j], scalar, b, cij);
+        });
+      });
     }
   };