Implement matrix-matrix multiplication for FieldMatrix

CHANGELOG.md

@@ -40,6 +40,7 @@
   - Matrix = Scalar * Matrix
   - Matrix = Matrix * Scalar
   - Matrix = Matrix / Scalar
+  - Matrix = Matrix * Matrix
   Note that the operators
   - Vector = Vector + Vector
   - Vector = Vector - Vector

dune/common/fmatrix.hh

@@ -193,6 +193,25 @@ namespace Dune
       return result;
     }
 
+    /** \brief Matrix-matrix multiplication
+     */
+    template <class OtherScalar, int otherCols>
+    friend auto operator* ( const FieldMatrix& matrixA,
+                            const FieldMatrix<OtherScalar, COLS, otherCols>& matrixB)
+    {
+      FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,ROWS,otherCols> result;
+
+      for (size_type i = 0; i < matrixA.mat_rows(); ++i)
+        for (size_type j = 0; j < matrixB.mat_cols(); ++j)
+        {
+          result[i][j] = 0;
+          for (size_type k = 0; k < matrixA.mat_cols(); ++k)
+            result[i][j] += matrixA[i][k] * matrixB[k][j];
+        }
+
+      return result;
+    }
+
     //! Multiplies M from the left to this matrix, this matrix is not modified
     template<int l>
     FieldMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
@@ -398,6 +417,20 @@ namespace Dune
 
     //===== solve
 
+        /** \brief Matrix-matrix multiplication
+     */
+    template <class OtherScalar, int otherCols>
+    friend auto operator* ( const FieldMatrix& matrixA,
+                            const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
+    {
+      FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
+
+      for (size_type j = 0; j < matrixB.mat_cols(); ++j)
+        result[0][j] = matrixA[0][0] * matrixB[0][j];
+
+      return result;
+    }
+
     //! Multiplies M from the left to this matrix, this matrix is not modified
     template<int l>
     FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const

dune/common/test/fmatrixtest.cc

@@ -469,6 +469,21 @@ void test_matrix()
         DUNE_THROW(FMatrixError, "Return value of operator-(matrix) incorrect!");
     }
 
+    // Matrix * Matrix
+    {
+      auto transposed = [](const FM& A)
+      {
+        FieldMatrix<typename FM::field_type,FM::cols,FM::rows> AT;
+        for (int i=0; i<AT.rows; i++)
+          for (int j=0; j<AT.cols; j++)
+            AT[i][j] = A[j][i];
+
+        return AT;
+      };
+
+      DUNE_UNUSED auto product = transposed(A) * A;
+    }
+
   }
   {
     using std::abs;