From b9ec02a9a6e1de54b03b76410c979f541698741e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 12 Dec 2013 21:35:28 +0000
Subject: [PATCH] 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 | 12 +++++-------
dune/gfe/targetspacertrsolver.hh | 8 +++++---
2 files changed, 10 insertions(+), 10 deletions(-)
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index ca3d9e05..2e9a2d3b 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 9ffcc002..10e2d87f 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
--
GitLab