diff --git a/dune/matrix-vector/algorithm.hh b/dune/matrix-vector/algorithm.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b9817e04eff80307b0a0db69f2e77ef6c0f1f0a4
--- /dev/null
+++ b/dune/matrix-vector/algorithm.hh
@@ -0,0 +1,53 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_MATRIX_VECTOR_ALGORITHM_HH
+#define DUNE_MATRIX_VECTOR_ALGORITHM_HH
+
+#include <type_traits>
+
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+
+namespace Dune {
+namespace MatrixVector {
+
+template <class C>
+struct IsMultiTypeBlockContainer : public std::false_type {};
+
+template <class... Args>
+struct IsMultiTypeBlockContainer<MultiTypeBlockMatrix<Args...>>
+  : public std::true_type {};
+
+template <class... Args>
+struct IsMultiTypeBlockContainer<MultiTypeBlockVector<Args...>>
+  : public std::true_type {};
+
+/**
+ * \brief Hybrid for loop over sparse range
+ */
+template <class Range, class F,
+          typename = std::enable_if_t<IsMultiTypeBlockContainer<Range>::value>>
+void rangeForEach(const Range& range, F&& f) {
+  using namespace Dune::Hybrid;
+  forEach(integralRange(size(range)), [&](auto&& i) { f(range[i], i); });
+}
+
+/**
+ * \brief Hybrid for loop over sparse range
+ */
+template<class Range, class F>
+void rangeForEach(Range&& range, F&& f)
+{
+  for (auto it = range.begin(); it != range.end(); ++it)
+    f(*it, it.index());
+}
+
+
+
+} // namespace MatrixVector
+} // namespace Dune
+
+
+#endif // DUNE_MATRIX_VECTOR_ALGORITHM_HH
diff --git a/dune/matrix-vector/axpy.hh b/dune/matrix-vector/axpy.hh
index 792dc034b71d522634ca3aab203f1b54abb6e29e..72442be7a86437c61c924d0cee14f85e2dd1dad7 100644
--- a/dune/matrix-vector/axpy.hh
+++ b/dune/matrix-vector/axpy.hh
@@ -8,11 +8,96 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/scaledidmatrix.hh>
 
+#include "algorithm.hh"
 #include "matrixtraits.hh"
 #include "scalartraits.hh"
 
 namespace Dune {
 namespace MatrixVector {
+
+  /// forward declarations of internal helper classes for product operations
+  template <class A, class B, class C, bool AisScalar, bool BisScalar,
+            bool CisScalar>
+  struct ProductHelper;
+
+  template <class A, class Scalar, class B, class C, bool AisScalar,
+            bool BisScalar, bool CisScalar>
+  struct ScaledProductHelper;
+
+
+  /** \brief Add a product to some matrix or vector
+   *
+   * This function computes a+=b*c.
+   *
+   * This function should tolerate all meaningful
+   * combinations of scalars, vectors, and matrices.
+   *
+   * a,b,c could be matrices with appropriate
+   * dimensions. But b can also always be a scalar
+   * represented by a 1-dim vector or a 1 by 1 matrix.
+   */
+  template <class A, class B, class C>
+  void addProduct(A& a, const B& b, const C& c) {
+    ProductHelper<A, B, C, ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar,
+                  ScalarTraits<C>::isScalar>::addProduct(a, b, c);
+  }
+
+  /** \brief Subtract a product from some matrix or vector
+   *
+   * This function computes a-=b*c.
+   *
+   * This function should tolerate all meaningful
+   * combinations of scalars, vectors, and matrices.
+   *
+   * a,b,c could be matrices with appropriate
+   * dimensions. But b can also always be a scalar
+   * represented by a 1-dim vector or a 1 by 1 matrix.
+   */
+  template <class A, class B, class C>
+  void subtractProduct(A& a, const B& b, const C& c) {
+    ScaledProductHelper<A, int, B, C, ScalarTraits<A>::isScalar,
+                        ScalarTraits<B>::isScalar,
+                        ScalarTraits<C>::isScalar>::addProduct(a, -1, b, c);
+  }
+
+  /** \brief Add a scaled product to some matrix or vector
+   *
+   * This function computes a+=b*c*d.
+   *
+   * This function should tolerate all meaningful
+   * combinations of scalars, vectors, and matrices.
+   *
+   * a,c,d could be matrices with appropriate dimensions. But b must
+   * (currently) and c can also always be a scalar represented by a
+   * 1-dim vector or a 1 by 1 matrix.
+   */
+  template <class A, class B, class C, class D>
+  typename std::enable_if_t<ScalarTraits<B>::isScalar, void> addProduct(
+      A& a, const B& b, const C& c, const D& d) {
+    ScaledProductHelper<A, B, C, D, ScalarTraits<A>::isScalar,
+                        ScalarTraits<C>::isScalar,
+                        ScalarTraits<D>::isScalar>::addProduct(a, b, c, d);
+  }
+
+  /** \brief Subtract a scaled product from some matrix or vector
+   *
+   * This function computes a-=b*c*d.
+   *
+   * This function should tolerate all meaningful
+   * combinations of scalars, vectors, and matrices.
+   *
+   * a,c,d could be matrices with appropriate dimensions. But b must
+   * (currently) and c can also always be a scalar represented by a
+   * 1-dim vector or a 1 by 1 matrix.
+   */
+  template <class A, class B, class C, class D>
+  typename std::enable_if_t<ScalarTraits<B>::isScalar, void>
+  subtractProduct(A& a, const B& b, const C& c, const D& d) {
+    ScaledProductHelper<A, B, C, D, ScalarTraits<A>::isScalar,
+                        ScalarTraits<C>::isScalar,
+                        ScalarTraits<D>::isScalar>::addProduct(a, -b, c, d);
+  }
+
   /** \brief Internal helper class for product operations
    *
    */
@@ -30,18 +115,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 +141,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);
+          });
+        });
+      });
     }
   };
 
