diff --git a/dune/solvers/solvers/CMakeLists.txt b/dune/solvers/solvers/CMakeLists.txt
index a2e7d1c6b629ea16fb66094512abd83940f8d40f..119244ed6efad8b2bd6d43c9a9d3d937fbf178f9 100644
--- a/dune/solvers/solvers/CMakeLists.txt
+++ b/dune/solvers/solvers/CMakeLists.txt
@@ -9,4 +9,6 @@ install(FILES
solver.hh
tcgsolver.cc
tcgsolver.hh
+ trustregionsolver.cc
+ trustregionsolver.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers/solvers)
diff --git a/dune/solvers/solvers/Makefile.am b/dune/solvers/solvers/Makefile.am
index fd48820fe5cefad57388951cd1f5e1f2802c88e9..c7d0879fc36b0d41eae4527735482eff033b343a 100644
--- a/dune/solvers/solvers/Makefile.am
+++ b/dune/solvers/solvers/Makefile.am
@@ -3,7 +3,7 @@ SUBDIRS =
solversdir = $(includedir)/dune/solvers/solvers
solvers_HEADERS = cgsolver.cc cgsolver.hh iterativesolver.cc iterativesolver.hh \
loopsolver.cc loopsolver.hh quadraticipopt.hh solver.hh \
- tcgsolver.cc tcgsolver.hh
+ tcgsolver.cc tcgsolver.hh trustregionsolver.cc trustregionsolver.hh
EXTRA_DIST = CMakeLists.txt
diff --git a/dune/solvers/solvers/trustregionsolver.cc b/dune/solvers/solvers/trustregionsolver.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2feb30593812e6a80449af8c02e3d44f9819dae3
--- /dev/null
+++ b/dune/solvers/solvers/trustregionsolver.cc
@@ -0,0 +1,179 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/io.hh>
+
+// Using a monotone multigrid as the inner solver
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+template <class ProblemType, class VectorType, class MatrixType>
+void TrustRegionSolver<ProblemType,VectorType,MatrixType>::solve()
+{
+ MonotoneMGStep<MatrixType,VectorType>* mgStep = nullptr;
+
+ // if the inner solver is a monotone multigrid set up a max-norm trust-region
+ if (dynamic_cast<::LoopSolver<VectorType>*>(innerSolver_.get())) {
+ mgStep = dynamic_cast<MonotoneMGStep<MatrixType,VectorType>*>(dynamic_cast<::LoopSolver<VectorType>*>(innerSolver_.get())->iterationStep_);
+
+ }
+
+ // setup the trust-region in the maximums norm
+ std::vector<BoxConstraint<field_type,blocksize> > trustRegionObstacles(x_.size());
+ for (size_t i=0; i<trustRegionObstacles.size(); i++)
+ for (int j=0; j<blocksize; j++) {
+ trustRegionObstacles[i].lower(j) = -initialRadius_;
+ trustRegionObstacles[i].upper(j) = initialRadius_;
+ }
+
+ mgStep->obstacles_ = &trustRegionObstacles;
+
+ // /////////////////////////////////////////////////////
+ // Trust-Region Solver
+ // /////////////////////////////////////////////////////
+
+ double oldEnergy = problem_->energy(x_);
+
+ for (int i=0; i<this->maxIterations_; i++) {
+
+ Dune::Timer totalTimer;
+ if (this->verbosity_ == Solver::FULL) {
+ std::cout << "----------------------------------------------------" << std::endl;
+ std::cout << " Trust-Region Step Number: " << i
+ << ", radius: " << initialRadius_
+ << ", energy: " << oldEnergy << std::endl;
+ std::cout << "----------------------------------------------------" << std::endl;
+ }
+
+ // Assemble the quadratic problem
+ Dune::Timer gradientTimer;
+ problem_->assembleQP(x_);
+
+ // Compute gradient norm to monitor convergence
+ VectorType gradient = problem_->f_;
+ for (size_t j=0; j<gradient.size(); j++)
+ for (int k=0; k<gradient[j].size(); k++)
+ if ((*mgStep->ignoreNodes_)[j][k])
+ gradient[j][k] = 0;
+
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Gradient norm: " << this->errorNorm_->operator()(gradient) << std::endl;
+
+ if (this->verbosity_ == Solver::FULL)
+ std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
+ VectorType corr(x_.size());
+ corr = 0;
+
+ mgStep->setProblem(problem_->A_, corr, problem_->f_);
+
+ innerSolver_->preprocess();
+
+ ///////////////////////////////
+ // Solve !
+ ///////////////////////////////
+
+ // Solve the subproblem until we find an acceptable iterate
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Solve quadratic problem..." << std::endl;
+ while(true) {
+
+ Dune::Timer solutionTimer;
+ innerSolver_->solve();
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+
+ corr = mgStep->getSol();
+
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+ if (corr.infinity_norm() < this->tolerance_) {
+ if (this->verbosity_ != NumProc::QUIET)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ return;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether trust-region step can be accepted
+ // ////////////////////////////////////////////////////
+
+ VectorType newIterate = x_;
+ newIterate += corr;
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of new iterate : " << newIterate.infinity_norm() << std::endl;
+
+ field_type energy = problem_->energy(newIterate);
+
+ // compute the model decrease
+ field_type modelDecrease = problem_->modelDecrease(corr);
+
+ field_type relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+ if (this->verbosity_ == NumProc::FULL) {
+ std::cout << "Absolute model decrease: " << modelDecrease
+ << ", functional decrease: " << oldEnergy - energy << std::endl;
+ std::cout << "Relative model decrease: " << relativeModelDecrease
+ << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+ }
+
+ // assure that the inner solver works correctly
+ assert(modelDecrease >= 0);
+
+ if (energy >= oldEnergy) {
+ if (this->verbosity_ == NumProc::FULL)
+ printf("No energy descent!\n");
+ }
+
+ if (energy >= oldEnergy &&
+ (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Suspecting rounding problems, stopping here!" << std::endl;
+
+ if (this->verbosity_ != NumProc::QUIET)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+ x_ = newIterate;
+ return;
+ }
+
+ // //////////////////////////////////////////////
+ // Check for acceptance of the step
+ // //////////////////////////////////////////////
+ if ( (oldEnergy-energy) / modelDecrease > 0.9) {
+ // very successful iteration
+
+ x_ = newIterate;
+ scaleTrustRegionRadius(trustRegionObstacles,enlargeRadius_);
+ std::cout << "Very Successful iteration, enlarging radius! " <<initialRadius_<<std::endl;
+
+ // current energy becomes 'oldEnergy' for the next iteration
+ oldEnergy = energy;
+ break;
+
+ } else if ( (oldEnergy-energy) / modelDecrease > 0.01
+ || std::abs(oldEnergy-energy) < 1e-12) {
+ // successful iteration
+ x_ = newIterate;
+
+ // current energy becomes 'oldEnergy' for the next iteration
+ oldEnergy = energy;
+ break;
+
+ } else {
+
+ // unsuccessful iteration
+
+ // Decrease the trust-region radius
+ scaleTrustRegionRadius(trustRegionObstacles,shrinkRadius_);
+
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Unsuccessful iteration, shrinking radius! " <<initialRadius_<<std::endl;
+ }
+ }
+ std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
+ }
+
+}
diff --git a/dune/solvers/solvers/trustregionsolver.hh b/dune/solvers/solvers/trustregionsolver.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9438b23ea2ad07aae79bfbdbbece0cd393817345
--- /dev/null
+++ b/dune/solvers/solvers/trustregionsolver.hh
@@ -0,0 +1,147 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_TRUST_REGION_SOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_TRUST_REGION_SOLVER_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+/** \brief A standard Trust-region solver with infinity-norm for the trust-region constraints, to solve an unconstrained nonlinear (possibly nonconvex) minimization problem.
+ *
+ * The energy functional is provided by the ProblemType class. It has to provide the following methods:
+ *
+ * void assembleQP(const VectorType& iterate)
+ * field_type energy(const VectorType& iterate)
+ * field_type modelDecrease(const VectorType& iterate)
+ *
+ * and member:
+ * MatrixType A_
+ * VectorType f_
+ *
+ * The method 'assembleQP' should assemble the local quadratic approximation at 'iterate'
+ * The local model is to be given by the quadratic part A_ and the linear part f_
+ *
+ * The method 'energy' should evaluate the nonlinear energy at 'iterate'.
+ *
+ * The method 'modelDecrease' should compute the decrease in the model-energy.
+ *
+ * The inner solver that has to be provided should be able to solve non-convex box-constrained quadratic problems,
+ * if the nonlinear functional is non-convex itself.
+ *
+ * \tparam ProblemType A class representing the energy functional
+ *
+ * */
+template <class ProblemType, class VectorType, class MatrixType>
+class TrustRegionSolver
+ : public IterativeSolver<VectorType>
+{
+ const static int blocksize = VectorType::value_type::dimension;
+ typedef IterativeSolver<VectorType> Base;
+ typedef typename VectorType::field_type field_type;
+
+public:
+
+ TrustRegionSolver()
+ : Base(0,100,NumProc::FULL), enlargeRadius_(2), shrinkRadius_(0.5),
+ goodIteration_(0.9), badIteration_(0.01)
+ {}
+
+ TrustRegionSolver(field_type tolerance, int maxIterations, NumProc::VerbosityMode verbosity=NumProc::FULL)
+ : Base(tolerance, maxIterations, verbosity), enlargeRadius_(2), shrinkRadius_(0.5),
+ goodIteration_(0.9), badIteration_(0.01)
+ {}
+
+ TrustRegionSolver(const Dune::ParameterTree& trConfig)
+ : Base(0,100,NumProc::FULL)
+ {
+ setupTrParameter(trConfig);
+ }
+
+ //! Setup the standard trust-region parameter
+ void setupTrParameter(const Dune::ParameterTree& trConfig) {
+
+ this->tolerance_ = trConfig.get<field_type>("tolerance");
+ this->maxIterations_ = trConfig.get<int>("maxIterations");
+ this->verbosity_ = trConfig.get("verbosity",NumProc::FULL);
+ initialRadius_ = trConfig.get<field_type>("initialRadius");
+ maxRadius_ = trConfig.get<field_type>("maxRadius");
+ enlargeRadius_ = trConfig.get("enlargeRadius",2.0);
+ shrinkRadius_ = trConfig.get("shrinkRadius",0.5);
+ goodIteration_ = trConfig.get("goodIteration",0.9);
+ badIteration_ = trConfig.get("badIteration",0.01);
+
+ }
+
+
+ /** \brief Set up the solver. */
+ void setup(const Dune::ParameterTree& trConfig,
+ Norm<VectorType>& errorNorm,
+ ProblemType& problem, Solver& innerSolver) {
+
+ setupTrParameter(trConfig);
+ this->errorNorm_ = &errorNorm;
+ problem_ = Dune::stackobject_to_shared_ptr(problem);
+ innerSolver_ = Dune::stackobject_to_shared_ptr(innerSolver);
+ }
+
+ //! Solve the problem.
+ void solve();
+
+ //! Set the initial iterate
+ void setInitialIterate(const VectorType& x) {x_ = x;}
+
+ VectorType getSol() const {return x_;}
+
+protected:
+
+ //! Adjust the trust-region obstacles
+ void scaleTrustRegionRadius(std::vector<BoxConstraint<field_type,blocksize> >& obs, const field_type& scaling)
+ {
+ if (scaling*initialRadius_>maxRadius_)
+ initialRadius_ = maxRadius_;
+ else
+ initialRadius_ *= scaling;
+
+ for (size_t i=0; i<obs.size(); i++)
+ for (int j=0; j<blocksize; j++) {
+ obs[i].lower(j) = -initialRadius_;
+ obs[i].upper(j) = initialRadius_;
+ }
+ }
+
+ /** \brief The problem type. */
+ std::shared_ptr<ProblemType> problem_;
+
+ /** \brief The solution vector */
+ VectorType x_;
+
+ /** \brief The solver for the quadratic inner problems */
+ std::shared_ptr<Solver> innerSolver_;
+
+ /** \brief The initial trust-region radius in the maximum-norm */
+ field_type initialRadius_;
+
+ /** \brief The maximal trust-region radius in the maximum-norm */
+ field_type maxRadius_;
+
+ /** \brief If the iteration is very successfull we enlarge the radius by this factor.*/
+ field_type enlargeRadius_;
+
+ /** \brief If the iteration is unsuccessfull we shrink the radius by this factor.*/
+ field_type shrinkRadius_;
+
+ /** \brief If the ratio between predicted and achieved decrease is above this number, the iteration is very successful.*/
+ field_type goodIteration_;
+
+ /** \brief If the ratio between predicted and achieved decrease is lower than this number, the iteration is unsuccessful.*/
+ field_type badIteration_;
+};
+
+#include "trustregionsolver.cc"
+
+#endif