diff --git a/dune/matrix-vector/axy.hh b/dune/matrix-vector/axy.hh
index c378c89b86c49b5a3f5e3d28cb356b52e673d42c..c626b128dd023d6a2cdfc0f06a98dd8f038cce87 100644
--- a/dune/matrix-vector/axy.hh
+++ b/dune/matrix-vector/axy.hh
@@ -5,6 +5,7 @@
 
 #include "axpy.hh"
 #include "matrixtraits.hh"
+#include "algorithm.hh"
 
 namespace Dune {
 namespace MatrixVector {
@@ -40,22 +41,23 @@ namespace MatrixVector {
     static typename VectorType::field_type Axy(const MatrixType &A,
                                                const VectorType &x,
                                                const VectorType2 &y) {
-      assert(x.N() == A.M());
-      assert(y.N() == A.N());
+      assert(x.size() == A.M());
+      assert(y.size() == A.N());
 
-      typename VectorType::field_type outer = 0.0;
-      typename VectorType2::block_type inner;
-      typename MatrixType::ConstIterator endi = A.end();
-      for (typename MatrixType::ConstIterator i = A.begin(); i != endi; ++i) {
-        inner = 0.0;
-        const size_t iindex = i.index();
-        typename MatrixType::row_type::ConstIterator endj = i->end();
-        for (typename MatrixType::row_type::ConstIterator j = i->begin();
-             j != endj; ++j)
-          addProduct(inner, *j, x[j.index()]);
-        outer += inner * y[iindex];
-      }
-      return outer;
+      using Result = std::decay_t<decltype(y*y)>;
+
+      Result result = 0;
+      sparseRangeFor(A, [&](auto&& Ai, auto i) {
+        // ToDo: Provide a specialization using a single temporary
+        // for the case that each y[i] has the same type.
+        auto Aix = y[i];
+        Aix = 0;
+        sparseRangeFor(Ai, [&](auto&& Aij, auto j) {
+          addProduct(Aix, Aij, x[j]);
+        });
+        result += Aix * y[i];
+      });
+      return result;
     }
 
     template <class VectorType, class VectorType2>
@@ -63,23 +65,23 @@ namespace MatrixVector {
                                                  const VectorType2 &b,
                                                  const VectorType &x,
                                                  const VectorType2 &y) {
-      assert(x.N() == A.M());
-      assert(y.N() == A.N());
-      assert(y.N() == b.N());
+      assert(x.size() == A.M());
+      assert(y.size() == A.N());
+      assert(y.size() == b.size());
+
+      using Result = std::decay_t<decltype(y*y)>;
 
-      typename VectorType::field_type outer = 0.0;
-      typename VectorType2::block_type inner;
-      typename MatrixType::ConstIterator endi = A.end();
-      for (typename MatrixType::ConstIterator i = A.begin(); i != endi; ++i) {
-        const size_t iindex = i.index();
-        inner = b[iindex];
-        typename MatrixType::row_type::ConstIterator endj = i->end();
-        for (typename MatrixType::row_type::ConstIterator j = i->begin();
-             j != endj; ++j)
-          subtractProduct(inner, *j, x[j.index()]);
-        outer += inner * y[iindex];
-      }
-      return outer;
+      Result result = 0;
+      sparseRangeFor(A, [&](auto&& Ai, auto i) {
+        // ToDo: Provide a specialization using a single temporary
+        // for the case that each b[i] has the same type.
+        auto Aix = b[i];
+        sparseRangeFor(Ai, [&](auto&& Aij, auto j) {
+          subtractProduct(Aix, Aij, x[j]);
+        });
+        result += Aix * y[i];
+      });
+      return result;
     }
   };