Really use a SymmetricMatrix for the Hesse matrix

Assuming that the AverageDistanceAssembler can provide it.

This is for the Gram-Schmidt solver only.  The old Gauss-Seidel
code will still use a full matrix.

[[Imported from SVN: r9585]]

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

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