diff --git a/dune/gfe/gramschmidtsolver.hh b/dune/gfe/gramschmidtsolver.hh
index 7e9be4d676f6a1ea9de20b47d5f14dd95393c838..3623ff6b067607e5eb0d657f89ab272d5eb7ff78 100644
--- a/dune/gfe/gramschmidtsolver.hh
+++ b/dune/gfe/gramschmidtsolver.hh
@@ -2,8 +2,8 @@
 #define DUNE_GFE_GRAMSCHMIDTSOLVER_HH
 
 #include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
 
+#include <dune/gfe/symmetricmatrix.hh>
 
 /** \brief Direct solver for a dense symmetric linear system, using an orthonormal basis
  *
@@ -24,16 +24,10 @@ class GramSchmidtSolver
    * \param matrix The matrix inducing the matrix norm
    * \param[in,out] v The vector to normalize
    */
-  static void normalize(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+  static void normalize(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
                         Dune::FieldVector<field_type,embeddedDim>& v)
   {
-    field_type energyNormSquared = 0;
-
-    for (int i=0; i<v.size(); i++)
-      for (int j=0; j<v.size(); j++)
-        energyNormSquared += matrix[i][j]*v[i]*v[j];
-
-    v /= std::sqrt(energyNormSquared);
+    v /= std::sqrt(matrix.energyScalarProduct(v,v));
   }
 
 
@@ -41,16 +35,12 @@ class GramSchmidtSolver
    *
    * \param matrix The matrix the defines the scalar product
    */
-  static void project(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+  static void project(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
                       const Dune::FieldVector<field_type,embeddedDim>& vi,
                       Dune::FieldVector<field_type,embeddedDim>& vj)
   {
 
-    field_type energyScalarProduct = 0;
-
-    for (int i=0; i<vi.size(); i++)
-      for (int j=0; j<vj.size(); j++)
-        energyScalarProduct += matrix[i][j]*vi[i]*vj[j];
+    field_type energyScalarProduct = matrix.energyScalarProduct(vi,vj);
 
     for (int i=0; i<vj.size(); i++)
       vj[i] -= energyScalarProduct * vi[i];
@@ -59,7 +49,7 @@ class GramSchmidtSolver
 
 public:
   /** Solve linear system by constructing an energy-orthonormal basis */
-  static void solve(const Dune::FieldMatrix<field_type,embeddedDim,embeddedDim>& matrix,
+  static void solve(const Dune::SymmetricMatrix<field_type,embeddedDim>& matrix,
                     Dune::FieldVector<field_type,embeddedDim>& x,
                     const Dune::FieldVector<field_type,embeddedDim>& rhs,
                     const Dune::FieldMatrix<field_type,dim,embeddedDim>& basis)
diff --git a/dune/gfe/symmetricmatrix.hh b/dune/gfe/symmetricmatrix.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c46112315289bdef4e316d56ac074d3b8d4d5a27
--- /dev/null
+++ b/dune/gfe/symmetricmatrix.hh
@@ -0,0 +1,68 @@
+#ifndef DUNE_GFE_SYMMETRICMATRIX_HH
+#define DUNE_GFE_SYMMETRICMATRIX_HH
+
+#include <dune/common/fvector.hh>
+
+namespace Dune {
+
+/** \brief  A class implementing a symmetric matrix with compile-time size
+ *
+ *  A \f$ dim\times dim \f$ matrix is stored internally as a <tt>Dune::FieldVector<double, dim*(dim+1)/2></tt>
+ *  The components are assumed to be \f$ [ E(1,1),\  E(2,1),\ E(2,2),\ E(3,1),\ E(3,2),\ E(3,3) ]\f$
+ *  and analogous for other dimensions
+ *  \tparam  dim number of lines/columns of the matrix
+ */
+template <class T, int N>
+class SymmetricMatrix
+{
+
+public:
+    /** \brief Default constructor
+     *
+     *  Tensor is initialized containing zeros if no argument is given.
+     *  \param eye if true tensor is initialized as identity
+     */
+    SymmetricMatrix()
+    {}
+
+    /** \brief Matrix style random read/write access to components 
+     *  \param i line index
+     *  \param j column index
+     * \note You need to know what you are doing:  You can only access the lower
+     * left triangular submatrix using this methods.  It requires i<=j.
+     */
+    T& operator() (int i, int j)
+    {
+      assert(i>=j);
+      return data_[i*(i+1)/2+j];
+    }
+
+    /** \brief Matrix style random read access to components
+     *  \param i line index
+     *  \param j column index
+     * \note You need to know what you are doing:  You can only access the lower
+     * left triangular submatrix using this methods.  It requires i<=j.
+     */
+    const T& operator() (int i, int j) const
+    {
+      assert(i>=j);
+      return data_[i*(i+1)/2+j];
+    }
+
+    T energyScalarProduct (const FieldVector<T,N>& v1, const FieldVector<T,N>& v2) const
+    {
+      T result = 0;
+      for (size_t i=0; i<N; i++)
+        for (size_t j=0; j<=i; j++)
+          result += (1+(i!=j)) * operator()(i,j) * v1[i] * v2[j];
+
+      return result;
+    }
+
+
+private:
+    Dune::FieldVector<T,N*(N+1)/2> data_;
+};
+
+}
+#endif
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index 25909b441d8dd8522108cf79d2689f4aa32a66ca..ca3d9e05a60415ab5d87be1ac65cbf858671fb91 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -107,8 +107,14 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
 #ifdef USE_GAUSS_SEIDEL_SOLVER
         innerSolver_->solve();
 #else
+        Dune::SymmetricMatrix<field_type, embeddedBlocksize> symmetricHessian;
+        for (size_t j=0; j<embeddedBlocksize; j++)
+          for (size_t k=0; k<=j; k++)
+            symmetricHessian(j,k) = hesseMatrix[0][0][j][k];
+
+
         Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> basis = x_.orthonormalFrame();
-        GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix[0][0], corr[0], rhs[0], basis);
+        GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(symmetricHessian, corr[0], rhs[0], basis);
 #endif
 
 #ifdef USE_GAUSS_SEIDEL_SOLVER