Move matrix-matrix multiplication methods into a separate file

[[Imported from SVN: r7547]]

dune/gfe/Makefile.am

@@ -11,6 +11,7 @@ srcinclude_HEADERS = averagedistanceassembler.hh \
                      geodesicfeassembler.hh \
                      geodesicfefunctionadaptor.hh \
                      harmonicenergystiffness.hh \
+                     linearalgebra.hh \
                      localgeodesicfefunction.hh \
                      localgeodesicfestiffness.hh \
                      maxnormtrustregion.hh \

dune/gfe/linearalgebra.hh

+#ifndef LINEAR_ALGEBRA_HH
+#define LINEAR_ALGEBRA_HH
+
+#include <dune/common/fmatrix.hh>
+
+//! calculates ret = A * B
+template< class K, int m, int n, int p >
+Dune::FieldMatrix< K, m, p > operator* ( const Dune::FieldMatrix< K, m, n > &A, const Dune::FieldMatrix< K, n, p > &B)
+{
+    typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
+    Dune::FieldMatrix< K, m, p > ret;
+        
+    for( size_type i = 0; i < m; ++i ) {
+        
+        for( size_type j = 0; j < p; ++j ) {
+            ret[ i ][ j ] = K( 0 );
+            for( size_type k = 0; k < n; ++k )
+                ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
+        }
+    }
+    return ret;
+}
+
+//! calculates ret = A - B
+template< class K, int m, int n>
+Dune::FieldMatrix<K,m,n> operator- ( const Dune::FieldMatrix<K, m, n> &A, const Dune::FieldMatrix<K,m,n> &B)
+{
+    typedef typename Dune::FieldMatrix<K,m,n> :: size_type size_type;
+    Dune::FieldMatrix<K,m,n> ret;
+        
+    for( size_type i = 0; i < m; ++i )
+        for( size_type j = 0; j < n; ++j )
+            ret[i][j] = A[i][j] - B[i][j];
+
+    return ret;
+}
+
+
+#endif

dune/gfe/localgeodesicfefunction.hh

@@ -10,39 +10,7 @@
 #include <dune/gfe/targetspacertrsolver.hh>
 
 #include <dune/gfe/tensor3.hh>
-
-//! calculates ret = A * B
-template< class K, int m, int n, int p >
-Dune::FieldMatrix< K, m, p > operator* ( const Dune::FieldMatrix< K, m, n > &A, const Dune::FieldMatrix< K, n, p > &B)
-{
-    typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
-    Dune::FieldMatrix< K, m, p > ret;
-        
-    for( size_type i = 0; i < m; ++i ) {
-        
-        for( size_type j = 0; j < p; ++j ) {
-            ret[ i ][ j ] = K( 0 );
-            for( size_type k = 0; k < n; ++k )
-                ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
-        }
-    }
-    return ret;
-}
-
-//! calculates ret = A - B
-template< class K, int m, int n>
-Dune::FieldMatrix<K,m,n> operator- ( const Dune::FieldMatrix<K, m, n> &A, const Dune::FieldMatrix<K,m,n> &B)
-{
-    typedef typename Dune::FieldMatrix<K,m,n> :: size_type size_type;
-    Dune::FieldMatrix<K,m,n> ret;
-        
-    for( size_type i = 0; i < m; ++i )
-        for( size_type j = 0; j < n; ++j )
-            ret[i][j] = A[i][j] - B[i][j];
-
-    return ret;
-}
-
+#include <dune/gfe/linearalgebra.hh>
 
 /** \brief A function defined by simplicial geodesic interpolation 
            from the reference element to a Riemannian manifold.