diff --git a/dune/matrix-vector/CMakeLists.txt b/dune/matrix-vector/CMakeLists.txt
index 4d090e0b4e56ebdb8612cbcc657f23116db8c1e6..07a1c2cbffa0a361f3f64fcac212a81a3a0d2b2e 100644
--- a/dune/matrix-vector/CMakeLists.txt
+++ b/dune/matrix-vector/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory(test)
+
 #install headers
 install(FILES
   addtodiagonal.hh
diff --git a/dune/matrix-vector/test/CMakeLists.txt b/dune/matrix-vector/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ab5d8968bb8072bf3a64118f9ed22ac91fd8431b
--- /dev/null
+++ b/dune/matrix-vector/test/CMakeLists.txt
@@ -0,0 +1,2 @@
+dune_add_test(SOURCES arithmetictest.cc)
+dune_add_test(SOURCES staticmatrixtoolstest.cc)
diff --git a/dune/matrix-vector/test/arithmetictest.cc b/dune/matrix-vector/test/arithmetictest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e319d122eb0fdf6a9fa970ff602410be446cd8a6
--- /dev/null
+++ b/dune/matrix-vector/test/arithmetictest.cc
@@ -0,0 +1,1834 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <limits>
+
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+#include "common.hh"
+
+#include "../axpy.hh"
+#include "../axy.hh"
+#include "../crossproduct.hh"
+#include "../transpose.hh"
+
+using namespace Dune::MatrixVector;
+
+template <class FT>
+struct ArithmeticTestSuite {
+  bool verbose;
+
+  const FT tol = 10 * std::numeric_limits<FT>::epsilon();
+
+  ArithmeticTestSuite(bool verbose = false) : verbose(verbose) {}
+
+  bool check() {
+    bool passed = true;
+
+    passed = passed and checkProductHelperMethods();
+    passed = passed and checkScaledProductHelperMethods();
+
+    passed = passed and checkCrossProduct();
+
+    passed = passed and checkAxy();
+    passed = passed and checkBmAxy();
+
+    passed = passed and checkTranspose();
+
+    passed = passed and checkBCRSSubPattern();
+
+    return passed;
+  }
+
+private:
+  bool checkProductHelperMethods() {
+    if (verbose)
+      std::cout << "checkProductHelperMethods<" << Dune::className<FT>() << ">"
+                << std::endl;
+
+    bool passed = true;
+
+    passed = passed and checkProductHelperMethods_MatrixMatrixMatrix();
+    passed = passed and checkProductHelperMethods_MatrixScalarMatrix();
+
+    passed = passed and checkProductHelperMethods_VectorMatrixVector();
+    passed = passed and checkProductHelperMethods_VectorScalarVector();
+
+    passed = passed and checkProductHelperMethods_ScalarScalarScalar();
+
+    return passed;
+  }
+
+  bool checkScaledProductHelperMethods() {
+    if (verbose)
+      std::cout << "checkScaledProductHelperMethods<" << Dune::className<FT>()
+                << ">" << std::endl;
+
+    bool passed = true;
+
+    passed = passed and checkScaledProductHelperMethods_MatrixMatrixMatrix();
+    passed = passed and checkScaledProductHelperMethods_MatrixScalarMatrix();
+
+    passed = passed and checkScaledProductHelperMethods_VectorMatrixVector();
+    passed = passed and checkScaledProductHelperMethods_VectorScalarVector();
+
+    passed = passed and checkScaledProductHelperMethods_ScalarScalarScalar();
+
+    return passed;
+  }
+
+  bool checkCrossProduct() {
+    if (verbose)
+      std::cout << "checkCrossProduct<" << Dune::className<FT>() << ">"
+                << std::endl;
+
+    bool passed = true;
+
+    typedef Dune::FieldVector<FT, 3> VectorType;
+    const VectorType fieldVector_x = {1, 2, 3}, fieldVector_y = {2, 3, 1};
+
+    passed = passed and isCloseDune(crossProduct(fieldVector_x, fieldVector_y),
+                                    VectorType{-7, 5, -1});
+
+    return passed;
+  }
+
+  bool checkAxy() {
+    if (verbose)
+      std::cout << "checkAxy<" << Dune::className<FT>() << ">" << std::endl;
+
+    bool passed = true;
+
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_A = {{1, 2}, {3, 4}};
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_A = {3, 2};
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_A(2);
+
+    const Dune::FieldVector<FT, 2> fieldVector_x = {1, 2},
+                                   fieldVector_y = {2, 3};
+
+    passed =
+        passed and
+        isCloseAbs(Axy(squareFieldMatrix_A, fieldVector_x, fieldVector_y), 43);
+    passed =
+        passed and
+        isCloseAbs(Axy(diagonalMatrix_A, fieldVector_x, fieldVector_y), 18);
+    passed = passed and
+             isCloseAbs(
+                 Axy(scaledIdentityMatrix_A, fieldVector_x, fieldVector_y), 16);
+
+    return passed;
+  }
+
+public:
+  bool checkBmAxy() {
+    if (verbose)
+      std::cout << "checkBmAxy<" << Dune::className<FT>() << ">" << std::endl;
+
+    bool passed = true;
+
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_A = {{1, 2}, {3, 4}};
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_A = {3, 2};
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_A(2);
+
+    const Dune::FieldVector<FT, 2> fieldVector_x = {1, 2},
+                                   fieldVector_y = {2, 3},
+                                   fieldVector_b = {10, 13};
+
+    passed = passed and isCloseAbs(bmAxy(squareFieldMatrix_A, fieldVector_b,
+                                         fieldVector_x, fieldVector_y),
+                                   16);
+    passed = passed and isCloseAbs(bmAxy(diagonalMatrix_A, fieldVector_b,
+                                         fieldVector_x, fieldVector_y),
+                                   41);
+    passed = passed and isCloseAbs(bmAxy(scaledIdentityMatrix_A, fieldVector_b,
+                                         fieldVector_x, fieldVector_y),
+                                   43);
+
+    return passed;
+  }
+
+  bool checkTranspose() {
+    if (verbose)
+      std::cout << "checkTranspose<" << Dune::className<FT>() << ">"
+                << std::endl;
+
+    bool passed = true;
+
+    typedef Dune::FieldMatrix<FT, 3, 3> SquareFieldMatrix;
+    typedef Dune::FieldMatrix<FT, 2, 4> RectFieldMatrix;
+    typedef Dune::DiagonalMatrix<FT, 2> DiagonalMatrix;
+    typedef Dune::ScaledIdentityMatrix<FT, 2> ScaledIdentityMatrix;
+
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<FT, 2, 4>> BCRSofRectFieldMatrix;
+
+    const SquareFieldMatrix squareFieldMatrix_A = {
+        {1, 2, 3}, {3, 4, 5}, {5, 6, 7}};
+    const Transposed<SquareFieldMatrix> squareFieldMatrix_check = {
+        {1, 3, 5}, {2, 4, 6}, {3, 5, 7}};
+    Transposed<SquareFieldMatrix> squareFieldMatrix_At(0);
+
+    const RectFieldMatrix rectFieldMatrix_A = {{1, 2, 3, 4}, {3, 4, 5, 6}};
+    const Transposed<RectFieldMatrix> rectFieldMatrix_check = {
+        {1, 3}, {2, 4}, {3, 5}, {4, 6}};
+    Transposed<RectFieldMatrix> rectFieldMatrix_At(0);
+
+    const DiagonalMatrix diagonalMatrix_A = {3, 2};
+    const Transposed<DiagonalMatrix> diagonalMatrix_check = {3, 2};
+    Transposed<DiagonalMatrix> diagonalMatrix_At(0);
+
+    const ScaledIdentityMatrix scaledIdentityMatrix_A(2);
+    const Transposed<ScaledIdentityMatrix> scaledIdentityMatrix_check(2);
+    Transposed<ScaledIdentityMatrix> scaledIdentityMatrix_At(0);
+
+    BCRSofRectFieldMatrix bcrsRectFieldMatrix_A(2, 3,
+                                                BCRSofRectFieldMatrix::random);
+    bcrsRectFieldMatrix_A.setrowsize(0, 2);
+    bcrsRectFieldMatrix_A.setrowsize(1, 1);
+    bcrsRectFieldMatrix_A.endrowsizes();
+
+    bcrsRectFieldMatrix_A.addindex(0, 0);
+    bcrsRectFieldMatrix_A.addindex(0, 2);
+    bcrsRectFieldMatrix_A.addindex(1, 0);
+    bcrsRectFieldMatrix_A.endindices();
+
+    bcrsRectFieldMatrix_A[0][0] = RectFieldMatrix{{1, 2, 3, 4}, {3, 4, 5, 6}};
+    bcrsRectFieldMatrix_A[0][2] = RectFieldMatrix{{2, 2, 3, 3}, {1, 1, 4, 2}};
+    bcrsRectFieldMatrix_A[1][0] = RectFieldMatrix{{3, 1, 2, 5}, {4, 3, 2, 3}};
+
+    Transposed<BCRSofRectFieldMatrix> bcrsRectFieldMatrix_check(
+        3, 2, Transposed<BCRSofRectFieldMatrix>::random);
+    bcrsRectFieldMatrix_check.setrowsize(0, 2);
+    bcrsRectFieldMatrix_check.setrowsize(1, 0);
+    bcrsRectFieldMatrix_check.setrowsize(2, 1);
+    bcrsRectFieldMatrix_check.endrowsizes();
+
+    bcrsRectFieldMatrix_check.addindex(0, 0);
+    bcrsRectFieldMatrix_check.addindex(0, 1);
+    bcrsRectFieldMatrix_check.addindex(2, 0);
+    bcrsRectFieldMatrix_check.endindices();
+
+    bcrsRectFieldMatrix_check[0][0] =
+        Transposed<RectFieldMatrix>{{1, 3}, {2, 4}, {3, 5}, {4, 6}};
+    bcrsRectFieldMatrix_check[2][0] =
+        Transposed<RectFieldMatrix>{{2, 1}, {2, 1}, {3, 4}, {3, 2}};
+    bcrsRectFieldMatrix_check[0][1] =
+        Transposed<RectFieldMatrix>{{3, 4}, {1, 3}, {2, 2}, {5, 3}};
+
+    Transposed<BCRSofRectFieldMatrix> bcrsRectFieldMatrix_At;
+
+    transpose(squareFieldMatrix_A, squareFieldMatrix_At);
+    passed =
+        passed and isCloseDune(squareFieldMatrix_At, squareFieldMatrix_check);
+
+    transpose(rectFieldMatrix_A, rectFieldMatrix_At);
+    passed = passed and isCloseDune(rectFieldMatrix_At, rectFieldMatrix_check);
+
+    transpose(diagonalMatrix_A, diagonalMatrix_At);
+    passed = passed and isCloseDune(diagonalMatrix_At, diagonalMatrix_check);
+
+    transpose(scaledIdentityMatrix_A, scaledIdentityMatrix_At);
+    passed = passed and
+             isCloseDune(scaledIdentityMatrix_At, scaledIdentityMatrix_check);
+
+    transpose(bcrsRectFieldMatrix_A, bcrsRectFieldMatrix_At);
+    passed = passed and
+             isCloseDune(bcrsRectFieldMatrix_At, bcrsRectFieldMatrix_check);
+
+    return passed;
+  }
+
+  bool checkProductHelperMethods_MatrixMatrixMatrix() {
+    bool passed = true;
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_b = {{1, 2}, {3, 4}},
+
+                                      squareFieldMatrix_c = {{2, 4}, {3, 5}};
+
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_b = {3, 2},
+                                      diagonalMatrix_c = {2, 1};
+
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_b(2),
+        scaledIdentityMatrix_c(4);
+    // case FM += FM*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 14}, {18, 32}};
+
+      addProduct(squareFieldMatrix_a, squareFieldMatrix_b, squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, squareFieldMatrix_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += DM*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{6, 12}, {6, 10}};
+
+      addProduct(squareFieldMatrix_a, diagonalMatrix_b, squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, diagonalMatrix_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += FM*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{2, 2}, {6, 4}};
+
+      addProduct(squareFieldMatrix_a, squareFieldMatrix_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, squareFieldMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += DM*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{6, 0}, {0, 2}};
+
+      addProduct(squareFieldMatrix_a, diagonalMatrix_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, diagonalMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += SM*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 8}, {6, 10}};
+
+      addProduct(squareFieldMatrix_a, scaledIdentityMatrix_b,
+                 squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaledIdentityMatrix_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += FM*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 8}, {12, 16}};
+
+      addProduct(squareFieldMatrix_a, squareFieldMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, squareFieldMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += SM*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 0}, {0, 8}};
+
+      addProduct(squareFieldMatrix_a, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaledIdentityMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += DM*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{12, 0}, {0, 8}};
+
+      addProduct(squareFieldMatrix_a, diagonalMatrix_b, scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, diagonalMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += SM*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 0}, {0, 2}};
+
+      addProduct(squareFieldMatrix_a, scaledIdentityMatrix_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaledIdentityMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case DM += DM*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {6, 2};
+
+      addProduct(diagonalMatrix_a, diagonalMatrix_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, diagonalMatrix_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += DM*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {12, 8};
+
+      addProduct(diagonalMatrix_a, diagonalMatrix_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, diagonalMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += SM*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {4, 2};
+
+      addProduct(diagonalMatrix_a, scaledIdentityMatrix_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaledIdentityMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += SM*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {8, 8};
+
+      addProduct(diagonalMatrix_a, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaledIdentityMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+
+    // case SM += SM*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(8);
+
+      addProduct(scaledIdentityMatrix_a, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, scaledIdentityMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // Extra test: SM += a*SM*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(8);
+
+      addProduct(scaledIdentityMatrix_a, 2, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      addProduct(scaledIdentityMatrix_a, -1, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, 2, scaledIdentityMatrix_b,
+                      scaledIdentityMatrix_c);
+      addProduct(scaledIdentityMatrix_check, 1, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // RECTANGULAR MATRICES
+
+    const Dune::FieldMatrix<FT, 3, 2> rectFieldMatrix_b = {
+        {1, 2}, {3, 4}, {5, 6}};
+
+    const Dune::FieldMatrix<FT, 2, 4> rectFieldMatrix_c = {{2, 3, 2, 4},
+                                                           {1, 4, 3, 5}};
+    // case FM<3,4> += FM<3,2>*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 3, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {
+              {4, 11, 8, 14}, {10, 25, 18, 32}, {16, 39, 28, 50}};
+
+      addProduct(rectFieldMatrix_a, rectFieldMatrix_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, rectFieldMatrix_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<2,4> += DM<2>*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{6, 9, 6, 12}, {2, 8, 6, 10}};
+
+      addProduct(rectFieldMatrix_a, diagonalMatrix_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, diagonalMatrix_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<3,2> += FM<3,2>*DM<2>
+    {
+      typedef Dune::FieldMatrix<FT, 3, 2> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{2, 2}, {6, 4}, {10, 6}};
+
+      addProduct(rectFieldMatrix_a, rectFieldMatrix_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, rectFieldMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<2,4> += SM<2>*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{4, 6, 4, 8}, {2, 8, 6, 10}};
+
+      addProduct(rectFieldMatrix_a, scaledIdentityMatrix_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaledIdentityMatrix_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<3,2> += FM<3,2>*SM<2>
+    {
+      typedef Dune::FieldMatrix<FT, 3, 2> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{4, 8}, {12, 16}, {20, 24}};
+
+      addProduct(rectFieldMatrix_a, rectFieldMatrix_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, rectFieldMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    return passed;
+  }
+
+  bool checkProductHelperMethods_MatrixScalarMatrix() {
+    bool passed = true;
+
+    const FT scalar_b = 2;
+    const double double_b = 2.0;
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_c = {{2, 4}, {3, 5}};
+
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_c = {2, 1};
+
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_c(4);
+
+    // case FM += scalar*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 8}, {6, 10}};
+
+      addProduct(squareFieldMatrix_a, scalar_b, squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scalar_b, squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += scalar*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 0}, {0, 2}};
+
+      addProduct(squareFieldMatrix_a, scalar_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scalar_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += scalar*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 0}, {0, 8}};
+
+      addProduct(squareFieldMatrix_a, scalar_b, scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scalar_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case DM += scalar*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {4, 2};
+
+      addProduct(diagonalMatrix_a, scalar_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scalar_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += scalar*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {8, 8};
+
+      addProduct(diagonalMatrix_a, scalar_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scalar_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case SM += scalar*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(8);
+
+      addProduct(scaledIdentityMatrix_a, scalar_b, scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, scalar_b,
+                      scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // RECTANGULAR MATRICES
+
+    const Dune::FieldMatrix<FT, 2, 4> rectFieldMatrix_c = {{2, 3, 2, 4},
+                                                           {1, 4, 3, 5}};
+    // case FM<2,4> += scalar*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{4, 6, 4, 8}, {2, 8, 6, 10}};
+
+      addProduct(rectFieldMatrix_a, scalar_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scalar_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+
+    // and again for fixed scalar type double
+    // SQUARE MATRICES
+    // case FM += double*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 8}, {6, 10}};
+
+      addProduct(squareFieldMatrix_a, double_b, squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, double_b, squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += double*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 0}, {0, 2}};
+
+      addProduct(squareFieldMatrix_a, double_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, double_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += double*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 0}, {0, 8}};
+
+      addProduct(squareFieldMatrix_a, double_b, scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, double_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case DM += double*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {4, 2};
+
+      addProduct(diagonalMatrix_a, double_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, double_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += double*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {8, 8};
+
+      addProduct(diagonalMatrix_a, double_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, double_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case SM += double*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(8);
+
+      addProduct(scaledIdentityMatrix_a, double_b, scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, double_b,
+                      scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // RECTANGULAR MATRICES
+    // case FM<2,4> += double*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{4, 6, 4, 8}, {2, 8, 6, 10}};
+
+      addProduct(rectFieldMatrix_a, double_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, double_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+
+    // BLOCKED MATRICES of RECTANGULAR MATRICES
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> RectFieldMatrix;
+      typedef Dune::BCRSMatrix<RectFieldMatrix> BCRSofRectFieldMatrix;
+
+      BCRSofRectFieldMatrix bcrsRectFieldMatrix_c(
+          2, 3, BCRSofRectFieldMatrix::random);
+      bcrsRectFieldMatrix_c.setrowsize(0, 2);
+      bcrsRectFieldMatrix_c.setrowsize(1, 1);
+      bcrsRectFieldMatrix_c.endrowsizes();
+
+      bcrsRectFieldMatrix_c.addindex(0, 0);
+      bcrsRectFieldMatrix_c.addindex(0, 2);
+      bcrsRectFieldMatrix_c.addindex(1, 0);
+      bcrsRectFieldMatrix_c.endindices();
+
+      bcrsRectFieldMatrix_c[0][0] = RectFieldMatrix{{1, 2, 3, 4}, {3, 4, 5, 6}};
+      bcrsRectFieldMatrix_c[0][2] = RectFieldMatrix{{2, 2, 3, 3}, {1, 1, 4, 2}};
+      bcrsRectFieldMatrix_c[1][0] = RectFieldMatrix{{3, 1, 2, 5}, {4, 3, 2, 3}};
+
+      BCRSofRectFieldMatrix bcrsRectFieldMatrix_check(bcrsRectFieldMatrix_c);
+      bcrsRectFieldMatrix_check *= scalar_b;
+
+      BCRSofRectFieldMatrix bcrsRectFieldMatrix_a(bcrsRectFieldMatrix_c),
+          bcrsRectFieldMatrix_zero(bcrsRectFieldMatrix_c);
+      bcrsRectFieldMatrix_a = bcrsRectFieldMatrix_zero = 0;
+
+      addProduct(bcrsRectFieldMatrix_a, scalar_b, bcrsRectFieldMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsRectFieldMatrix_check, bcrsRectFieldMatrix_a);
+
+      subtractProduct(bcrsRectFieldMatrix_check, scalar_b,
+                      bcrsRectFieldMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsRectFieldMatrix_check, bcrsRectFieldMatrix_zero);
+    }
+
+    // BLOCKED MATRICES of DIAGONAL MATRICES
+    {
+      typedef Dune::DiagonalMatrix<FT, 3> DiagMatrix;
+      typedef Dune::BCRSMatrix<DiagMatrix> BCRSofDiagonalMatrix;
+
+      BCRSofDiagonalMatrix bcrsDiagonalMatrix_c(2, 3,
+                                                BCRSofDiagonalMatrix::random);
+      bcrsDiagonalMatrix_c.setrowsize(0, 2);
+      bcrsDiagonalMatrix_c.setrowsize(1, 1);
+      bcrsDiagonalMatrix_c.endrowsizes();
+
+      bcrsDiagonalMatrix_c.addindex(0, 0);
+      bcrsDiagonalMatrix_c.addindex(0, 2);
+      bcrsDiagonalMatrix_c.addindex(1, 0);
+      bcrsDiagonalMatrix_c.endindices();
+
+      bcrsDiagonalMatrix_c[0][0] = DiagMatrix{1, 4, 5};
+      bcrsDiagonalMatrix_c[0][2] = DiagMatrix{.1, .4, .5};
+      bcrsDiagonalMatrix_c[1][0] = DiagMatrix{-1, -.1, 99};
+
+      BCRSofDiagonalMatrix bcrsDiagonalMatrix_check(bcrsDiagonalMatrix_c);
+      bcrsDiagonalMatrix_check *= scalar_b;
+
+      BCRSofDiagonalMatrix bcrsDiagonalMatrix_a(bcrsDiagonalMatrix_c);
+      bcrsDiagonalMatrix_a = 0;
+
+      addProduct(bcrsDiagonalMatrix_a, scalar_b, bcrsDiagonalMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsDiagonalMatrix_check, bcrsDiagonalMatrix_a);
+
+      subtractProduct(bcrsDiagonalMatrix_a, scalar_b, bcrsDiagonalMatrix_c);
+      addProduct(bcrsDiagonalMatrix_a, scalar_b, bcrsDiagonalMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsDiagonalMatrix_check, bcrsDiagonalMatrix_a);
+
+      addProduct(bcrsDiagonalMatrix_a, 2.0, scalar_b, bcrsDiagonalMatrix_c);
+      subtractProduct(bcrsDiagonalMatrix_a, scalar_b, bcrsDiagonalMatrix_c);
+      subtractProduct(bcrsDiagonalMatrix_a, scalar_b, bcrsDiagonalMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsDiagonalMatrix_check, bcrsDiagonalMatrix_a);
+
+      addProduct(bcrsDiagonalMatrix_a, 2.0, scalar_b, bcrsDiagonalMatrix_c);
+      subtractProduct(bcrsDiagonalMatrix_a, 2.0, scalar_b,
+                      bcrsDiagonalMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsDiagonalMatrix_check, bcrsDiagonalMatrix_a);
+    }
+
+    // BLOCKED MATRICES of ARBITRARY MATRIX but with differing index sets - what
+    // do we expect?!
+
+    return passed;
+  }
+
+  bool checkProductHelperMethods_VectorMatrixVector() {
+    bool passed = true;
+
+    const Dune::FieldVector<FT, 2> fieldVector_c = {4, 2};
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_b = {{1, 2}, {3, 4}};
+
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_b = {3, 2};
+
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_b(2);
+
+    {
+      // case FV += FM*FV
+      typedef Dune::FieldVector<FT, 2> ResultType;
+      {
+        ResultType fieldVector_a(0), fieldVector_check = {8, 20};
+
+        addProduct(fieldVector_a, squareFieldMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+        subtractProduct(fieldVector_check, squareFieldMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+      }
+      // case FV += DM*FV
+      {
+        ResultType fieldVector_a(0), fieldVector_check = {12, 4};
+
+        addProduct(fieldVector_a, diagonalMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+        subtractProduct(fieldVector_check, diagonalMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+      }
+      // case FV += SM*FV
+      {
+        ResultType fieldVector_a(0), fieldVector_check = {8, 4};
+
+        addProduct(fieldVector_a, scaledIdentityMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+        subtractProduct(fieldVector_check, scaledIdentityMatrix_b,
+                        fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+      }
+    }
+
+    // RECTANGULAR MATRICES
+
+    const Dune::FieldMatrix<FT, 3, 2> rectFieldMatrix_b = {
+        {1, 2}, {3, 4}, {5, 6}};
+
+    // case FV<3> += FM<3,2>*FV<2>
+    {
+      Dune::FieldVector<FT, 3> fieldVector_a(0),
+          fieldVector_check = {8, 20, 32};
+
+      addProduct(fieldVector_a, rectFieldMatrix_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+      subtractProduct(fieldVector_check, rectFieldMatrix_b, fieldVector_c);
+      passed = passed and
+               isCloseDune(fieldVector_check, Dune::FieldVector<FT, 3>(0));
+    }
+
+    return passed;
+  }
+
+  bool checkProductHelperMethods_VectorScalarVector() {
+    bool passed = true;
+
+    const FT scalar_b = 2;
+    const double double_b = 2;
+
+    const Dune::FieldVector<FT, 2> fieldVector_c = {4, 2};
+
+    // SQUARE MATRICES
+    // case FV += scalar*FV
+    {
+      typedef Dune::FieldVector<FT, 2> ResultType;
+      ResultType fieldVector_a(0), fieldVector_check = {8, 4};
+
+      addProduct(fieldVector_a, scalar_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+      subtractProduct(fieldVector_check, scalar_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+    }
+    // case FV += double*FV
+    {
+      typedef Dune::FieldVector<FT, 2> ResultType;
+      ResultType fieldVector_a(0), fieldVector_check = {8, 4};
+
+      addProduct(fieldVector_a, double_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+      subtractProduct(fieldVector_check, double_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+    }
+
+    return passed;
+  }
+
+  bool checkProductHelperMethods_ScalarScalarScalar() {
+    bool passed = true;
+
+    const FT scalar_b = 2, scalar_c = 3;
+
+    const double double_b = 2.0, double_c = 3.0;
+
+    // SQUARE MATRICES
+    // case scalar += scalar*scalar
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 6;
+
+      addProduct(scalar_a, scalar_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, scalar_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case scalar += double*scalar
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 6;
+
+      addProduct(scalar_a, double_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, double_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case scalar += scalar*double
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 6;
+
+      addProduct(scalar_a, scalar_b, double_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, scalar_b, double_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case scalar += double*double
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 6;
+
+      addProduct(scalar_a, double_b, double_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, double_b, double_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case double += double*double
+    {
+      typedef double ResultType;
+      ResultType double_a(0), double_check = 6;
+
+      addProduct(double_a, double_b, double_c);
+      passed = passed and isCloseAbs(double_a, double_check);
+
+      subtractProduct(double_check, double_b, double_c);
+      passed = passed and isCloseAbs(double_check, ResultType(0));
+    }
+
+    return passed;
+  }
+
+  bool checkScaledProductHelperMethods_MatrixMatrixMatrix() {
+    bool passed = true;
+
+    const FT scaling = 2;
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_b = {{1, 2}, {3, 4}},
+
+                                      squareFieldMatrix_c = {{2, 4}, {3, 5}};
+
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_b = {3, 2},
+                                      diagonalMatrix_c = {2, 1};
+
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_b(2),
+        scaledIdentityMatrix_c(4);
+    // case FM += FM*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{16, 28}, {36, 64}};
+
+      addProduct(squareFieldMatrix_a, scaling, squareFieldMatrix_b,
+                 squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, squareFieldMatrix_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += DM*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{12, 24}, {12, 20}};
+
+      addProduct(squareFieldMatrix_a, scaling, diagonalMatrix_b,
+                 squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, diagonalMatrix_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += FM*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{4, 4}, {12, 8}};
+
+      addProduct(squareFieldMatrix_a, scaling, squareFieldMatrix_b,
+                 diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, squareFieldMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += DM*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{12, 0}, {0, 4}};
+
+      addProduct(squareFieldMatrix_a, scaling, diagonalMatrix_b,
+                 diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, diagonalMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += SM*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 16}, {12, 20}};
+
+      addProduct(squareFieldMatrix_a, scaling, scaledIdentityMatrix_b,
+                 squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, scaledIdentityMatrix_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += FM*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 16}, {24, 32}};
+
+      addProduct(squareFieldMatrix_a, scaling, squareFieldMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, squareFieldMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += SM*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{16, 0}, {0, 16}};
+
+      addProduct(squareFieldMatrix_a, scaling, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, scaledIdentityMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += DM*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{24, 0}, {0, 16}};
+
+      addProduct(squareFieldMatrix_a, scaling, diagonalMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, diagonalMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += SM*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 0}, {0, 4}};
+
+      addProduct(squareFieldMatrix_a, scaling, scaledIdentityMatrix_b,
+                 diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, scaledIdentityMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case DM += DM*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {12, 4};
+
+      addProduct(diagonalMatrix_a, scaling, diagonalMatrix_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, diagonalMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += DM*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {24, 16};
+
+      addProduct(diagonalMatrix_a, scaling, diagonalMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, diagonalMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += SM*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {8, 4};
+
+      addProduct(diagonalMatrix_a, scaling, scaledIdentityMatrix_b,
+                 diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, scaledIdentityMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += SM*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {16, 16};
+
+      addProduct(diagonalMatrix_a, scaling, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, scaledIdentityMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case SM += SM*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(16);
+
+      addProduct(scaledIdentityMatrix_a, scaling, scaledIdentityMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, scaling,
+                      scaledIdentityMatrix_b, scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // RECTANGULAR MATRICES
+
+    const Dune::FieldMatrix<FT, 3, 2> rectFieldMatrix_b = {
+        {1, 2}, {3, 4}, {5, 6}};
+
+    const Dune::FieldMatrix<FT, 2, 4> rectFieldMatrix_c = {{2, 3, 2, 4},
+                                                           {1, 4, 3, 5}};
+    // case FM<3,4> += FM<3,2>*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 3, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {
+              {8, 22, 16, 28}, {20, 50, 36, 64}, {32, 78, 56, 100}};
+
+      addProduct(rectFieldMatrix_a, scaling, rectFieldMatrix_b,
+                 rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, rectFieldMatrix_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<2,4> += DM<2>*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{12, 18, 12, 24}, {4, 16, 12, 20}};
+
+      addProduct(rectFieldMatrix_a, scaling, diagonalMatrix_b,
+                 rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, diagonalMatrix_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<3,2> += FM<3,2>*DM<2>
+    {
+      typedef Dune::FieldMatrix<FT, 3, 2> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{4, 4}, {12, 8}, {20, 12}};
+
+      addProduct(rectFieldMatrix_a, scaling, rectFieldMatrix_b,
+                 diagonalMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, rectFieldMatrix_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<2,4> += SM<2>*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{8, 12, 8, 16}, {4, 16, 12, 20}};
+
+      addProduct(rectFieldMatrix_a, scaling, scaledIdentityMatrix_b,
+                 rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, scaledIdentityMatrix_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    // case FM<3,2> += FM<3,2>*SM<2>
+    {
+      typedef Dune::FieldMatrix<FT, 3, 2> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{8, 16}, {24, 32}, {40, 48}};
+
+      addProduct(rectFieldMatrix_a, scaling, rectFieldMatrix_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, rectFieldMatrix_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+    return passed;
+  }
+
+  bool checkScaledProductHelperMethods_MatrixScalarMatrix() {
+    bool passed = true;
+
+    const FT scalar_b = 2;
+    const double double_b = 2.0;
+    const FT scaling = 2;
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_c = {{2, 4}, {3, 5}};
+
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_c = {2, 1};
+
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_c(4);
+
+    // case FM += scalar*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 16}, {12, 20}};
+
+      addProduct(squareFieldMatrix_a, scaling, scalar_b, squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, scalar_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += scalar*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 0}, {0, 4}};
+
+      addProduct(squareFieldMatrix_a, scaling, scalar_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, scalar_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += scalar*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{16, 0}, {0, 16}};
+
+      addProduct(squareFieldMatrix_a, scaling, scalar_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, scalar_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case DM += scalar*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {8, 4};
+
+      addProduct(diagonalMatrix_a, scaling, scalar_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, scalar_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += scalar*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {16, 16};
+
+      addProduct(diagonalMatrix_a, scaling, scalar_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, scalar_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case SM += scalar*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(16);
+
+      addProduct(scaledIdentityMatrix_a, scaling, scalar_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, scaling, scalar_b,
+                      scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // RECTANGULAR MATRICES
+
+    const Dune::FieldMatrix<FT, 2, 4> rectFieldMatrix_c = {{2, 3, 2, 4},
+                                                           {1, 4, 3, 5}};
+    // case FM<2,4> += scalar*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{8, 12, 8, 16}, {4, 16, 12, 20}};
+
+      addProduct(rectFieldMatrix_a, scaling, scalar_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, scalar_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+
+    // and again for fixed scalar type double
+    // SQUARE MATRICES
+    // case FM += double*FM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 16}, {12, 20}};
+
+      addProduct(squareFieldMatrix_a, scaling, double_b, squareFieldMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, double_b,
+                      squareFieldMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += double*DM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{8, 0}, {0, 4}};
+
+      addProduct(squareFieldMatrix_a, scaling, double_b, diagonalMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, double_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case FM += double*SM
+    {
+      typedef Dune::FieldMatrix<FT, 2, 2> ResultType;
+      ResultType squareFieldMatrix_a(0),
+          squareFieldMatrix_check = {{16, 0}, {0, 16}};
+
+      addProduct(squareFieldMatrix_a, scaling, double_b,
+                 scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(squareFieldMatrix_a, squareFieldMatrix_check);
+
+      subtractProduct(squareFieldMatrix_check, scaling, double_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(squareFieldMatrix_check, ResultType(0));
+    }
+    // case DM += double*DM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {8, 4};
+
+      addProduct(diagonalMatrix_a, scaling, double_b, diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, double_b,
+                      diagonalMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case DM += double*SM
+    {
+      typedef Dune::DiagonalMatrix<FT, 2> ResultType;
+      ResultType diagonalMatrix_a(0), diagonalMatrix_check = {16, 16};
+
+      addProduct(diagonalMatrix_a, scaling, double_b, scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_a, diagonalMatrix_check);
+
+      subtractProduct(diagonalMatrix_check, scaling, double_b,
+                      scaledIdentityMatrix_c);
+      passed = passed and isCloseDune(diagonalMatrix_check, ResultType(0));
+    }
+    // case SM += double*SM
+    {
+      typedef Dune::ScaledIdentityMatrix<FT, 2> ResultType;
+      ResultType scaledIdentityMatrix_a(0), scaledIdentityMatrix_check(16);
+
+      addProduct(scaledIdentityMatrix_a, scaling, double_b,
+                 scaledIdentityMatrix_c);
+      passed = passed and
+               isCloseDune(scaledIdentityMatrix_a, scaledIdentityMatrix_check);
+
+      subtractProduct(scaledIdentityMatrix_check, scaling, double_b,
+                      scaledIdentityMatrix_c);
+      passed =
+          passed and isCloseDune(scaledIdentityMatrix_check, ResultType(0));
+    }
+
+    // RECTANGULAR MATRICES
+    // case FM<2,4> += double*FM<2,4>
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> ResultType;
+      ResultType rectFieldMatrix_a(0),
+          rectFieldMatrix_check = {{8, 12, 8, 16}, {4, 16, 12, 20}};
+
+      addProduct(rectFieldMatrix_a, scaling, double_b, rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_a, rectFieldMatrix_check);
+
+      subtractProduct(rectFieldMatrix_check, scaling, double_b,
+                      rectFieldMatrix_c);
+      passed = passed and isCloseDune(rectFieldMatrix_check, ResultType(0));
+    }
+
+    // BLOCKED MATRICES
+    {
+      typedef Dune::FieldMatrix<FT, 2, 4> RectFieldMatrix;
+      typedef Dune::BCRSMatrix<RectFieldMatrix> BCRSofRectFieldMatrix;
+
+      BCRSofRectFieldMatrix bcrsRectFieldMatrix_c(
+          2, 3, BCRSofRectFieldMatrix::random);
+      bcrsRectFieldMatrix_c.setrowsize(0, 2);
+      bcrsRectFieldMatrix_c.setrowsize(1, 1);
+      bcrsRectFieldMatrix_c.endrowsizes();
+
+      bcrsRectFieldMatrix_c.addindex(0, 0);
+      bcrsRectFieldMatrix_c.addindex(0, 2);
+      bcrsRectFieldMatrix_c.addindex(1, 0);
+      bcrsRectFieldMatrix_c.endindices();
+
+      bcrsRectFieldMatrix_c[0][0] = RectFieldMatrix{{1, 2, 3, 4}, {3, 4, 5, 6}};
+      bcrsRectFieldMatrix_c[0][2] = RectFieldMatrix{{2, 2, 3, 3}, {1, 1, 4, 2}};
+      bcrsRectFieldMatrix_c[1][0] = RectFieldMatrix{{3, 1, 2, 5}, {4, 3, 2, 3}};
+
+      BCRSofRectFieldMatrix bcrsRectFieldMatrix_check(bcrsRectFieldMatrix_c);
+      bcrsRectFieldMatrix_check *= scaling * scalar_b;
+
+      BCRSofRectFieldMatrix bcrsRectFieldMatrix_a(bcrsRectFieldMatrix_c),
+          bcrsRectFieldMatrix_zero(bcrsRectFieldMatrix_c);
+      bcrsRectFieldMatrix_a = bcrsRectFieldMatrix_zero = 0;
+
+      addProduct(bcrsRectFieldMatrix_a, scaling, scalar_b,
+                 bcrsRectFieldMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsRectFieldMatrix_check, bcrsRectFieldMatrix_a);
+
+      subtractProduct(bcrsRectFieldMatrix_check, scaling, scalar_b,
+                      bcrsRectFieldMatrix_c);
+      passed = passed and
+               isCloseDune(bcrsRectFieldMatrix_check, bcrsRectFieldMatrix_zero);
+    }
+    return passed;
+  }
+
+  bool checkScaledProductHelperMethods_VectorMatrixVector() {
+    bool passed = true;
+
+    const FT scaling = 2;
+
+    const Dune::FieldVector<FT, 2> fieldVector_c = {4, 2};
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FT, 2, 2> squareFieldMatrix_b = {{1, 2}, {3, 4}};
+
+    const Dune::DiagonalMatrix<FT, 2> diagonalMatrix_b = {3, 2};
+
+    const Dune::ScaledIdentityMatrix<FT, 2> scaledIdentityMatrix_b(2);
+
+    {
+      // case FV += FM*FV
+      typedef Dune::FieldVector<FT, 2> ResultType;
+      {
+        ResultType fieldVector_a(0), fieldVector_check = {16, 40};
+
+        addProduct(fieldVector_a, scaling, squareFieldMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+        subtractProduct(fieldVector_check, scaling, squareFieldMatrix_b,
+                        fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+      }
+      // case FV += DM*FV
+      {
+        ResultType fieldVector_a(0), fieldVector_check = {24, 8};
+
+        addProduct(fieldVector_a, scaling, diagonalMatrix_b, fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+        subtractProduct(fieldVector_check, scaling, diagonalMatrix_b,
+                        fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+      }
+      // case FV += SM*FV
+      {
+        ResultType fieldVector_a(0), fieldVector_check = {16, 8};
+
+        addProduct(fieldVector_a, scaling, scaledIdentityMatrix_b,
+                   fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+        subtractProduct(fieldVector_check, scaling, scaledIdentityMatrix_b,
+                        fieldVector_c);
+        passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+      }
+    }
+
+    // RECTANGULAR MATRICES
+
+    const Dune::FieldMatrix<FT, 3, 2> rectFieldMatrix_b = {
+        {1, 2}, {3, 4}, {5, 6}};
+
+    // case FV<3> += FM<3,2>*FV<2>
+    {
+      Dune::FieldVector<FT, 3> fieldVector_a(0),
+          fieldVector_check = {16, 40, 64};
+
+      addProduct(fieldVector_a, scaling, rectFieldMatrix_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+      subtractProduct(fieldVector_check, scaling, rectFieldMatrix_b,
+                      fieldVector_c);
+      passed = passed and
+               isCloseDune(fieldVector_check, Dune::FieldVector<FT, 3>(0));
+    }
+
+    return passed;
+  }
+
+  bool checkScaledProductHelperMethods_VectorScalarVector() {
+    bool passed = true;
+
+    const FT scalar_b = 2;
+    const double double_b = 2;
+    const FT scaling = 2;
+
+    const Dune::FieldVector<FT, 2> fieldVector_c = {4, 2};
+
+    // SQUARE MATRICES
+    // case FV += scalar*FV
+    {
+      typedef Dune::FieldVector<FT, 2> ResultType;
+      ResultType fieldVector_a(0), fieldVector_check = {16, 8};
+
+      addProduct(fieldVector_a, scaling, scalar_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+      subtractProduct(fieldVector_check, scaling, scalar_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+    }
+    // case FV += double*FV
+    {
+      typedef Dune::FieldVector<FT, 2> ResultType;
+      ResultType fieldVector_a(0), fieldVector_check = {16, 8};
+
+      addProduct(fieldVector_a, scaling, double_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_a, fieldVector_check);
+
+      subtractProduct(fieldVector_check, scaling, double_b, fieldVector_c);
+      passed = passed and isCloseDune(fieldVector_check, ResultType(0));
+    }
+
+    return passed;
+  }
+
+  bool checkScaledProductHelperMethods_ScalarScalarScalar() {
+    bool passed = true;
+
+    const FT scalar_b = 2, scalar_c = 3;
+
+    const double double_b = 2.0, double_c = 3.0;
+
+    const FT scaling = 2;
+
+    // SQUARE MATRICES
+    // case scalar += scalar*scalar
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 12;
+
+      addProduct(scalar_a, scaling, scalar_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, scaling, scalar_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case scalar += double*scalar
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 12;
+
+      addProduct(scalar_a, scaling, double_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, scaling, double_b, scalar_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case scalar += scalar*double
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 12;
+
+      addProduct(scalar_a, scaling, scalar_b, double_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, scaling, scalar_b, double_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case scalar += double*double
+    {
+      typedef FT ResultType;
+      ResultType scalar_a(0), scalar_check = 12;
+
+      addProduct(scalar_a, scaling, double_b, double_c);
+      passed = passed and isCloseAbs(scalar_a, scalar_check);
+
+      subtractProduct(scalar_check, scaling, double_b, double_c);
+      passed = passed and isCloseAbs(scalar_check, ResultType(0));
+    }
+    // case double += double*double
+    {
+      typedef double ResultType;
+      ResultType double_a(0), double_check = 12;
+
+      addProduct(double_a, scaling, double_b, double_c);
+      passed = passed and isCloseAbs(double_a, double_check);
+
+      subtractProduct(double_check, scaling, double_b, double_c);
+      passed = passed and isCloseAbs(double_check, ResultType(0));
+    }
+
+    return passed;
+  }
+
+  bool checkBCRSSubPattern() {
+    typedef Dune::FieldMatrix<FT, 1, 1> M;
+
+    Dune::BCRSMatrix<M> B(3, 3, Dune::BCRSMatrix<M>::random);
+    {
+      B.setrowsize(0, 3);
+      B.setrowsize(1, 2);
+      B.setrowsize(2, 1);
+      B.endrowsizes();
+
+      B.addindex(0, 0);
+      B.addindex(0, 1);
+      B.addindex(0, 2);
+      B.addindex(1, 1);
+      B.addindex(1, 2);
+      B.addindex(2, 2);
+      B.endindices();
+
+      B[0][0] = 10;
+      B[0][1] = 20;
+      B[0][2] = 30;
+      B[1][1] = 4000;
+      B[1][2] = 5000;
+      B[2][2] = 60;
+    }
+
+    Dune::BCRSMatrix<M> A(3, 3, Dune::BCRSMatrix<M>::random);
+    {
+      A.setrowsize(0, 3);
+      A.setrowsize(1, 0);
+      A.setrowsize(2, 1);
+      A.endrowsizes();
+
+      A.addindex(0, 0);
+      A.addindex(0, 1);
+      A.addindex(0, 2);
+      A.addindex(2, 2);
+      A.endindices();
+
+      A[0][0] = B[0][0];
+      A[0][1] = B[0][1];
+      A[0][2] = B[0][2];
+      A[2][2] = B[2][2];
+    }
+
+    FT b(2);
+    Dune::BCRSMatrix<M> C(3, 3, Dune::BCRSMatrix<M>::random);
+    {
+      C.setrowsize(0, 3);
+      C.setrowsize(1, 2);
+      C.setrowsize(2, 1);
+      C.endrowsizes();
+
+      C.addindex(0, 0);
+      C.addindex(0, 1);
+      C.addindex(0, 2);
+      C.addindex(1, 1);
+      C.addindex(1, 2);
+      C.addindex(2, 2);
+      C.endindices();
+
+      C[0][0] = (b + 1) * B[0][0];
+      C[0][1] = (b + 1) * B[0][1];
+      C[0][2] = (b + 1) * B[0][2];
+      C[1][1] = B[1][1];
+      C[1][2] = B[1][2];
+      C[2][2] = (b + 1) * B[2][2];
+    }
+
+    addProduct(B, b, A);
+
+    return isCloseDune(B, C);
+  }
+};
+
+int main(int argc, char* argv[]) {
+  Dune::MPIHelper::instance(argc, argv);
+  bool passed = true;
+  {
+    ArithmeticTestSuite<double> testSuite;
+    passed &= testSuite.check();
+  }
+  {
+    ArithmeticTestSuite<long double> testSuite;
+    passed &= testSuite.check();
+  }
+  {
+    ArithmeticTestSuite<float> testSuite;
+    passed &= testSuite.check();
+  }
+  return passed ? 0 : 1;
+}
diff --git a/dune/matrix-vector/test/common.hh b/dune/matrix-vector/test/common.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1036d57f7b2223101bc7cd27ab7e869ab462379e
--- /dev/null
+++ b/dune/matrix-vector/test/common.hh
@@ -0,0 +1,24 @@
+#ifndef DUNE_MATRIX_VECTOR_TEST_COMMON_HH
+#define DUNE_MATRIX_VECTOR_TEST_COMMON_HH
+
+template <typename FT>
+struct ToleranceTraits {
+  static constexpr FT tol = 50 * std::numeric_limits<FT>::epsilon();
+};
+
+template <class T>
+typename Dune::FieldTraits<typename T::field_type>::real_type diffDune(
+    T t1, const T& t2) {
+  t1 -= t2;
+  return t1.infinity_norm();
+}
+template <class TX>
+bool isCloseDune(TX t1, const TX& t2) {
+  return diffDune(t1, t2) < ToleranceTraits<typename TX::field_type>::tol;
+}
+template <class FT, class Other>
+bool isCloseAbs(FT t1, Other t2) {
+  return std::abs<FT>(t1 - t2) < ToleranceTraits<FT>::tol;
+}
+
+#endif
diff --git a/dune/matrix-vector/test/staticmatrixtoolstest.cc b/dune/matrix-vector/test/staticmatrixtoolstest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..368f523445f45f9a1b3f9baca94d29252c7fc7f4
--- /dev/null
+++ b/dune/matrix-vector/test/staticmatrixtoolstest.cc
@@ -0,0 +1,884 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/diagonalmatrix.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/matrix.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+#include "../addtodiagonal.hh"
+#include "../transformmatrix.hh"
+
+using namespace Dune::MatrixVector;
+
+template <class FieldType>
+class StaticMatrixToolsTestSuite {
+public:
+  bool check() {
+    bool passed = true;
+
+    passed = passed and checkAddToDiagonal();
+    passed = passed and checkAddTransformedMatrix();
+    passed = passed and checkTransformMatrix();
+
+    return passed;
+  }
+
+private:
+  bool checkAddToDiagonal() {
+    bool passed = true;
+
+    const FieldType scalar_a = 3;
+    const double double_a = 3;
+
+    // SQUARE MATRICES
+
+    // case FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_x = {{1, 2}, {3, 4}},
+
+                                         squareFieldMatrix_check = {{4, 2},
+                                                                    {3, 7}};
+
+      addToDiagonal(squareFieldMatrix_x, scalar_a);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_x, squareFieldMatrix_check) < 1e-12;
+    }
+    // case DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_x = {3, 2},
+                                         diagonalMatrix_check = {6, 5};
+
+      addToDiagonal(diagonalMatrix_x, scalar_a);
+
+      passed =
+          passed and myDiff(diagonalMatrix_x, diagonalMatrix_check) < 1e-12;
+    }
+    // case SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_x(2),
+          scaledIdentityMatrix_check(5);
+
+      addToDiagonal(scaledIdentityMatrix_x, scalar_a);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_x, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case with hardwired scalar type as double (as occurring in
+    // localassemblers, see also the commment in staticmatrixtools.hh)
+    // why does this even compile?
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_x(2),
+          scaledIdentityMatrix_check(5);
+
+      addToDiagonal(scaledIdentityMatrix_x, double_a);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_x, scaledIdentityMatrix_check) < 1e-12;
+    }
+
+    // RECTANGULAR MATRICES
+    {
+      Dune::FieldMatrix<FieldType, 2, 4> rectFieldMatrix_x = {{1, 2, 3, 4},
+                                                              {5, 6, 7, 8}},
+
+                                         rectFieldMatrix_check = {{4, 2, 3, 4},
+                                                                  {5, 9, 7, 8}};
+
+      addToDiagonal(rectFieldMatrix_x, scalar_a);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_x, rectFieldMatrix_check) < 1e-12;
+    }
+
+    return passed;
+  }
+
+  bool checkAddTransformedMatrix() {
+    bool passed = true;
+
+    const FieldType scalar_T1 = 3, scalar_T2 = 4, scalar_B = 2;
+
+    const double double_T1 = 3, double_T2 = 4, double_B = 2;
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_T1 = {{1, 2},
+                                                                     {3, 4}},
+                                             squareFieldMatrix_T2 = {{2, 3},
+                                                                     {1, 4}},
+                                             squareFieldMatrix_B = {{3, 4},
+                                                                    {5, 2}};
+
+    const Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_T1 = {1, 2},
+                                             diagonalMatrix_T2 = {2, 3},
+                                             diagonalMatrix_B = {3, 2};
+
+    const Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_T1(3),
+        scaledIdentityMatrix_T2(4), scaledIdentityMatrix_B(2);
+
+    // test default imp by case Matrix += FM^t*Matrix*FM
+    {
+      Dune::Matrix<Dune::FieldMatrix<FieldType, 1, 1>> squareMatrix_A(2, 2),
+          squareMatrix_B(2, 2), squareMatrix_check(2, 2);
+      squareMatrix_A = squareMatrix_B = squareMatrix_check = 0;
+
+      squareMatrix_B[0][0] = 3;
+      squareMatrix_B[0][1] = 4;
+      squareMatrix_B[1][0] = 5;
+      squareMatrix_B[1][1] = 2;
+
+      squareMatrix_check[0][0] = 46;
+      squareMatrix_check[0][1] = 94;
+      squareMatrix_check[1][0] = 68;
+      squareMatrix_check[1][1] = 142;
+
+      addTransformedMatrix(squareMatrix_A, squareFieldMatrix_T1, squareMatrix_B,
+                           squareFieldMatrix_T2);
+
+      passed = passed and myDiff(squareMatrix_A, squareMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*FM*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{46, 94}, {68, 142}};
+
+      addTransformedMatrix(squareFieldMatrix_A, squareFieldMatrix_T1,
+                           squareFieldMatrix_B, squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*DM*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{12, 33}, {20, 50}};
+
+      addTransformedMatrix(squareFieldMatrix_A, squareFieldMatrix_T1,
+                           diagonalMatrix_B, squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += DM^t*FM*DM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{6, 12}, {20, 12}};
+
+      addTransformedMatrix(squareFieldMatrix_A, diagonalMatrix_T1,
+                           squareFieldMatrix_B, diagonalMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += DM^t*DM*DM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{6, 0}, {0, 12}};
+
+      addTransformedMatrix(squareFieldMatrix_A, diagonalMatrix_T1,
+                           diagonalMatrix_B, diagonalMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*SM*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{10, 30}, {16, 44}};
+
+      addTransformedMatrix(squareFieldMatrix_A, squareFieldMatrix_T1,
+                           scaledIdentityMatrix_B, squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += SM^t*FM*SM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{36, 48}, {60, 24}};
+
+      addTransformedMatrix(squareFieldMatrix_A, scaledIdentityMatrix_T1,
+                           squareFieldMatrix_B, scaledIdentityMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*field_type*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{10, 30}, {16, 44}};
+
+      addTransformedMatrix(squareFieldMatrix_A, squareFieldMatrix_T1, scalar_B,
+                           squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += field_type*FM*field_type
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{36, 48}, {60, 24}};
+
+      addTransformedMatrix(squareFieldMatrix_A, scalar_T1, squareFieldMatrix_B,
+                           scalar_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*double*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{10, 30}, {16, 44}};
+
+      addTransformedMatrix(squareFieldMatrix_A, squareFieldMatrix_T1, double_B,
+                           squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += double*FM*double
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{36, 48}, {60, 24}};
+
+      addTransformedMatrix(squareFieldMatrix_A, double_T1, squareFieldMatrix_B,
+                           double_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*DM*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {6, 12};
+
+      addTransformedMatrix(diagonalMatrix_A, diagonalMatrix_T1,
+                           diagonalMatrix_B, diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*SM*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {4, 12};
+
+      addTransformedMatrix(diagonalMatrix_A, diagonalMatrix_T1,
+                           scaledIdentityMatrix_B, diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += SM^t*DM*SM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {36, 24};
+
+      addTransformedMatrix(diagonalMatrix_A, scaledIdentityMatrix_T1,
+                           diagonalMatrix_B, scaledIdentityMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*field_type*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {4, 12};
+
+      addTransformedMatrix(diagonalMatrix_A, diagonalMatrix_T1, scalar_B,
+                           diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += field_type*DM*field_type
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {36, 24};
+
+      addTransformedMatrix(diagonalMatrix_A, scalar_T1, diagonalMatrix_B,
+                           scalar_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*double*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {4, 12};
+
+      addTransformedMatrix(diagonalMatrix_A, diagonalMatrix_T1, double_B,
+                           diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += double*DM*double
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {36, 24};
+
+      addTransformedMatrix(diagonalMatrix_A, double_T1, diagonalMatrix_B,
+                           double_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case SM += SM^t*SM*SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      addTransformedMatrix(scaledIdentityMatrix_A, scaledIdentityMatrix_T1,
+                           scaledIdentityMatrix_B, scaledIdentityMatrix_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += SM^t*field_type*SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      addTransformedMatrix(scaledIdentityMatrix_A, scaledIdentityMatrix_T1,
+                           scalar_B, scaledIdentityMatrix_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += field_type*SM*field_type
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      addTransformedMatrix(scaledIdentityMatrix_A, scalar_T1,
+                           scaledIdentityMatrix_B, scalar_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += SM^t*double*SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      addTransformedMatrix(scaledIdentityMatrix_A, scaledIdentityMatrix_T1,
+                           double_B, scaledIdentityMatrix_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += double*SM*double
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      addTransformedMatrix(scaledIdentityMatrix_A, double_T1,
+                           scaledIdentityMatrix_B, double_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case field_type += field_type*field_type*field_type
+    {
+      FieldType scalar_A(0), scalar_check = 24;
+
+      addTransformedMatrix(scalar_A, scalar_T1, scalar_B, scalar_T2);
+
+      passed = passed and myDiff(scalar_A, scalar_check) < 1e-12;
+    }
+
+    // RECTANGULAR MATRICES
+    const Dune::FieldMatrix<FieldType, 2, 4> rectFieldMatrix_T1 = {{1, 2, 3, 4},
+                                                                   {5, 6, 7,
+                                                                    8}},
+                                             rectFieldMatrix_T2 = {
+                                                 {2, 3, 4, 5}, {6, 7, 8, 9}};
+
+    // case FM += FM^t*FM*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{140, 182, 224, 266},
+                                   {192, 248, 304, 360},
+                                   {244, 314, 384, 454},
+                                   {296, 380, 464, 548}};
+
+      addTransformedMatrix(rectFieldMatrix_A, rectFieldMatrix_T1,
+                           squareFieldMatrix_B, rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*DM*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{66, 79, 92, 105},
+                                   {84, 102, 120, 138},
+                                   {102, 125, 148, 171},
+                                   {120, 148, 176, 204}};
+
+      addTransformedMatrix(rectFieldMatrix_A, rectFieldMatrix_T1,
+                           diagonalMatrix_B, rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*SM*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{64, 76, 88, 100},
+                                   {80, 96, 112, 128},
+                                   {96, 116, 136, 156},
+                                   {112, 136, 160, 184}};
+
+      addTransformedMatrix(rectFieldMatrix_A, rectFieldMatrix_T1,
+                           scaledIdentityMatrix_B, rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*field_type*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{64, 76, 88, 100},
+                                   {80, 96, 112, 128},
+                                   {96, 116, 136, 156},
+                                   {112, 136, 160, 184}};
+
+      addTransformedMatrix(rectFieldMatrix_A, rectFieldMatrix_T1, scalar_B,
+                           rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*double*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{64, 76, 88, 100},
+                                   {80, 96, 112, 128},
+                                   {96, 116, 136, 156},
+                                   {112, 136, 160, 184}};
+
+      addTransformedMatrix(rectFieldMatrix_A, rectFieldMatrix_T1, double_B,
+                           rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+
+    return passed;
+  }
+
+  bool checkTransformMatrix() {
+    bool passed = true;
+
+    const FieldType scalar_T1 = 3, scalar_T2 = 4, scalar_B = 2;
+
+    const double double_T1 = 3, double_T2 = 4, double_B = 2;
+
+    // SQUARE MATRICES
+    const Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_T1 = {{1, 2},
+                                                                     {3, 4}},
+                                             squareFieldMatrix_T2 = {{2, 3},
+                                                                     {1, 4}},
+                                             squareFieldMatrix_B = {{3, 4},
+                                                                    {5, 2}};
+
+    const Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_T1 = {1, 2},
+                                             diagonalMatrix_T2 = {2, 3},
+                                             diagonalMatrix_B = {3, 2};
+
+    const Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_T1(3),
+        scaledIdentityMatrix_T2(4), scaledIdentityMatrix_B(2);
+
+    // test default imp by case Matrix += FM^t*Matrix*FM
+    {
+      Dune::Matrix<Dune::FieldMatrix<FieldType, 1, 1>> squareMatrix_A(2, 2),
+          squareMatrix_B(2, 2), squareMatrix_check(2, 2);
+      squareMatrix_A = squareMatrix_B = squareMatrix_check = 0;
+
+      squareMatrix_B[0][0] = 3;
+      squareMatrix_B[0][1] = 4;
+      squareMatrix_B[1][0] = 5;
+      squareMatrix_B[1][1] = 2;
+
+      squareMatrix_check[0][0] = 46;
+      squareMatrix_check[0][1] = 94;
+      squareMatrix_check[1][0] = 68;
+      squareMatrix_check[1][1] = 142;
+
+      transformMatrix(squareMatrix_A, squareFieldMatrix_T1, squareMatrix_B,
+                      squareFieldMatrix_T2);
+
+      passed = passed and myDiff(squareMatrix_A, squareMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*FM*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{46, 94}, {68, 142}};
+
+      transformMatrix(squareFieldMatrix_A, squareFieldMatrix_T1,
+                      squareFieldMatrix_B, squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*DM*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{12, 33}, {20, 50}};
+
+      transformMatrix(squareFieldMatrix_A, squareFieldMatrix_T1,
+                      diagonalMatrix_B, squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += DM^t*FM*DM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{6, 12}, {20, 12}};
+
+      transformMatrix(squareFieldMatrix_A, diagonalMatrix_T1,
+                      squareFieldMatrix_B, diagonalMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += DM^t*DM*DM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{6, 0}, {0, 12}};
+
+      transformMatrix(squareFieldMatrix_A, diagonalMatrix_T1, diagonalMatrix_B,
+                      diagonalMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*SM*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{10, 30}, {16, 44}};
+
+      transformMatrix(squareFieldMatrix_A, squareFieldMatrix_T1,
+                      scaledIdentityMatrix_B, squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += SM^t*FM*SM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{36, 48}, {60, 24}};
+
+      transformMatrix(squareFieldMatrix_A, scaledIdentityMatrix_T1,
+                      squareFieldMatrix_B, scaledIdentityMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*field_type*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{10, 30}, {16, 44}};
+
+      transformMatrix(squareFieldMatrix_A, squareFieldMatrix_T1, scalar_B,
+                      squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += field_type*FM*field_type
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{36, 48}, {60, 24}};
+
+      transformMatrix(squareFieldMatrix_A, scalar_T1, squareFieldMatrix_B,
+                      scalar_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*double*FM
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{10, 30}, {16, 44}};
+
+      transformMatrix(squareFieldMatrix_A, squareFieldMatrix_T1, double_B,
+                      squareFieldMatrix_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case FM += double*FM*double
+    {
+      Dune::FieldMatrix<FieldType, 2, 2> squareFieldMatrix_A(0),
+          squareFieldMatrix_check = {{36, 48}, {60, 24}};
+
+      transformMatrix(squareFieldMatrix_A, double_T1, squareFieldMatrix_B,
+                      double_T2);
+
+      passed = passed and
+               myDiff(squareFieldMatrix_A, squareFieldMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*DM*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {6, 12};
+
+      transformMatrix(diagonalMatrix_A, diagonalMatrix_T1, diagonalMatrix_B,
+                      diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*SM*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {4, 12};
+
+      transformMatrix(diagonalMatrix_A, diagonalMatrix_T1,
+                      scaledIdentityMatrix_B, diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += SM^t*DM*SM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {36, 24};
+
+      transformMatrix(diagonalMatrix_A, scaledIdentityMatrix_T1,
+                      diagonalMatrix_B, scaledIdentityMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*field_type*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {4, 12};
+
+      transformMatrix(diagonalMatrix_A, diagonalMatrix_T1, scalar_B,
+                      diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += field_type*DM*field_type
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {36, 24};
+
+      transformMatrix(diagonalMatrix_A, scalar_T1, diagonalMatrix_B, scalar_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += DM^t*double*DM
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {4, 12};
+
+      transformMatrix(diagonalMatrix_A, diagonalMatrix_T1, double_B,
+                      diagonalMatrix_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case DM += double*DM*double
+    {
+      Dune::DiagonalMatrix<FieldType, 2> diagonalMatrix_A(0),
+          diagonalMatrix_check = {36, 24};
+
+      transformMatrix(diagonalMatrix_A, double_T1, diagonalMatrix_B, double_T2);
+
+      passed =
+          passed and myDiff(diagonalMatrix_A, diagonalMatrix_check) < 1e-12;
+    }
+    // case SM += SM^t*SM*SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      transformMatrix(scaledIdentityMatrix_A, scaledIdentityMatrix_T1,
+                      scaledIdentityMatrix_B, scaledIdentityMatrix_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += SM^t*field_type*SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      transformMatrix(scaledIdentityMatrix_A, scaledIdentityMatrix_T1, scalar_B,
+                      scaledIdentityMatrix_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += field_type*SM*field_type
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      transformMatrix(scaledIdentityMatrix_A, scalar_T1, scaledIdentityMatrix_B,
+                      scalar_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += SM^t*double*SM
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      transformMatrix(scaledIdentityMatrix_A, scaledIdentityMatrix_T1, double_B,
+                      scaledIdentityMatrix_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case SM += double*SM*double
+    {
+      Dune::ScaledIdentityMatrix<FieldType, 2> scaledIdentityMatrix_A(0),
+          scaledIdentityMatrix_check(24);
+
+      transformMatrix(scaledIdentityMatrix_A, double_T1, scaledIdentityMatrix_B,
+                      double_T2);
+
+      passed =
+          passed and
+          myDiff(scaledIdentityMatrix_A, scaledIdentityMatrix_check) < 1e-12;
+    }
+    // case field_type += field_type*field_type*field_type
+    {
+      FieldType scalar_A(0), scalar_check = 24;
+
+      transformMatrix(scalar_A, scalar_T1, scalar_B, scalar_T2);
+
+      passed = passed and myDiff(scalar_A, scalar_check) < 1e-12;
+    }
+
+    // RECTANGULAR MATRICES
+    const Dune::FieldMatrix<FieldType, 2, 4> rectFieldMatrix_T1 = {{1, 2, 3, 4},
+                                                                   {5, 6, 7,
+                                                                    8}},
+                                             rectFieldMatrix_T2 = {
+                                                 {2, 3, 4, 5}, {6, 7, 8, 9}};
+
+    // case FM += FM^t*FM*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{140, 182, 224, 266},
+                                   {192, 248, 304, 360},
+                                   {244, 314, 384, 454},
+                                   {296, 380, 464, 548}};
+
+      transformMatrix(rectFieldMatrix_A, rectFieldMatrix_T1,
+                      squareFieldMatrix_B, rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*DM*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{66, 79, 92, 105},
+                                   {84, 102, 120, 138},
+                                   {102, 125, 148, 171},
+                                   {120, 148, 176, 204}};
+
+      transformMatrix(rectFieldMatrix_A, rectFieldMatrix_T1, diagonalMatrix_B,
+                      rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*SM*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{64, 76, 88, 100},
+                                   {80, 96, 112, 128},
+                                   {96, 116, 136, 156},
+                                   {112, 136, 160, 184}};
+
+      transformMatrix(rectFieldMatrix_A, rectFieldMatrix_T1,
+                      scaledIdentityMatrix_B, rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*field_type*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{64, 76, 88, 100},
+                                   {80, 96, 112, 128},
+                                   {96, 116, 136, 156},
+                                   {112, 136, 160, 184}};
+
+      transformMatrix(rectFieldMatrix_A, rectFieldMatrix_T1, scalar_B,
+                      rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+    // case FM += FM^t*double*FM
+    {
+      Dune::FieldMatrix<FieldType, 4, 4> rectFieldMatrix_A(0),
+          rectFieldMatrix_check = {{64, 76, 88, 100},
+                                   {80, 96, 112, 128},
+                                   {96, 116, 136, 156},
+                                   {112, 136, 160, 184}};
+
+      transformMatrix(rectFieldMatrix_A, rectFieldMatrix_T1, double_B,
+                      rectFieldMatrix_T2);
+
+      passed =
+          passed and myDiff(rectFieldMatrix_A, rectFieldMatrix_check) < 1e-12;
+    }
+
+    return passed;
+  }
+
+  template <class T>
+  double myDiff(T t1, const T& t2) {
+    t1 -= t2;
+    return t1.infinity_norm();
+  }
+
+  double myDiff(FieldType t1, const FieldType& t2) { return std::abs(t1 - t2); }
+};
+
+int main(int argc, char* argv[]) {
+  Dune::MPIHelper::instance(argc, argv);
+
+  {
+    StaticMatrixToolsTestSuite<double> testSuite;
+    testSuite.check();
+  }
+  {
+    StaticMatrixToolsTestSuite<float> testSuite;
+    testSuite.check();
+  }
+  {
+    StaticMatrixToolsTestSuite<long double> testSuite;
+    testSuite.check();
+  }
+}