Implement axy using sparseRangeFor()

This way it works for multi-type matrices again.
This helps to fix the dune-tnnmg test which currently
breaks.

Notice that this may be slightly less efficient, because
we use a new temporary for each loop instance. This is
necessary because we may need a different type each time.
This may be fixed by adding a specialization for the
case that the corresponding vector is a classic dynamic
container.

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;
     }
   };