From 26d119b201a066242b5851455562765ec889fb94 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 12 Dec 2013 19:11:51 +0000
Subject: [PATCH] Use the Gram-Schmidt solver instead of CG
CG seems to converge to quickly. If too few iterations are performed,
then ADOL-C has trouble determining the correct second derivatives.
Since the Gram-Schmidt solver is a direct solver, OTOH, the derivatives
should be perfect. It is also faster than the Gauss-Seidel solver I
was previously using.
[[Imported from SVN: r9581]]
---
dune/gfe/targetspacertrsolver.cc | 53 +++++++++++++-------------------
dune/gfe/targetspacertrsolver.hh | 17 +++++-----
2 files changed, 30 insertions(+), 40 deletions(-)
diff --git a/dune/gfe/targetspacertrsolver.cc b/dune/gfe/targetspacertrsolver.cc
index a3b91483..25909b44 100644
--- a/dune/gfe/targetspacertrsolver.cc
+++ b/dune/gfe/targetspacertrsolver.cc
@@ -6,6 +6,8 @@
#include "maxnormtrustregion.hh"
+#include <dune/gfe/gramschmidtsolver.hh>
+
template <class TargetSpace>
void TargetSpaceRiemannianTRSolver<TargetSpace>::
setup(const AverageDistanceAssembler<TargetSpace>* assembler,
@@ -26,21 +28,7 @@ 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
+#ifdef USE_GAUSS_SEIDEL_SOLVER
// ////////////////////////////////
// Create a projected gauss-seidel solver
// ////////////////////////////////
@@ -86,26 +74,27 @@ 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
+#ifdef USE_GAUSS_SEIDEL_SOLVER
corr = 0;
#endif
MatrixType hesseMatrix(1,1);
+#ifdef USE_GAUSS_SEIDEL_SOLVER
assembler_->assembleGradient(x_, rhs[0]);
assembler_->assembleHessian(x_, hesseMatrix[0][0]);
+#else
+ /** \todo Fix this sense copying */
+ typename TargetSpace::EmbeddedTangentVector foo;
+ assembler_->assembleEmbeddedGradient(x_, foo);
+ rhs[0] = foo;
+ assembler_->assembleEmbeddedHessian(x_, hesseMatrix[0][0]);
+#endif
// 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
+#ifdef USE_GAUSS_SEIDEL_SOLVER
dynamic_cast<LinearIterationStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->setProblem(hesseMatrix, corr, rhs);
dynamic_cast<TrustRegionGSStep<MatrixType,CorrectionType>*>(innerSolver_->iterationStep_)->obstacles_ = &trustRegion.obstacles();
@@ -115,14 +104,17 @@ void TargetSpaceRiemannianTRSolver<TargetSpace>::solve()
// /////////////////////////////
// Solve !
// /////////////////////////////
-
+#ifdef USE_GAUSS_SEIDEL_SOLVER
innerSolver_->solve();
-
-#ifdef USE_TCGSOLVER
- corr = *innerSolver_->x_;
#else
+ Dune::FieldMatrix<field_type,blocksize,embeddedBlocksize> basis = x_.orthonormalFrame();
+ GramSchmidtSolver<field_type, blocksize, embeddedBlocksize>::solve(hesseMatrix[0][0], corr[0], rhs[0], basis);
+#endif
+
+#ifdef USE_GAUSS_SEIDEL_SOLVER
corr = innerSolver_->iterationStep_->getSol();
#endif
+
//std::cout << "Corr: " << corr << std::endl;
if (this->verbosity_ == NumProc::FULL)
@@ -153,11 +145,8 @@ 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;
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index 99ba1646..9ffcc002 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -2,14 +2,12 @@
#define TARGET_SPACE_RIEMANNIAN_TRUST_REGION_SOLVER_HH
-//#define USE_TCGSOLVER
+//#define USE_GAUSS_SEIDEL_SOLVER
#include <dune/istl/matrix.hh>
#include <dune/solvers/common/boxconstraint.hh>
-#ifdef USE_TCGSOLVER
-#include <dune/solvers/solvers/tcgsolver.hh>
-#else
+#ifdef USE_GAUSS_SEIDEL_SOLVER
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
#endif
@@ -23,6 +21,7 @@ class TargetSpaceRiemannianTRSolver
: public NumProc
{
const static int blocksize = TargetSpace::TangentVector::dimension;
+ const static int embeddedBlocksize = TargetSpace::EmbeddedTangentVector::dimension;
// Centralize the field type here
typedef typename TargetSpace::ctype field_type;
@@ -30,8 +29,13 @@ class TargetSpaceRiemannianTRSolver
// Some types that I need
// The types have the dynamic outer type because the dune-solvers solvers expect
// this sort of type.
+#ifdef USE_GAUSS_SEIDEL_SOLVER
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::BlockVector<Dune::FieldVector<field_type, embeddedBlocksize> > CorrectionType;
+#endif
public:
@@ -79,10 +83,7 @@ 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
+#ifdef USE_GAUSS_SEIDEL_SOLVER
/** \brief The solver for the quadratic inner problems */
std::auto_ptr< ::LoopSolver<CorrectionType> > innerSolver_;
--
GitLab