@@ -83,12 +158,10 @@ namespace MatrixVector {
     typedef Dune::FieldMatrix<T, n, m> B;
     typedef Dune::FieldMatrix<T, m, k> C;
     static void addProduct(A& a, const B& b, const C& c) {
-      for (size_t row = 0; row < n; ++row) {
-        for (size_t col = 0; col < k; ++col) {
+      for (size_t row = 0; row < n; ++row)
+        for (size_t col = 0; col < k; ++col)
           for (size_t i = 0; i < m; ++i)
             a[row][col] += b[row][i] * c[i][col];
-        }
-      }
     }
   };
 
@@ -100,12 +173,10 @@ namespace MatrixVector {
     typedef Dune::FieldMatrix<T, n, m> C;
     typedef Dune::FieldMatrix<T, m, k> D;
     static void addProduct(A& a, const T& b, const C& c, const D& d) {
-      for (size_t row = 0; row < n; ++row) {
-        for (size_t col = 0; col < k; ++col) {
+      for (size_t row = 0; row < n; ++row)
+        for (size_t col = 0; col < k; ++col)
           for (size_t i = 0; i < m; ++i)
             a[row][col] += b * c[row][i] * d[i][col];
-        }
-      }
     }
   };
 
@@ -257,15 +328,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,65 +353,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);
-        }
-      }
-    }
-  };
-
-  template <class ABlock, class ScalarB, class CBlock>
-  struct ProductHelper<Dune::BCRSMatrix<ABlock>, ScalarB,
-                       Dune::BCRSMatrix<CBlock>, false, true, false> {
-    typedef Dune::BCRSMatrix<ABlock> A;
-    typedef ScalarB B;
-    typedef Dune::BCRSMatrix<CBlock> C;
-
-    static void addProduct(A& a, const B& b, const C& c) {
-      // Workaround for the fact that Dune::BCRSMatrix assumes
-      // its blocks to have an axpy() implementation, yet
-      // Dune::DiagonalMatrix does not.
-
-      // a.axpy(b, c);
-
-      for (typename C::ConstIterator i = c.begin(); i != c.end(); ++i) {
-        const size_t iindex = i.index();
-        for (typename C::row_type::ConstIterator j = i->begin(); j != i->end();
-             ++j)
-          ProductHelper<typename A::block_type, B, typename C::block_type,
-                        false, true, false>::addProduct(a[iindex][j.index()], b,
-                                                        *j);
-      }
-    }
-  };
-
-  template <class ABlock, class Scalar, class ScalarB, class CBlock>
-  struct ScaledProductHelper<Dune::BCRSMatrix<ABlock>, Scalar, ScalarB,
-                             Dune::BCRSMatrix<CBlock>, false, true, false> {
-    typedef Dune::BCRSMatrix<ABlock> A;
-    typedef ScalarB B;
-    typedef Dune::BCRSMatrix<CBlock> C;
-
-    static void addProduct(A& a, const Scalar& scalar, const B& b, const C& c) {
-      // Workaround for the fact that Dune::BCRSMatrix assumes
-      // its blocks to have an axpy() implementation, yet
-      // Dune::DiagonalMatrix does not.
-
-      // a.axpy(scalar*b, c);
-
-      for (typename C::ConstIterator i = c.begin(); i != c.end(); ++i) {
-        const size_t iindex = i.index();
-        for (typename C::row_type::ConstIterator j = i->begin(); j != i->end();
-             ++j)
-          ProductHelper<typename A::block_type, B, typename C::block_type,
-                        false, true, false>::addProduct(a[iindex][j.index()],
-                                                        scalar * b, *j);
-      }
+      rangeForEach(c, [&](auto&& ci, auto&& i) {
+        rangeForEach(ci, [&](auto&& cij, auto&& j) {
+          Dune::MatrixVector::addProduct(a[i][j], scalar, b, cij);
+        });
+      });
     }
   };
 
