diff --git a/dune/gfe/Makefile.am b/dune/gfe/Makefile.am
index 448e15dd44aecf48ce909db430d7f9a89ecd5d42..0bc5d92341f833a555cceef1d500b702d892befe 100644
--- a/dune/gfe/Makefile.am
+++ b/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 \
diff --git a/dune/gfe/linearalgebra.hh b/dune/gfe/linearalgebra.hh
new file mode 100644
index 0000000000000000000000000000000000000000..04f464fa76839b3c336c855d5826c05d77909c46
--- /dev/null
+++ b/dune/gfe/linearalgebra.hh
@@ -0,0 +1,39 @@
+#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
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 74765c7af9b972d2a31fc5a3f5bf08ff66498945..eeb2f42c6ead9bf0ccc83277c1f3cf0932808be3 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/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.