diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index ca3d9e05a60415ab5d87be1ac65cbf858671fb91..2e9a2d3bb80a064c3510ac41d3c38bf976c07d43 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -107,14 +107,8 @@ 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(symmetricHessian, corr[0], rhs[0], basis);
+        GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix[0][0], corr[0], rhs[0], basis);
 #endif
 
 #ifdef USE_GAUSS_SEIDEL_SOLVER
@@ -148,10 +142,14 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
         // compute the model decrease
         // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
         // Note that rhs = -g
+#ifdef USE_GAUSS_SEIDEL_SOLVER
         CorrectionType tmp(corr.size());
         tmp = 0;
         hesseMatrix.umv(corr, tmp);
         field_type modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+#else
+        field_type modelDecrease = (rhs*corr) - 0.5 * hesseMatrix[0][0].energyScalarProduct(corr[0],corr[0]);
+#endif
 
         if (this->verbosity_ == NumProc::FULL) {
             std::cout << "Absolute model decrease: " << modelDecrease
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index 9ffcc002576beffdd7c75adc487e9afb26308046..10e2d87f60f9b6e29bc7dca3b8d05b72156ed133 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -13,6 +13,8 @@
 #endif
 #include <dune/solvers/norms/energynorm.hh>
 
+#include <dune/gfe/symmetricmatrix.hh>
+
 /** \brief Riemannian trust-region solver for geodesic finite-element problems
  \tparam TargetSpace The manifold that our functions take values in
  */
@@ -33,7 +35,7 @@ class TargetSpaceRiemannianTRSolver
     typedef Dune::Matrix<Dune::FieldMatrix<field_type, blocksize, blocksize> > MatrixType;
     typedef Dune::BlockVector<Dune::FieldVector<field_type, blocksize> >       CorrectionType;
 #else
-    typedef Dune::Matrix<Dune::FieldMatrix<field_type, embeddedBlocksize, embeddedBlocksize> > MatrixType;
+    typedef Dune::Matrix<Dune::SymmetricMatrix<field_type, embeddedBlocksize> > MatrixType;
     typedef Dune::BlockVector<Dune::FieldVector<field_type, embeddedBlocksize> >       CorrectionType;
 #endif
 
@@ -89,10 +91,10 @@ protected:
 
     /** \brief The iteration step for the quadratic inner problems */
     std::auto_ptr<TrustRegionGSStep<MatrixType, CorrectionType> > innerSolverStep_;
-#endif
+
     /** \brief Norm for the quadratic inner problems */
     std::auto_ptr<EnergyNorm<MatrixType, CorrectionType> > energyNorm_;
-
+#endif
     /** \brief Specify a minimal number of iterations the trust-region solver has to do
      *
      * This is needed when working with automatic differentiation.    While a very low