@@ -385,78 +398,6 @@ namespace MatrixVector {
     }
   };
 
-  /** \brief Add a product to some matrix or vector
-   *
-   * This function computes a+=b*c.
-   *
-   * This function should tolerate all meaningful
-   * combinations of scalars, vectors, and matrices.
-   *
-   * a,b,c could be matrices with appropriate
-   * dimensions. But b can also always be a scalar
-   * represented by a 1-dim vector or a 1 by 1 matrix.
-   */
-  template <class A, class B, class C>
-  void addProduct(A& a, const B& b, const C& c) {
-    ProductHelper<A, B, C, ScalarTraits<A>::isScalar, ScalarTraits<B>::isScalar,
-                  ScalarTraits<C>::isScalar>::addProduct(a, b, c);
-  }
-
-  /** \brief Subtract a product from some matrix or vector
-   *
-   * This function computes a-=b*c.
-   *
-   * This function should tolerate all meaningful
-   * combinations of scalars, vectors, and matrices.
-   *
-   * a,b,c could be matrices with appropriate
-   * dimensions. But b can also always be a scalar
-   * represented by a 1-dim vector or a 1 by 1 matrix.
-   */
-  template <class A, class B, class C>
-  void subtractProduct(A& a, const B& b, const C& c) {
-    ScaledProductHelper<A, int, B, C, ScalarTraits<A>::isScalar,
-                        ScalarTraits<B>::isScalar,
-                        ScalarTraits<C>::isScalar>::addProduct(a, -1, b, c);
-  }
-
-  /** \brief Add a scaled product to some matrix or vector
-   *
-   * This function computes a+=b*c*d.
-   *
-   * This function should tolerate all meaningful
-   * combinations of scalars, vectors, and matrices.
-   *
-   * a,c,d could be matrices with appropriate dimensions. But b must
-   * (currently) and c can also always be a scalar represented by a
-   * 1-dim vector or a 1 by 1 matrix.
-   */
-  template <class A, class B, class C, class D>
-  typename std::enable_if<ScalarTraits<B>::isScalar, void>::type addProduct(
-      A& a, const B& b, const C& c, const D& d) {
-    ScaledProductHelper<A, B, C, D, ScalarTraits<A>::isScalar,
-                        ScalarTraits<C>::isScalar,
-                        ScalarTraits<D>::isScalar>::addProduct(a, b, c, d);
-  }
-
-  /** \brief Subtract a scaled product from some matrix or vector
-   *
-   * This function computes a-=b*c*d.
-   *
-   * This function should tolerate all meaningful
-   * combinations of scalars, vectors, and matrices.
-   *
-   * a,c,d could be matrices with appropriate dimensions. But b must
-   * (currently) and c can also always be a scalar represented by a
-   * 1-dim vector or a 1 by 1 matrix.
-   */
-  template <class A, class B, class C, class D>
-  typename std::enable_if<ScalarTraits<B>::isScalar, void>::type
-  subtractProduct(A& a, const B& b, const C& c, const D& d) {
-    ScaledProductHelper<A, B, C, D, ScalarTraits<A>::isScalar,
-                        ScalarTraits<C>::isScalar,
-                        ScalarTraits<D>::isScalar>::addProduct(a, -b, c, d);
-  }
 }
 }
 #endif
