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