From f8f468aa810eea67d00f73e7ec04d518a9174b11 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 10 Dec 2013 08:53:16 +0000
Subject: [PATCH] Alternatively allow a CG as the quadratic solver for the
TargetSpace-RTR
I hoped that this would be faster. However, it seems to deterioate the
convergence of the outer trust region algorithm. Apparently, the steps
taken by a CG make it difficult for ADOL-C to compute the correct
derivatives.
[[Imported from SVN: r9576]]
---
dune/gfe/targetspacertrsolver.cc | 41 ++++++++++++++++++++++++++++----
dune/gfe/targetspacertrsolver.hh | 13 +++++++++-
2 files changed, 49 insertions(+), 5 deletions(-)
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index b197b73d..a3b91483 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -26,6 +26,21 @@ setup(const AverageDistanceAssembler<TargetSpace>* assembler,
this->verbosity_ = NumProc::QUIET;
minNumberOfIterations_ = 4;
+#ifdef USE_TCGSOLVER
+ ///////////////////////////////////////////////////
+ // Create a truncated CG solver
+ ///////////////////////////////////////////////////
+
+ innerSolver_ = std::auto_ptr<TruncatedCGSolver<MatrixType, CorrectionType> >
+ (new TruncatedCGSolver<MatrixType, CorrectionType>(NULL,
+ 3,//innerIterations,
+ -1,//innerTolerance,
+ NULL, //energyNorm_.get(),
+ NULL,
+ Solver::QUIET,
+ false));
+
+#else
// ////////////////////////////////
// Create a projected gauss-seidel solver
// ////////////////////////////////
@@ -42,7 +57,7 @@ setup(const AverageDistanceAssembler<TargetSpace>* assembler,
Solver::QUIET));
innerSolver_->useRelativeError_ = false;
-
+#endif
}
@@ -71,7 +86,13 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
CorrectionType rhs(1); // length is 1 _block_
CorrectionType corr(1); // length is 1 _block_
+#ifdef USE_TCGSOLVER
+ // CG stops to early when started from zero.
+ // ADOLC needs at least a few iterations to pick up the correct derivatives
+ corr = 1;
+#else
corr = 0;
+#endif
MatrixType hesseMatrix(1,1);
@@ -81,19 +102,28 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
// The right hand side is the _negative_ gradient
rhs *= -1;
+#ifdef USE_TCGSOLVER
+ CorrectionType rhsBackup = rhs; // TruncatedCGSolver overwrites rhs
+ innerSolver_->setProblem(hesseMatrix, &corr, &rhs, trustRegion.obstacles()[0].upper(0));
+#else
dynamic_cast<LinearIterationStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->setProblem(hesseMatrix, corr, rhs);
dynamic_cast<TrustRegionGSStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->obstacles_ = &trustRegion.obstacles();
innerSolver_->preprocess();
-
+#endif
// /////////////////////////////
// Solve !
// /////////////////////////////
innerSolver_->solve();
+#ifdef USE_TCGSOLVER
+ corr = *innerSolver_->x_;
+#else
corr = innerSolver_->iterationStep_->getSol();
+#endif
+ //std::cout << "Corr: " << corr << std::endl;
if (this->verbosity_ == NumProc::FULL)
std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
@@ -123,8 +153,11 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
CorrectionType tmp(corr.size());
tmp = 0;
hesseMatrix.umv(corr, tmp);
+#ifdef USE_TCGSOLVER
+ field_type modelDecrease = (rhsBackup*corr) - 0.5 * (corr*tmp);
+#else
field_type modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-
+#endif
if (this->verbosity_ == NumProc::FULL) {
std::cout << "Absolute model decrease: " << modelDecrease
<< ", functional decrease: " << oldEnergy - energy << std::endl;
@@ -132,7 +165,7 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
<< ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
}
- assert(modelDecrease >= 0);
+ assert(modelDecrease >= -1e-15);
if (energy >= oldEnergy) {
if (this->verbosity_ == NumProc::FULL)
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index bd3baea6..99ba1646 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -1,11 +1,18 @@
#ifndef TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
#define TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
+
+//#define USE_TCGSOLVER
+
#include <dune/istl/matrix.hh>
#include <dune/solvers/common/boxconstraint.hh>
+#ifdef USE_TCGSOLVER
+#include <dune/solvers/solvers/tcgsolver.hh>
+#else
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
+#endif
#include <dune/solvers/norms/energynorm.hh>
/** \brief Riemannian trust-region solver for geodesic finite-element problems
@@ -72,12 +79,16 @@ protected:
/** \brief The assembler for the average-distance functional */
const AverageDistanceAssembler<TargetSpace>* assembler_;
+#ifdef USE_TCGSOLVER
+ /** \brief The solver for the quadratic inner problems */
+ std::auto_ptr<TruncatedCGSolver<MatrixType, CorrectionType> > innerSolver_;
+#else
/** \brief The solver for the quadratic inner problems */
std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;
/** \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_;
--
GitLab