diff --git a/dune/matrix-vector/matrixtraits.hh b/dune/matrix-vector/matrixtraits.hh
index da72f651f436099afa3f39ed82e41cfa260a7cdd..5b673ed6d66a37d523f6a8135b54a31fd2293e0f 100644
--- a/dune/matrix-vector/matrixtraits.hh
+++ b/dune/matrix-vector/matrixtraits.hh
@@ -4,6 +4,7 @@
 #include <dune/common/diagonalmatrix.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
 #include <dune/istl/scaledidmatrix.hh>
 
 namespace Dune {
@@ -50,6 +51,12 @@ namespace Dune {
     {
         enum { isMatrix=true };
     };
+
+    template<class... T>
+    struct MatrixTraits<MultiTypeBlockMatrix<T...> >
+    {
+        enum { isMatrix=true };
+    };
   }
 }
 #endif
diff --git a/dune/matrix-vector/scalartraits.hh b/dune/matrix-vector/scalartraits.hh
index 5fabfc15b132f6d55aa292e68e4997b8765ebf3a..0894730ad29590547ce0cb9883ad2338125c1c5a 100644
--- a/dune/matrix-vector/scalartraits.hh
+++ b/dune/matrix-vector/scalartraits.hh
@@ -3,6 +3,7 @@
 
 #include <dune/common/diagonalmatrix.hh>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/typetraits.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/scaledidmatrix.hh>
 
@@ -16,7 +17,7 @@ namespace MatrixVector {
   template <class T>
   struct ScalarTraits {
     enum {
-      isScalar = (std::is_scalar<T>::value and not std::is_pointer<T>::value)
+      isScalar = Dune::IsNumber<T>::value
     };
   };
 
diff --git a/dune/matrix-vector/test/arithmetictest.cc b/dune/matrix-vector/test/arithmetictest.cc
index e319d122eb0fdf6a9fa970ff602410be446cd8a6..be4ef793f1c977ab74ea0145e17c0292c835c790 100644
--- a/dune/matrix-vector/test/arithmetictest.cc
+++ b/dune/matrix-vector/test/arithmetictest.cc
@@ -9,7 +9,10 @@
 #include <dune/common/diagonalmatrix.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
+#include <dune/common/indices.hh>
 #include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
 #include <dune/istl/scaledidmatrix.hh>
 
 #include "common.hh"
@@ -636,6 +639,37 @@ public:
       passed =
           passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
     }
+    // case MM += scalar*MM
+    {
+        using ResultType = Dune::MultiTypeBlockMatrix<
+            Dune::MultiTypeBlockVector<
+                Dune::FieldMatrix<FT, 1, 1>, Dune::FieldMatrix<FT, 1, 2> >,
+        Dune::MultiTypeBlockVector<
+                Dune::FieldMatrix<FT, 2, 1>, Dune::FieldMatrix<FT, 2, 2> > >;
+        ResultType multiTypeBlockMatrix_a, multiTypeBlockMatrix_b;
+        using namespace Dune::Indices;
+        multiTypeBlockMatrix_a[_0][_0] = {{1}};
+        multiTypeBlockMatrix_a[_0][_1] = {{2, 3}};
+        multiTypeBlockMatrix_a[_1][_0] = {{4}, {5}};
+        multiTypeBlockMatrix_a[_1][_1] = {{6, 7}, {8, 9}};
+        multiTypeBlockMatrix_b[_0][_0] = {{9}};
+        multiTypeBlockMatrix_b[_0][_1] = {{8, 7}};
+        multiTypeBlockMatrix_b[_1][_0] = {{6}, {5}};
+        multiTypeBlockMatrix_b[_1][_1] = {{4, 3}, {2, 1}};
+        addProduct(multiTypeBlockMatrix_a, scalar_b, multiTypeBlockMatrix_b);
+        ResultType multiTypeBlockMatrix_check;
+        multiTypeBlockMatrix_check[_0][_0] = {{19}};
+        multiTypeBlockMatrix_check[_0][_1] = {{18, 17}};
+        multiTypeBlockMatrix_check[_1][_0] = {{16}, {15}};
+        multiTypeBlockMatrix_check[_1][_1] = {{14, 13}, {12, 11}};
+        using namespace Dune::Hybrid;
+        forEach(integralRange(size(multiTypeBlockMatrix_a)), [&](auto&& i) {
+            forEach(integralRange(size(multiTypeBlockMatrix_a[i])), [&](auto && j) {
+                passed = passed and isCloseDune(multiTypeBlockMatrix_a[i][j],
+                                                multiTypeBlockMatrix_check[i][j]);
+            });
+        });
+    }
 
     // RECTANGULAR MATRICES