dune/solvers/norms/sumnorm.hh

-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef SUMNORM_HH
 #define SUMNORM_HH
 
@@ -12,7 +12,10 @@ class SumNorm: public Norm<V>
 {
     public:
         typedef V VectorType;
-        typedef typename VectorType::field_type field_type;
+        using Base = Norm<V>;
+
+        /** \brief The type used for the result */
+        using typename Base::field_type;
 
         SumNorm(field_type _alpha1, const Norm<VectorType> &_norm1, field_type _alpha2, const Norm<VectorType> &_norm2) :
             alpha1(_alpha1),
@@ -22,13 +25,13 @@ class SumNorm: public Norm<V>
     {}
 
         //! Compute the norm of the given vector
-        field_type operator()(const VectorType &f) const
+        field_type operator()(const VectorType &f) const override
         {
             return std::sqrt(normSquared(f));
         }
 
         //! Compute the square of the norm of the given vector
-        virtual field_type normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
             field_type r1 = norm1.normSquared(f);
             field_type r2 = norm2.normSquared(f);
@@ -37,7 +40,7 @@ class SumNorm: public Norm<V>
         }
 
         //! Compute the norm of the difference of two vectors
-        field_type diff(const VectorType &f1, const VectorType &f2) const
+        field_type diff(const VectorType &f1, const VectorType &f2) const override
         {
             field_type r1 = norm1.diff(f1,f2);
             field_type r2 = norm2.diff(f1,f2);

dune/solvers/norms/twonorm.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef TWONORM_HH
 #define TWONORM_HH
 
@@ -5,43 +7,52 @@
 #include <cstring> // For size_t
 
 #include <dune/common/fvector.hh>
+#include <dune/common/hybridutilities.hh>
 
 #include "norm.hh"
 
 //! Wrapper around the two_norm() method of a vector class
-template <class VectorType>
-class TwoNorm : public Norm<VectorType>
+template <class V>
+class TwoNorm : public Norm<V>
 {
-
     public:
+       typedef V VectorType;
+       using Base = Norm<V>;
+
+       /** \brief The type used for the result */
+       using typename Base::field_type;
 
         /** \brief Destructor, doing nothing */
-        virtual ~TwoNorm() {};
+        virtual ~TwoNorm() {}
 
         //! Compute the norm of the given vector
-        virtual double operator()(const VectorType& f) const
+        virtual field_type operator()(const VectorType& f) const override
         {
             return f.two_norm();
         }
 
         //! Compute the square of the norm of the given vector
-        virtual double normSquared(const VectorType& f) const
+        virtual field_type normSquared(const VectorType& f) const override
         {
             return f.two_norm2();
         }
 
+        // TODO there is no reason to have this if f.two_norm2() can be
+        //      assumed (above) anyway. Then also the specialization for
+        //      FieldVector (below) can go.
         //! Compute the norm of the difference of two vectors
-        virtual double diff(const VectorType& f1, const VectorType& f2) const
+        virtual field_type diff(const VectorType& f1, const VectorType& f2) const override
         {
-            double r = 0.0;
+            assert(f1.size() == f2.size());
+            field_type r = 0.0;
 
-            for (size_t i=0; i<f1.size(); ++i)
+            Dune::Hybrid::forEach(Dune::Hybrid::integralRange(Dune::Hybrid::size(f1)), [&](auto&& i)
             {
-                typename VectorType::block_type block = f1[i];
+                auto block = f1[i];
                 block -= f2[i];
 
-                r += block.two_norm2();
-            }
+                r += Dune::Impl::asVector(block).two_norm2();
+            });
 
             return std::sqrt(r);
         }
@@ -60,19 +71,19 @@ class TwoNorm<Dune::FieldVector<T, dimension> >
         virtual ~TwoNorm() {};
 
         //! Compute the norm of the given vector
-        virtual double operator()(const VectorType& f) const
+        virtual T operator()(const VectorType& f) const override
         {
             return f.two_norm();
         }
 
         //! Compute the square of the norm of the given vector
-        virtual double normSquared(const VectorType& f) const
+        virtual T normSquared(const VectorType& f) const override
         {
             return f.two_norm2();
         }
 
         //! Compute the norm of the difference of two vectors
-        virtual double diff(const VectorType& f1, const VectorType& f2) const
+        virtual T diff(const VectorType& f1, const VectorType& f2) const override
         {
           VectorType tmp = f1;
           tmp -= f2;

dune/solvers/operators/Makefile.am

-SUBDIRS =
-
-operatorsdir = $(includedir)/dune/solvers/operators
-operators_HEADERS = lowrankoperator.hh \
-                    nulloperator.hh \
-                    sumoperator.hh
-
-EXTRA_DIST = CMakeLists.txt
-
-include $(top_srcdir)/am/global-rules

dune/solvers/operators/lowrankoperator.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef LOWRANKOPERATOR_HH
 #define LOWRANKOPERATOR_HH
 

dune/solvers/operators/nulloperator.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef NULLOPERATOR_HH
 #define NULLOPERATOR_HH
 
@@ -22,11 +24,11 @@ class NullOperator
             public:
                 RowDummy(): zero_(0){}
 
-                const BlockType& operator[](size_t i) const
+                const BlockType& operator[]([[maybe_unused]] size_t i) const
                 {
                     return zero_;
                 }
-                BlockType& operator[](size_t i)
+                BlockType& operator[]([[maybe_unused]] size_t i)
                 {
                     return zero_;
                 }
@@ -58,7 +60,7 @@ class NullOperator
          *  Implements b += Nx and hence does nothing (N=0 !)
          */
         template <class LVectorType, class RVectorType>
-        void umv(const LVectorType& x, RVectorType& b) const
+        void umv([[maybe_unused]] const LVectorType& x, [[maybe_unused]] RVectorType& b) const
         {}
 
         /** \brief transposed Matrix-Vector multiplication
@@ -66,7 +68,7 @@ class NullOperator
          *  Implements b += N^tx and hence does nothing (N=0 !)
          */
         template <class LVectorType, class RVectorType>
-        void umtv(const LVectorType& x, RVectorType& b) const
+        void umtv([[maybe_unused]] const LVectorType& x, [[maybe_unused]] RVectorType& b) const
         {}
 
         /** \brief Matrix-Vector multiplication with scalar multiplication
@@ -74,7 +76,7 @@ class NullOperator
          *  Implements b += a*Nx and hence does nothing (N=0 !)
          */
         template <class LVectorType, class RVectorType>
-        void usmv(const double a, const LVectorType& x, RVectorType& b) const
+        void usmv([[maybe_unused]] const double a, [[maybe_unused]] const LVectorType& x, [[maybe_unused]] RVectorType& b) const
         {}
 
         /** \brief transposed Matrix-Vector multiplication with scalar multiplication
@@ -82,7 +84,7 @@ class NullOperator
          *  Implements b += a*N^tx and hence does nothing (N=0 !)
          */
         template <class LVectorType, class RVectorType>
-        void usmtv(const double a, const LVectorType& x, RVectorType& b) const
+        void usmtv([[maybe_unused]] const double a, [[maybe_unused]] const LVectorType& x, [[maybe_unused]] RVectorType& b) const
         {}
 
         /** \brief Matrix-Vector multiplication
@@ -90,25 +92,25 @@ class NullOperator
          *  Implements b = Nx and hence does nothing but set b=0 (N=0 !)
          */
         template <class LVectorType, class RVectorType>
-        void mv(const LVectorType& x, RVectorType& b) const
+        void mv([[maybe_unused]] const LVectorType& x, RVectorType& b) const
         {
             b = 0.0;
         }
 
         //! random access operator
-        const RowDummy& operator[](size_t i) const
+        const RowDummy& operator[]([[maybe_unused]] size_t i) const
         {
             return rowDummy_;
         }
 
         //! random access operator
-        RowDummy& operator[](size_t i)
+        RowDummy& operator[]([[maybe_unused]] size_t i)
         {
             return rowDummy_;
         }
 
         //! return j-th diagonal entry
-        const block_type& diagonal(size_t i) const
+        const block_type& diagonal([[maybe_unused]] size_t i) const
         {
             return rowDummy_[0];
         }

dune/solvers/operators/sumoperator.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef SUMOPERATOR_HH
 #define SUMOPERATOR_HH
 
 #include <cstddef>
+#include <dune/matrix-vector/resize.hh>
 
 /** \brief represents the sum of two linear operators
  *
@@ -33,8 +36,8 @@ class SumOperator
         //! construct from summands
         SumOperator(double a, SparseMatrixType& A, double b, LowRankMatrixType& M):
             alpha_(a),
-            sparse_matrix_(&A),
             beta_(b),
+            sparse_matrix_(&A),
             lowrank_matrix_(&M),
             summands_allocated_internally_(false)
         {}
@@ -42,8 +45,8 @@ class SumOperator
         //! construct from summands
         SumOperator(SparseMatrixType& A, LowRankMatrixType& M):
             alpha_(1.0),
-            sparse_matrix_(&A),
             beta_(1.0),
+            sparse_matrix_(&A),
             lowrank_matrix_(&M),
             summands_allocated_internally_(false)
         {}
@@ -137,5 +140,14 @@ class SumOperator
         const bool summands_allocated_internally_;
 };
 
+//! inject the resize from SumOperator to generic vector tools
+namespace Dune { namespace MatrixVector {
+template <class Vector, class SparseMatrix, class LowRankMatrix>
+//! Resize vector from SumOperator. Note that we only consider the sparse part.
+void resize(Vector& v, const SumOperator<SparseMatrix, LowRankMatrix>& sumOp) {
+  resize(v, sumOp.sparseMatrix());
+}
+}}
+
 #endif
 

dune/solvers/solvers/CMakeLists.txt

 install(FILES
-    cgsolver.cc
-    cgsolver.hh
+    criterion.hh
     iterativesolver.cc
     iterativesolver.hh
+    linearipopt.hh
+    linearsolver.hh
     loopsolver.cc
     loopsolver.hh
+    proximalnewtonsolver.hh
     quadraticipopt.hh
     solver.hh
     tcgsolver.cc
     tcgsolver.hh
+    trustregionsolver.cc
+    trustregionsolver.hh
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers/solvers)

dune/solvers/solvers/Makefile.am

-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
-
-EXTRA_DIST = CMakeLists.txt
-
-include $(top_srcdir)/am/global-rules
\ No newline at end of file

dune/solvers/solvers/cgsolver.cc

-#include <cmath>
-#include <limits>
-#include <iomanip>
-
-template <class MatrixType, class VectorType>
-void CGSolver<MatrixType, VectorType>::check() const
-{
-    if (preconditioner_)
-        preconditioner_->check();
-
-    if (!errorNorm_)
-        DUNE_THROW(SolverError, "You need to supply a norm to a CG solver!");
-
-}
-
-template <class MatrixType, class VectorType>
-void CGSolver<MatrixType, VectorType>::solve()
-{
-
-    int i;
-
-    VectorType& b = *rhs_;
-
-    // Check whether the solver is set up properly
-    check();
-
-    // set up preconditioner
-    if (preconditioner_)
-        preconditioner_->preprocess();
-
-    if (this->verbosity_ != NumProc::QUIET)
-        std::cout << "--- CGSolver ---\n";
-
-    if (this->verbosity_ == NumProc::FULL)
-    {
-        std::cout << " iter";
-        std::cout << "          error";
-        std::cout << "     rate";
-        std::string header;
-        if (preconditioner_) {
-            header = preconditioner_->getOutput();
-            std::cout << header;
-        }
-        std::cout << std::endl;
-        std::cout << "-----";
-        std::cout << "---------------";
-        std::cout << "---------";
-        for(size_t i=0; i<header.size(); ++i)
-            std::cout << "-";
-        std::cout << std::endl;
-    }
-
-    double error = std::numeric_limits<double>::max();
-
-    double normOfOldCorrection = 0;
-    double totalConvRate = 1;
-    this->maxTotalConvRate_ = 0;
-    int convRateCounter = 0;
-
-    // overwrite b with residual
-    matrix_->mmv(*x_,b);
-
-    VectorType p(*x_);              // the search direction
-    VectorType q(*x_);              // a temporary vector
-    
-    // some local variables
-    double rho,rholast,lambda,alpha,beta;
-    
-    // determine initial search direction
-    p = 0;                          // clear correction
-    
-    if (preconditioner_) {
-        preconditioner_->setProblem(*matrix_,p,b);
-        preconditioner_->iterate();
-        p = preconditioner_->getSol();
-    } else
-        p = b;
-
-    rholast = p*b;         // orthogonalization
-    
-    // Loop until desired tolerance or maximum number of iterations is reached
-    for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++) {
-        
-        // Backup of the current solution for the error computation later on
-        VectorType oldSolution = *x_;
-        
-        // Perform one iteration step
-        // minimize in given search direction p
-        matrix_->mv(p,q);       // q=Ap
-        alpha = p*q;            // scalar product
-        lambda = rholast/alpha; // minimization
-        x_->axpy(lambda,p);     // update solution
-        b.axpy(-lambda,q);      // update residual
-
-        // determine new search direction
-        q = 0;                  // clear correction
-        
-        // apply preconditioner
-        if (preconditioner_) {
-            preconditioner_->setProblem(*matrix_,q,b);
-            preconditioner_->iterate();
-            q = preconditioner_->getSol();
-        } else
-            q = b;
-        
-        rho = q*b;              // orthogonalization
-        beta = rho/rholast;     // scaling factor
-        p *= beta;              // scale old search direction
-        p += q;                 // orthogonalization with correction
-        rholast = rho;          // remember rho for recurrence
-
-        // write iteration to file, if requested
-        if (this->historyBuffer_!="") 
-            this->writeIterate(*x_, i);
-
-        // Compute error        
-        double oldNorm = errorNorm_->operator()(oldSolution);
-        oldSolution -= *x_;
-
-        double normOfCorrection = errorNorm_->operator()(oldSolution);
-
-        if (this->useRelativeError_)
-            error = normOfCorrection / oldNorm;
-        else
-            error = normOfCorrection;
-
-        double convRate = normOfCorrection / normOfOldCorrection;
-
-        normOfOldCorrection = normOfCorrection;
-
-        if (!isinf(convRate) && !isnan(convRate)) {
-            totalConvRate *= convRate;
-            this->maxTotalConvRate_ = std::max(this->maxTotalConvRate_, std::pow(totalConvRate, 1/((double)convRateCounter+1)));
-            convRateCounter++;
-        }
-
-        // Output
-        if (this->verbosity_ == NumProc::FULL) {
-            std::cout << std::setw(5) << i;
-
-            std::cout << std::setiosflags(std::ios::scientific);
-            std::cout << std::setw(15) << std::setprecision(7) << error;
-            std::cout << std::resetiosflags(std::ios::scientific);
-
-            std::cout << std::setiosflags(std::ios::fixed);
-            std::cout << std::setw(9) << std::setprecision(5) << convRate;
-            std::cout << std::resetiosflags(std::ios::fixed);
-
-            std::cout << std::endl;
-        }
-
-    }
-
-
-    if (this->verbosity_ != NumProc::QUIET) {
-      std::cout << "maxTotalConvRate: " << this->maxTotalConvRate_ << ",   "
-                  << i << " iterations performed\n";
-        std::cout << "--------------------\n";
-    }
-
-}

dune/solvers/solvers/cgsolver.hh

-#ifndef CG_SOLVER_HH
-#define CG_SOLVER_HH
-
-#include <dune/solvers/solvers/iterativesolver.hh>
-#include <dune/solvers/iterationsteps/lineariterationstep.hh>
-#include <dune/solvers/norms/norm.hh>
-
-/** \brief A conjugate gradient solver 
- *
- */
-template <class MatrixType, class VectorType>
-class CGSolver : public IterativeSolver<VectorType>
-{
-    static const int blocksize = VectorType::block_type::dimension;
-
-public:     
-
-    /** \brief Constructor taking all relevant data */
-        
-    CGSolver(const MatrixType* matrix,
-             VectorType* x,
-             VectorType* rhs,
-             LinearIterationStep<MatrixType,VectorType>* preconditioner,
-             int maxIterations,
-             double tolerance,
-             Norm<VectorType>* errorNorm,
-             NumProc::VerbosityMode verbosity,
-             bool useRelativeError=true)
-      : IterativeSolver<VectorType>(tolerance,maxIterations,verbosity,useRelativeError),
-          matrix_(matrix), x_(x), rhs_(rhs),
-          preconditioner_(preconditioner), errorNorm_(errorNorm)
-    {}
-    
-    /** \brief Loop, call the iteration procedure
-     * and monitor convergence */
-    virtual void solve();
-    
-    /** \brief Checks whether all relevant member variables are set
-     * \exception SolverError if the iteration step is not set up properly
-     */
-    virtual void check() const;
-
-    const MatrixType* matrix_;
-
-    VectorType* x_;
-
-    VectorType* rhs_;
-    
-    //! The iteration step used by the algorithm
-    LinearIterationStep<MatrixType,VectorType>* preconditioner_;
-    
-    //! The norm used to measure convergence
-    Norm<VectorType>* errorNorm_;
-    
-};
-
-#include "cgsolver.cc"
-
-#endif

dune/solvers/solvers/cholmodsolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_CHOLMODSOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_CHOLMODSOLVER_HH
+
+/** \file
+    \brief A wrapper for the CHOLMOD sparse direct solver
+*/
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/solver.hh>
+#include <dune/istl/cholmod.hh>
+#include <dune/istl/foreach.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/solvers/linearsolver.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
+
+namespace Dune
+{
+
+namespace Solvers
+{
+
+/** \brief Wraps the CHOLMOD sparse symmetric direct solver
+*
+*  CHOLMOD uses the Cholesky decomposition A = LL^T for sparse
+*  symmetric matrices to solve linear systems Ax = b efficiently.
+*
+*  This class wraps the dune-istl CHOLMOD solver for dune-solvers to provide
+*  a Solvers::LinearSolver interface.
+*/
+template <class MatrixType, class VectorType, class BitVectorType = DefaultBitVector_t<VectorType>>
+class CholmodSolver
+: public LinearSolver<MatrixType,VectorType>, public CanIgnore<BitVectorType>
+{
+  using field_type = typename VectorType::field_type;
+
+public:
+
+  /** \brief Default constructor */
+  CholmodSolver ()
+  : LinearSolver<MatrixType,VectorType>(NumProc::FULL)
+  {
+    // Bug in LibScotchMetis:
+    // metis_dswitch: see cholmod.h: LibScotchMetis throws segmentation faults for large problems.
+    // This caused a hard to understand CI problems and therefore, METIS is disabled for problems
+    // larger than defined in metis_nswitch for now (default 3000).
+    // Bug report: https://gitlab.inria.fr/scotch/scotch/issues/8
+    cholmodCommonObject().metis_dswitch = 0.0;
+  }
+
+  /** \brief Constructor for a linear problem */
+  CholmodSolver (const MatrixType& matrix,
+                 VectorType& x,
+                 const VectorType& rhs,
+                 NumProc::VerbosityMode verbosity=NumProc::FULL)
+  : LinearSolver<MatrixType,VectorType>(verbosity),
+    matrix_(&matrix),
+    x_(&x),
+    rhs_(&rhs)
+  {
+    // Bug in LibScotchMetis:
+    // metis_dswitch: see cholmod.h: LibScotchMetis throws segmentation faults for large problems.
+    // This caused a hard to understand CI problems and therefore, METIS is disabled for problems
+    // larger than defined in metis_nswitch for now (default 3000).
+    // Bug report: https://gitlab.inria.fr/scotch/scotch/issues/8
+    cholmodCommonObject().metis_dswitch = 0.0;
+  }
+
+  /** \brief Set the linear problem to solve
+   *
+   * The matrix is assumed to be symmetric! CHOLMOD uses only the
+   * upper part of the matrix, the lower part is ignored and can be
+   * empty for optimal memory use.
+   */
+  void setProblem(const MatrixType& matrix,
+                  VectorType& x,
+                  const VectorType& rhs) override
+  {
+    setMatrix(matrix);
+    setSolutionVector(x);
+    setRhs(rhs);
+  }
+
+  /** \brief Set the matrix for the linear linear problem to solve
+   *
+   * The matrix is assumed to be symmetric! CHOLMOD uses only the
+   * upper part of the matrix, the lower part is ignored and can be
+   * empty for optimal memory usage.
+   *
+   * \warning Unlike the setMatrix method of the dune-istl Cholmod class,
+   * the method here does not factorize the matrix!   Call the factorize
+   * method for that.
+   */
+  void setMatrix(const MatrixType& matrix)
+  {
+    matrix_ = &matrix;
+    isFactorized_ = false;
+  }
+
+  /** \brief Set the rhs vector for the linear problem
+   */
+  void setRhs(const VectorType& rhs)
+  {
+    rhs_ = &rhs;
+  }
+
+  /** \brief Set the vector where to write the solution of the linear problem
+   */
+  void setSolutionVector(VectorType& x)
+  {
+    x_ = &x;
+  }
+
+  /** \brief Compute the Cholesky decomposition of the matrix
+   *
+   * You must have set a matrix to be able to call this,
+   * either by calling setProblem or the appropriate constructor.
+   *
+   * The degrees of freedom to ignore must have been set before
+   * calling the `factorize` method.
+   */
+  void factorize()
+  {
+    if (not this->hasIgnore())
+    {
+      // setMatrix will decompose the matrix
+      istlCholmodSolver_.setMatrix(*matrix_);
+
+      // get the error code from Cholmod in case something happened
+      errorCode_ = cholmodCommonObject().status;
+    }
+    else
+    {
+      const auto& ignore = this->ignore();
+
+      // The setMatrix method stores only the selected entries of the matrix (determined by the ignore field)
+      istlCholmodSolver_.setMatrix(*matrix_,&ignore);
+      // get the error code from Cholmod in case something happened
+      errorCode_ = cholmodCommonObject().status;
+    }
+
+    if (errorCode_ >= 0)  // okay or warning, but no error
+      isFactorized_ = true;
+  }
+
+  /** \brief Solve the linear system
+   *
+   * This factorizes the matrix if it hasn't been factorized earlier
+   */
+  virtual void solve() override
+  {
+    // prepare the solver
+    Dune::InverseOperatorResult statistics;
+
+    // Factorize the matrix now, if the caller hasn't done so explicitly earlier.
+    if (!isFactorized_)
+      factorize();
+
+    if (not this->hasIgnore())
+    {
+      // the apply() method doesn't like constant values
+      auto mutableRhs = *rhs_;
+
+      // The back-substitution
+      istlCholmodSolver_.apply(*x_, mutableRhs, statistics);
+    }
+    else
+    {
+      const auto& ignore = this->ignore();
+
+      // create flat vectors
+      using T = typename VectorType::field_type;
+      std::size_t flatSize = flatVectorForEach(ignore, [](...){});
+      std::vector<bool> flatIgnore(flatSize);
+      std::vector<T> flatX(flatSize), flatRhsModifier(flatSize);
+
+      // define a lambda that copies the entries of a blocked vector into a flat one
+      auto copyToFlat = [&](auto&& blocked, auto&& flat)
+      {
+        flatVectorForEach(blocked, [&](auto&& entry, auto&& offset)
+        {
+          flat[offset] = entry;
+        });
+      };
+
+      // copy x and the ignore field for the modification of the rhs
+      copyToFlat(ignore,flatIgnore);
+      copyToFlat(*x_,flatX);
+
+      flatMatrixForEach( *matrix_, [&]( auto&& entry, auto&& rowOffset, auto&& colOffset){
+
+        // we assume entries only in the upper part and do nothing on the diagonal
+        if ( rowOffset >= colOffset )
+          return;
+
+        // upper part: if either col or row is ignored: move the entry or the symmetric duplicate to the rhs
+        if ( flatIgnore[colOffset] and not flatIgnore[rowOffset] )
+            flatRhsModifier[rowOffset] += entry * flatX[colOffset];
+
+        else if ( not flatIgnore[colOffset] and flatIgnore[rowOffset] )
+            flatRhsModifier[colOffset] += entry * flatX[rowOffset];
+      });
+
+      // update the rhs
+      auto modifiedRhs = *rhs_;
+      flatVectorForEach(modifiedRhs, [&](auto&& entry, auto&& offset){
+        if ( not flatIgnore[offset] )
+          entry -= flatRhsModifier[offset];
+      });
+
+      // Solve the modified system
+      istlCholmodSolver_.apply(*x_, modifiedRhs, statistics);
+    }
+  }
+
+  cholmod_common& cholmodCommonObject()
+  {
+    return istlCholmodSolver_.cholmodCommonObject();
+  }
+
+  /** \brief return the error code of the Cholmod factorize call
+   *
+   * In setMatrix() the matrix factorization takes place and Cholmod is
+   * able to communicate error codes of the factorization in the status
+   * field of the cholmodCommonObject.
+   * The return value 0 means "good" and the other values can be found
+   * in the Cholmod User Guide.
+   */
+  int errorCode() const
+  {
+    return errorCode_;
+  }
+
+  //! dune-istl Cholmod solver object
+  Dune::Cholmod<VectorType> istlCholmodSolver_;
+
+  //! Pointer to the system matrix
+  const MatrixType* matrix_;
+
+  //! Pointer to the solution vector
+  VectorType* x_;
+
+  //! Pointer to the right hand side
+  const VectorType* rhs_;
+
+  //! error code of Cholmod factorize call
+  int errorCode_ = 0;
+
+  /** \brief Whether istlCholmodSolver_ currently holds a factorized matrix.
+   *
+   * isFactorized==true is equivalent to istlCholmodSolver_->L_ != nullptr,
+   * but the latter information is private to the dune-istl Cholmod class.
+   */
+  bool isFactorized_ = false;
+};
+
+}   // namespace Solvers
+}   // namespace Dune
+
+#endif

dune/solvers/solvers/criterion.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=8 sw=2 sts=2:
+
+#include <config.h>
+
+#include <tuple>
+
+#include <dune/solvers/solvers/criterion.hh>
+
+
+Dune::Solvers::Criterion::Result
+Dune::Solvers::Criterion::operator()() const
+{
+  return f_();
+}
+
+std::string
+Dune::Solvers::Criterion::header() const
+{
+  return header_;
+}
+
+Dune::Solvers::Criterion
+Dune::Solvers::operator| (const Criterion& c1, const Criterion& c2)
+{
+  return Criterion(
+      [=]() {
+        auto r1 = c1();
+        auto r2 = c2();
+        return std::make_tuple(std::get<0>(r1) or std::get<0>(r2), std::get<1>(r1).append(std::get<1>(r2)));
+      },
+      c1.header().append(c2.header()));
+}
+
+Dune::Solvers::Criterion
+Dune::Solvers::operator& (const Criterion& c1, const Criterion& c2)
+{
+  return Criterion(
+      [=]() {
+        auto r1 = c1();
+        auto r2 = c2();
+        return std::make_tuple(std::get<0>(r1) and std::get<0>(r2), std::get<1>(r1).append(std::get<1>(r2)));
+      },
+      c1.header().append(c2.header()));
+}
+
+Dune::Solvers::Criterion
+Dune::Solvers::operator~ (const Criterion& c)
+{
+  return Criterion(
+      [=]() {
+        auto r = c();
+        return std::make_tuple(not(std::get<0>(r)), std::get<1>(r));
+      },
+      c.header());
+}
+

dune/solvers/solvers/criterion.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=8 sw=2 sts=2:
+#ifndef DUNE_SOLVERS_SOLVERS_CRITERION_HH
+#define DUNE_SOLVERS_SOLVERS_CRITERION_HH
+
+
+#include <functional>
+#include <string>
+#include <tuple>
+#include <memory>
+#include <type_traits>
+
+#include <dune/common/stringutility.hh>
+
+
+/**
+ * \attention Everything contained in this headere is experimental and may change in the near future.
+ */
+
+namespace Dune {
+
+  namespace Solvers {
+
+    /**
+     * \brief Class to wrap termination criteria
+     *
+     * Each criterion is a functor returning a std::tuple<bool,string>
+     * where the bool encodes if the criterion is satisfied whereas
+     * the string will be written to the logging output.
+     */
+    class Criterion
+    {
+      typedef std::tuple<bool, std::string> Result;
+
+    public:
+
+      /**
+       * \brief Create Criterion from function and header string
+       *
+       * \param f A functor returning a tuple<bool,string>
+       * \param header The header string describing the output of this criterion
+       *
+       * The bool of the tuple returned by the fuctor is interpreted
+       * as 'criterion satisfied' whereas the string will be written to
+       * the log output.
+       *
+       * Be aware that this stores a copy of the provided functor.
+       * Hence it should be lightweight and avoid holding large
+       * data sets or wrapped using std::cref.
+       *
+       * Take care that the length of returned strings and the length
+       * of the header string coincide to get a readable log.
+       */
+      template<class F,
+        std::enable_if_t<std::is_convertible<std::invoke_result_t<F>, Result>::value, int> = 0>
+      Criterion(F&& f, const std::string& header) :
+        f_(std::forward<F>(f)),
+        header_(header)
+      {}
+
+      /**
+       * \brief Create always false Criterion for logging strings
+       *
+       * \param f A functor returning a string
+       * \param header The header string describing the output of this criterion
+       *
+       * The string returned by the functor will be written to
+       * the log output while criterion is interpreted as always
+       * false.
+       *
+       * Be aware that this stores a copy of the provided functor.
+       * Hence it should be lightweight and avoid holding large
+       * data sets or wrapped using std::cref.
+       *
+       * Take care that the length of returned strings and the length
+       * of the header string coincide to get a readable log.
+       */
+      template<class F,
+        std::enable_if_t<std::is_convertible<std::invoke_result_t<F>, std::string>::value, int> = 0>
+      Criterion(F&& f, const std::string& header) :
+        f_( [=]() mutable {return std::make_tuple(false, f());} ),
+        header_(header)
+      {}
+
+      /**
+       * \brief Create always false Criterion for logging any printf-like type
+       *
+       * \param f A functor returning some type
+       * \param format A printf format string
+       * \param header The header string describing the output of this criterion
+       *
+       * The valued returned by the functor will be converted using the
+       * format string and written to the log output while criterion is
+       * interpreted as always false.
+       *
+       * Be aware that this stores a copy of the provided functor.
+       * Hence it should be lightweight and avoid holding large
+       * data sets or wrapped using std::cref.
+       *
+       * Take care that the length of returned strings and the length
+       * of the header string coincide to get a readable log.
+       */
+      template<class F>
+      Criterion(F&& f, const std::string& format, const std::string& header) :
+        f_( [=]() {return std::make_tuple(false, formatString(format, f()));} ),
+        header_(header)
+      {}
+
+      /**
+       * \brief Check criterion
+       *
+       * \returns A tuple<bool,string>
+       */
+      Result operator()() const;
+
+      /**
+       * \brief Obtain header string
+       */
+      std::string header() const;
+
+    protected:
+      std::function<Result()> f_;
+      std::string header_;
+    }; // end of class Criterion
+
+
+
+    // free helper functions for logical operations on criteria
+
+    /**
+     * \brief Create Criterion that combines two others using 'or'
+     *
+     * The log output of both criteria will be concatenated.
+     *
+     * Notice that this will store copies of both provided criteria.
+     */
+    Criterion operator| (const Criterion& c1, const Criterion& c2);
+
+    /**
+     * \brief Create Criterion that combines two others using 'and'
+     *
+     * The log output of both criteria will be concatenated.
+     *
+     * Notice that this will store copies of both provided criteria.
+     */
+    Criterion operator& (const Criterion& c1, const Criterion& c2);
+
+    /**
+     * \brief Create Criterion that negates another
+     *
+     * The log output will be forwarded.
+     *
+     * Notice that this will store a copy of the provided criterion.
+     */
+    Criterion operator~ (const Criterion& c);
+
+
+
+
+    // free helper functions for generating some classic criteria
+
+    /**
+     * \brief Create a Criterion for limiting the number of iteration steps
+     *
+     * \param solver The criterion will query this solver for its current iteration step.
+     * \param maxIterations The criterion will be true once this number of iteration steps is reached.
+     */
+    template<class SolverType>
+    Criterion maxIterCriterion(const SolverType& solver, int maxIterations)
+    {
+        return Criterion(
+                [&solver, maxIterations](){
+                    int i = solver.iterationCount();
+                    return std::make_tuple( i>=maxIterations, Dune::formatString("% 7d", i) );
+                },
+                "   iter");
+    }
+
+    /**
+     * \brief Create a Criterion for checking the correction norm
+     *
+     * \param iterationStep The IterationStep to be monitored
+     * \param norm          Norm to be evaluated for the correction
+     * \param tolerance     Tolerance necessary to satisfy the criterion
+     * \param correctionNorms A container to store the computed correction norms
+     */
+    template<class IterationStepType, class NormType, class CorrectionNormContainer>
+    Criterion correctionNormCriterion(
+          const IterationStepType& iterationStep,
+          const NormType& norm,
+          double tolerance,
+          CorrectionNormContainer& correctionNorms)
+    {
+      auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+      return Criterion(
+          [&, lastIterate, tolerance] () mutable {
+            double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+            correctionNorms.push_back(normOfCorrection);
+            *lastIterate = *iterationStep.getIterate();
+            return std::make_tuple(normOfCorrection < tolerance, Dune::formatString(" % 14.7e", normOfCorrection));
+          },
+          " abs correction");
+    }
+
+    /**
+     * \brief Create a Criterion for checking the correction norm
+     *
+     * \param iterationStep The IterationStep to be monitored
+     * \param norm          Norm to be evaluated for the correction
+     * \param tolerance     Tolerance necessary to satisfy the criterion
+     */
+    template<class IterationStepType, class NormType>
+    Criterion correctionNormCriterion(
+          const IterationStepType& iterationStep,
+          const NormType& norm,
+          double tolerance)
+    {
+      auto lastIterate = std::make_shared<typename IterationStepType::Vector>(*iterationStep.getIterate());
+      return Criterion(
+          [&, lastIterate, tolerance] () mutable {
+            double normOfCorrection = norm.diff(*lastIterate, *iterationStep.getIterate());
+            *lastIterate = *iterationStep.getIterate();
+            return std::make_tuple(normOfCorrection < tolerance, Dune::formatString(" % 14.7e", normOfCorrection));
+          },
+          " abs correction");
+    }
+
+  } // end namespace Solvers
+
+} // end namespace Dune
+
+
+#endif
+
+

dune/solvers/solvers/iterativesolver.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+
 #include <cmath>
 #include <limits>
 #include <iostream>
 #include <iomanip>
 #include <fstream>
 
-#include <dune/solvers/common/genericvectortools.hh>
+#include <dune/matrix-vector/genericvectortools.hh>
+
+namespace Dune {
+namespace Solvers {
 
 template <class VectorType, class BitVectorType>
 void IterativeSolver<VectorType, BitVectorType>::check() const
 {
+    if (!iterationStep_)
+        DUNE_THROW(SolverError, "You need to supply an iteration step to an iterative solver!");
+
+    iterationStep_->check();
+
     if (!errorNorm_)
         DUNE_THROW(SolverError, "You need to supply a norm-computing routine to an iterative solver!");
 }
@@ -18,14 +29,17 @@ void IterativeSolver<VectorType, BitVectorType>::writeIterate(const VectorType&
                                                               int iterationNumber) const
 {
     std::stringstream iSolFilename;
-    iSolFilename << historyBuffer_ << "/intermediatesolution_" 
+    iSolFilename << historyBuffer_ << "/intermediatesolution_"
                  << std::setw(4) << std::setfill('0') << iterationNumber;
-    
+
     std::ofstream file(iSolFilename.str().c_str(), std::ios::out|std::ios::binary);
     if (not(file))
         DUNE_THROW(SolverError, "Couldn't open file " << iSolFilename.str() << " for writing");
-    
-    GenericVector::writeBinary(file, iterate);
-    
+
+    Dune::MatrixVector::Generic::writeBinary(file, iterate);
+
     file.close();
 }
+
+} /* namespace Solvers */
+} /* namespace Dune */

dune/solvers/solvers/iterativesolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef DUNE_SOLVERS_ITERATIVE_SOLVER_HH
 #define DUNE_SOLVERS_ITERATIVE_SOLVER_HH
 
 #include <dune/common/ftraits.hh>
+#include <dune/common/shared_ptr.hh>
 
+#include <dune/solvers/common/defaultbitvector.hh>
+#include <dune/solvers/common/wrapownshare.hh>
 #include <dune/solvers/solvers/solver.hh>
 #include <dune/solvers/iterationsteps/iterationstep.hh>
 #include <dune/solvers/norms/norm.hh>
 
+namespace Dune {
+
+  namespace Solvers {
+
     /** \brief Abstract base class for iterative solvers */
-    template <class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
+    template <class VectorType, class BitVectorType = DefaultBitVector_t<VectorType> >
     class IterativeSolver : public Solver
     {
-        typedef typename VectorType::value_type::field_type field_type;
-
-        // For complex-valued data
-        typedef typename Dune::FieldTraits<field_type>::real_type real_type;
+        // For norms and convergence rates
+        using real_type =  typename Dune::FieldTraits<VectorType>::real_type;
+        using ItStep = IterationStep<VectorType, BitVectorType>;
+        using ErrorNorm = Norm<VectorType>;
 
     public:
 
@@ -25,12 +34,13 @@
                         bool useRelativeError=true)
             : Solver(verbosity),
               tolerance_(tolerance),
-              maxIterations_(maxIterations), errorNorm_(NULL),
+              maxIterations_(maxIterations), errorNorm_(nullptr),
               historyBuffer_(""),
               useRelativeError_(useRelativeError)
         {}
 
         /** \brief Constructor taking all relevant data */
+        [[deprecated("Handing over raw pointer in the constructor is deprecated!")]]
         IterativeSolver(int maxIterations,
                         double tolerance,
                         const Norm<VectorType>* errorNorm,
@@ -38,7 +48,23 @@
                         bool useRelativeError=true)
             : Solver(verbosity),
               tolerance_(tolerance),
-              maxIterations_(maxIterations), errorNorm_(errorNorm),
+              maxIterations_(maxIterations),
+              errorNorm_(Dune::stackobject_to_shared_ptr(*errorNorm)),
+              historyBuffer_(""),
+              useRelativeError_(useRelativeError)
+        {}
+
+        /** \brief Constructor taking all relevant data */
+        template <class NormT, class Enable = std::enable_if_t<not std::is_pointer<NormT>::value> >
+        IterativeSolver(int maxIterations,
+                        double tolerance,
+                        NormT&& errorNorm,
+                        VerbosityMode verbosity,
+                        bool useRelativeError=true)
+            : Solver(verbosity),
+              tolerance_(tolerance),
+              maxIterations_(maxIterations),
+              errorNorm_(wrap_own_share<const ErrorNorm>(std::forward<NormT>(errorNorm))),
               historyBuffer_(""),
               useRelativeError_(useRelativeError)
         {}
@@ -54,15 +80,52 @@
         /** \brief Write the current iterate to disk (for convergence measurements) */
         void writeIterate(const VectorType& iterate, int iterationNumber) const;
 
+        /** \brief Set iteration step */
+        template <class Step>
+        void setIterationStep(Step&& iterationStep)
+        {
+            iterationStep_ = wrap_own_share<ItStep>(std::forward<Step>(iterationStep));
+        }
+
+        /** \brief Get iteration step */
+        const ItStep& getIterationStep() const
+        {
+          return *iterationStep_;
+        }
+
+        /** \brief Get iteration step */
+        ItStep& getIterationStep()
+        {
+          return *iterationStep_;
+        }
+
+        /** \brief Set the error norm */
+        template <class NormT>
+        void setErrorNorm(NormT&& errorNorm)
+        {
+            errorNorm_ = wrap_own_share<ErrorNorm>(std::forward<NormT>(errorNorm));
+        }
+
+        /** \brief Get error norm */
+        const ErrorNorm& getErrorNorm() const
+        {
+          return *errorNorm_;
+        }
+
         /** \brief The requested tolerance of the solver */
         double tolerance_;
 
         //! The maximum number of iterations
         int maxIterations_;
 
+    protected:
+        //! The iteration step used by the algorithm
+        std::shared_ptr<ItStep> iterationStep_;
+
         //! The norm used to measure convergence
-        const Norm<VectorType>* errorNorm_;
+        std::shared_ptr<const ErrorNorm> errorNorm_;
 
+    public:
         /** \brief If this string is nonempty it is expected to contain a valid
             directory name.  All intermediate iterates are then written there. */
         std::string historyBuffer_;
@@ -73,6 +136,13 @@
 
     };
 
+  }   // namespace Solvers
+
+}   // namespace Dune
+
+// For backward compatibility: will be removed eventually
+using Dune::Solvers::IterativeSolver;
+
 #include "iterativesolver.cc"
 
 #endif

dune/solvers/solvers/linearipopt.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_LINEAR_IPOPT_SOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_LINEAR_IPOPT_SOLVER_HH
+
+/** \file
+    \brief A wrapper for the IPOpt solver for linear optimization problem
+    with box constraints and linear constraints
+*/
+#include <algorithm>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/shared_ptr.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+
+#include "IpTNLP.hpp"
+#include "IpIpoptApplication.hpp"
+
+/** \brief Implementation class used to pipe linear problems to
+    the IPOpt interior-point solver
+
+    \tparam JacobianType The type of the jacobian of the linear constraints
+*/
+template <class VectorType, class JacobianType>
+class LinearIPOptProblem : public Ipopt::TNLP
+{
+    enum {blocksize = VectorType::block_type::dimension};
+    enum {constraintBlocksize = JacobianType::block_type::rows};
+
+    typedef typename VectorType::field_type field_type;
+
+public:
+    /** Setup problem if no linear constraints are present */
+    LinearIPOptProblem( VectorType* x, const VectorType* rhs,
+                          const Dune::BitSetVector<blocksize>* dirichletNodes,
+                          std::vector<BoxConstraint<field_type, blocksize> >* obstacles)
+        : x_(x), rhs_(rhs),
+          dirichletNodes_(dirichletNodes), obstacles_(obstacles),
+          constraintMatrix_(nullptr), constraintObstacles_(nullptr), jacobianCutoff_(1e-8)
+    {}
+
+    /** Setup problem full problem */
+    LinearIPOptProblem(VectorType* x, const VectorType* rhs,
+                          const Dune::BitSetVector<blocksize>* dirichletNodes,
+                          std::vector<BoxConstraint<field_type, blocksize> >* boundObstacles,
+                          std::shared_ptr<const JacobianType> constraintMatrix,
+                          std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles)
+        : x_(x), rhs_(rhs),
+          dirichletNodes_(dirichletNodes), obstacles_(boundObstacles),
+          constraintMatrix_(constraintMatrix), constraintObstacles_(constraintObstacles),
+          jacobianCutoff_(1e-8)
+    {}
+
+  /** default destructor */
+    virtual ~LinearIPOptProblem() {}
+
+  /**@name Overloaded from TNLP */
+  //@{
+  /** Method to return some info about the nlp */
+  virtual bool get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+                            Ipopt::Index& nnz_h_lag, IndexStyleEnum& index_style);
+
+  /** Method to return the bounds for my problem */
+  virtual bool get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+                               Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u);
+
+  /** Method to return the starting point for the algorithm */
+  virtual bool get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+                                  bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+                                  Ipopt::Index m, bool init_lambda,
+                                  Ipopt::Number* lambda);
+
+  /** Method to return the objective value */
+  virtual bool eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value);
+
+  /** Method to return the gradient of the objective */
+  virtual bool eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f);
+
+  /** Method to return the constraint residuals */
+  virtual bool eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g);
+
+  /** Method to return:
+   *   1) The structure of the jacobian (if "values" is NULL)
+   *   2) The values of the jacobian (if "values" is not NULL)
+   */
+  virtual bool eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+                          Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+                          Ipopt::Number* values);
+
+  /** Method to return:
+   *   1) The structure of the hessian of the lagrangian (if "values" is NULL)
+   *   2) The values of the hessian of the lagrangian (if "values" is not NULL)
+   */
+  virtual bool eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+                      Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+                      bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+                      Ipopt::Index* jCol, Ipopt::Number* values);
+
+  //@}
+
+  /** @name Solution Methods */
+  //@{
+  /** This method is called when the algorithm is complete so the TNLP can store/write the solution */
+    virtual void finalize_solution(Ipopt::SolverReturn status,
+                                 Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+                                 Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+                                 Ipopt::Number obj_value,
+                                   const Ipopt::IpoptData* ip_data,
+                                   Ipopt::IpoptCalculatedQuantities* ip_cq);
+  //@}
+
+    // /////////////////////////////////
+    //   Data
+    // /////////////////////////////////
+
+    //! Vector to store the solution
+    VectorType* x_;
+
+    //! The linear term of the linear energy
+    const VectorType* rhs_;
+
+    const Dune::BitSetVector<blocksize>* dirichletNodes_;
+
+    //! Vector containing the bound constraints of the variables
+    // Can stay unset when no bound constraints exist
+    std::vector<BoxConstraint<field_type, blocksize> >* obstacles_;
+
+    //! The matrix corresponding to linear constraints
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const JacobianType> constraintMatrix_;
+
+    //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const std::vector<BoxConstraint<field_type, constraintBlocksize> > > constraintObstacles_;
+
+private:
+    const field_type jacobianCutoff_;
+    /**@name Methods to block default compiler methods.*/
+    //@{
+    LinearIPOptProblem(const LinearIPOptProblem&);
+    LinearIPOptProblem& operator=(const LinearIPOptProblem&);
+    //@}
+};
+
+// returns the size of the problem
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType, JacobianType>::
+get_nlp_info(Ipopt::Index& n, Ipopt::Index& m, Ipopt::Index& nnz_jac_g,
+             Ipopt::Index& nnz_h_lag, Ipopt::TNLP::IndexStyleEnum& index_style)
+{
+    // use the C style indexing (0-based)
+    index_style = Ipopt::TNLP::C_STYLE;
+
+    // Number of variables
+    n = x_->dim();
+
+    // Number of constraints
+    // Dirichlet conditions are actually bound inequality constraints
+    // bound constraints are handled differently by IpOpt
+    m = 0;
+
+    // hence the constraint jacobian is empty
+    nnz_jac_g = 0;
+
+    // energy is linear, thus the hessian is zero
+    nnz_h_lag = 0;
+
+    // We only need the lower left corner of the Hessian (since it is symmetric)
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    assert(constraintMatrix_->N()==constraintObstacles_->size());
+
+    // only change the constraint information
+    m = constraintMatrix_->N()*constraintBlocksize;
+
+    int numJacobianNonzeros = 0;
+    for (size_t i=0; i<constraintMatrix_->N(); i++) {
+
+        const typename JacobianType::row_type& row = (*constraintMatrix_)[i];
+
+        // dim for  each entry in the constraint matrix
+        typename JacobianType::row_type::ConstIterator cIt    = row.begin();
+        typename JacobianType::row_type::ConstIterator cEndIt = row.end();
+        for (; cIt!=cEndIt; cIt++)
+            for (int k=0; k<constraintBlocksize; k++)
+                for (int l=0; l<blocksize; l++)
+                    if (std::abs((*cIt)[k][l]) > jacobianCutoff_ )
+                        numJacobianNonzeros++;
+    }
+
+    nnz_jac_g = numJacobianNonzeros;
+
+    return true;
+}
+
+
+// returns the variable bounds
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType, JacobianType>::
+get_bounds_info(Ipopt::Index n, Ipopt::Number* x_l, Ipopt::Number* x_u,
+                Ipopt::Index m, Ipopt::Number* g_l, Ipopt::Number* g_u)
+{
+    // here, the n and m we gave IPOPT in get_nlp_info are passed back to us.
+    // If desired, we could assert to make sure they are what we think they are.
+    assert(n == (Ipopt::Index)x_->dim());
+
+    if (obstacles_) {
+
+        for (size_t i=0; i<x_->size(); i++) {
+
+            for (size_t j=0; j<blocksize; j++) {
+
+                // Simply copy the values.  If they exceed a certain limit they are
+                // considered 'off'.
+                x_l[i*blocksize+j] = (*obstacles_)[i].lower(j);
+                x_u[i*blocksize+j] = (*obstacles_)[i].upper(j);
+
+            }
+
+        }
+
+    } else {
+
+        // unset all bounds
+        for (size_t i=0; i<x_->dim(); i++) {
+            x_l[i] = -std::numeric_limits<double>::max();
+            x_u[i] =  std::numeric_limits<double>::max();
+        }
+
+    }
+
+    if (dirichletNodes_) {
+
+        for (size_t i=0; i<x_->size(); i++) {
+
+            for (size_t j=0; j<blocksize; j++) {
+
+                if ( (*dirichletNodes_)[i][j])
+                    x_l[i*blocksize+j] = x_u[i*blocksize+j] = (*this->x_)[i][j];
+
+            }
+
+        }
+
+    }
+
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    // initialize with large numbers
+    for (int i=0; i<m; i++) {
+        g_l[i] = -std::numeric_limits<field_type>::max();
+        g_u[i] =  std::numeric_limits<field_type>::max();
+    }
+
+    for (size_t i=0; i<constraintObstacles_->size(); i++) {
+        for (int k=0; k<constraintBlocksize; k++) {
+
+                g_l[i*constraintBlocksize + k] = (*constraintObstacles_)[i].lower(k);
+                g_u[i*constraintBlocksize + k] = (*constraintObstacles_)[i].upper(k);
+        }
+    }
+
+    return true;
+}
+
+// returns the initial point for the problem
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+get_starting_point(Ipopt::Index n, bool init_x, Ipopt::Number* x,
+                   bool init_z, Ipopt::Number* z_L, Ipopt::Number* z_U,
+                   Ipopt::Index m, bool init_lambda, Ipopt::Number* lambda)
+{
+    // Here, we assume we only have starting values for x, if you code
+    // your own NLP, you can provide starting values for the dual variables
+    // if you wish
+    assert(init_x == true);
+    assert(init_z == false);
+    assert(init_lambda == false);
+
+    // initialize to the given starting point
+    for (size_t i=0; i<x_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            x[i*blocksize+j] = (*x_)[i][j];
+
+    return true;
+}
+
+// returns the value of the objective function
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number& obj_value)
+{
+    // Init return value
+    obj_value = 0;
+
+    ////////////////////////////////////
+    // Compute rhs*x
+    ////////////////////////////////////
+
+    for (size_t i=0; i<rhs_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            obj_value += x[i*blocksize+j] * (*rhs_)[i][j];
+
+  return true;
+}
+
+// return the gradient of the objective function grad_{x} f(x)
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_grad_f(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Number* grad_f)
+{
+    // g = g - b
+    for (size_t i=0; i<rhs_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            grad_f[i*blocksize+j] = (*rhs_)[i][j];
+
+  return true;
+}
+
+// return the value of the constraints: g(x)
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x, Ipopt::Index m, Ipopt::Number* g)
+{
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    for (int i=0; i<m; i++)
+        g[i] = 0;
+
+    int rowCounter = 0;
+    for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+        const auto& row = (*constraintMatrix_)[rowIdx];
+        auto cIt = row.begin();
+        auto cEnd = row.end();
+
+        // constraintMatrix times x
+        for (; cIt != cEnd; cIt++)
+            for (int l=0; l<constraintBlocksize; l++)
+                for (int k=0; k<blocksize; k++)
+                    if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+                        g[rowCounter+l] += (*cIt)[l][k] * x[cIt.index()*blocksize+k];
+
+        // increment constraint counter
+        rowCounter += constraintBlocksize;
+    }
+
+    return true;
+
+}
+
+// return the structure or values of the jacobian
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_jac_g(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+           Ipopt::Index m, Ipopt::Index nele_jac, Ipopt::Index* iRow, Ipopt::Index *jCol,
+           Ipopt::Number* values)
+{
+    if (constraintMatrix_==nullptr)
+        return true;
+
+    int idx = 0;
+    int rowCounter = 0;
+
+    // If the values are null then IpOpt wants the sparsity pattern
+    if (values==NULL) {
+
+        for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+            const auto& row = (*constraintMatrix_)[rowIdx];
+            auto cIt = row.begin();
+            auto cEnd = row.end();
+
+            for (; cIt != cEnd; ++cIt) {
+
+                for (int l=0; l<constraintBlocksize; l++)
+                    for (int k=0; k<blocksize; k++)
+                        if ( std::abs((*cIt)[l][k]) > jacobianCutoff_) {
+
+                            iRow[idx] = rowCounter + l;
+                            jCol[idx] = cIt.index()*blocksize+k;
+                            idx++;
+                        }
+            }
+
+            rowCounter += constraintBlocksize;
+        }
+
+    // If the values are not null IpOpt wants the evaluated constraints gradient
+    } else {
+
+        for (size_t rowIdx=0; rowIdx<constraintMatrix_->N(); rowIdx++) {
+
+            const auto& row = (*constraintMatrix_)[rowIdx];
+            auto cIt = row.begin();
+            auto cEnd = row.end();
+
+            for (; cIt != cEnd; ++cIt) {
+
+                for (int l=0; l<constraintBlocksize; l++)
+                    for (int k=0; k<blocksize; k++)
+                        if ( std::abs((*cIt)[l][k]) > jacobianCutoff_)
+                            values[idx++] = (*cIt)[l][k];
+            }
+
+            rowCounter += constraintBlocksize;
+
+        }
+
+    }
+
+    return true;
+
+}
+
+//return the structure or values of the hessian
+template <class VectorType, class JacobianType>
+bool LinearIPOptProblem<VectorType,JacobianType>::
+eval_h(Ipopt::Index n, const Ipopt::Number* x, bool new_x,
+       Ipopt::Number obj_factor, Ipopt::Index m, const Ipopt::Number* lambda,
+       bool new_lambda, Ipopt::Index nele_hess, Ipopt::Index* iRow,
+       Ipopt::Index* jCol, Ipopt::Number* values)
+{
+    return true;
+}
+
+template <class VectorType, class JacobianType>
+void LinearIPOptProblem<VectorType,JacobianType>::
+finalize_solution(Ipopt::SolverReturn status,
+                  Ipopt::Index n, const Ipopt::Number* x, const Ipopt::Number* z_L, const Ipopt::Number* z_U,
+                  Ipopt::Index m, const Ipopt::Number* g, const Ipopt::Number* lambda,
+                  Ipopt::Number obj_value,
+                  const Ipopt::IpoptData* ip_data,
+                  Ipopt::IpoptCalculatedQuantities* ip_cq)
+{
+  // here is where we would store the solution to variables, or write to a file, etc
+  // so we could use the solution.
+
+    for (size_t i=0; i<x_->size(); i++)
+        for (int j=0; j<blocksize; j++)
+            (*x_)[i][j] = x[i*blocksize+j];
+
+}
+
+/** \brief Wraps the IPOpt interior-point solver for linear functionals with
+ *    linear constraints and bound constraints.
+ *
+ *   The problems that can be solved are of the form
+ *   \f[
+ *       \min_x \frac{1}{2}x^T A x - b^T x
+ *               x_l \leq x \leq x_u
+ *               g_l \leq G x \leq g_u
+ *   \f]
+ *
+ * \note This is not a 'LinearSolver' in the dune-istl sense of the word,
+ *  because it minimizes a linear functional, rather than solving
+ *  a linear equation.
+ *
+ *  \tparam JacobianType The type of the jacobian of the linear constraints
+*/
+template <class VectorType, class JacobianType>
+class LinearIPOptSolver
+    : public IterativeSolver<VectorType>, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> >
+{
+
+    enum {blocksize = VectorType::block_type::dimension};
+    enum {constraintBlocksize = JacobianType::block_type::rows};
+
+    typedef typename VectorType::field_type field_type;
+
+    typedef Dune::FieldMatrix<field_type, blocksize, blocksize> MatrixBlock;
+    typedef Dune::FieldVector<field_type, blocksize>            VectorBlock;
+
+
+public:
+
+    /** \brief Default constructor */
+    LinearIPOptSolver () : IterativeSolver<VectorType>(1e-8, 100, NumProc::FULL),
+                              rhs_(NULL),
+                              obstacles_(NULL),
+                              linearSolverType_(""),
+                              constraintObstacles_(nullptr),
+                              constraintMatrix_(nullptr)
+    {}
+
+    /** \brief Constructor for a linear problem */
+    LinearIPOptSolver (VectorType& x,
+                       const VectorType& rhs,
+                       NumProc::VerbosityMode verbosity=NumProc::FULL,
+                       std::string linearSolverType = "")
+        : IterativeSolver<VectorType>(1e-8, 100, verbosity),
+          rhs_(&rhs), obstacles_(NULL),
+          linearSolverType_(linearSolverType),
+          constraintObstacles_(nullptr),
+          constraintMatrix_(nullptr)
+    {
+        this->x_ = &x;
+    }
+
+    void setProblem(VectorType& x,
+                    const VectorType& rhs) {
+        x_ = &x;
+        rhs_ = &rhs;
+    }
+
+    void setLinearConstraints(const JacobianType& constraintMatrix,
+                     const std::vector<BoxConstraint<field_type,constraintBlocksize> >& obstacles) {
+        constraintMatrix_ = Dune::stackobject_to_shared_ptr(constraintMatrix);
+        constraintObstacles_ = Dune::stackobject_to_shared_ptr(obstacles);
+    }
+
+    virtual void solve();
+
+    ///////////////////////////////////////////////////////
+    // Data
+    ///////////////////////////////////////////////////////
+
+    //! Vector to store the solution
+    VectorType* x_;
+
+    //! The linear term in the linear energy
+    const VectorType* rhs_;
+
+    //! Vector containing the bound constraints of the variables
+    // Can stay unset when no bound constraints exist
+    std::vector<BoxConstraint<field_type,blocksize> >* obstacles_;
+
+    //! The type of the linear solver to be used within IpOpt, default is ""
+    std::string linearSolverType_;
+
+    //! IpOpt expects constraints to be of the type g_l <= g(x) <= g_u
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const std::vector<BoxConstraint<field_type,constraintBlocksize> > > constraintObstacles_;
+
+    //! The matrix corresponding to linear constraints
+    // Can stay unset when no linear constraints exist
+    std::shared_ptr<const JacobianType> constraintMatrix_;
+};
+
+template <class VectorType, class JacobianType>
+void LinearIPOptSolver<VectorType,JacobianType>::solve()
+{
+  // Create a new instance of your nlp
+    Ipopt::SmartPtr<Ipopt::TNLP> mynlp = new LinearIPOptProblem<VectorType,JacobianType>(x_, rhs_,
+                                                           this->ignoreNodes_, obstacles_,
+                                                            constraintMatrix_,constraintObstacles_);
+
+  // Create a new instance of IpoptApplication
+    Ipopt::SmartPtr<Ipopt::IpoptApplication> app = new Ipopt::IpoptApplication();
+
+  // Change some options
+  app->Options()->SetNumericValue("tol", this->tolerance_);
+  app->Options()->SetIntegerValue("max_iter", this->maxIterations_);
+  app->Options()->SetStringValue("mu_strategy", "adaptive");
+  if (linearSolverType_ != "")
+    app->Options()->SetStringValue("linear_solver",linearSolverType_);
+  app->Options()->SetStringValue("output_file", "ipopt.out");
+  app->Options()->SetStringValue("hessian_constant", "yes");
+
+  // We only support linear constraints
+  app->Options()->SetStringValue("jac_d_constant","yes");
+  app->Options()->SetStringValue("jac_c_constant","yes");
+
+  switch (this->verbosity_) {
+  case NumProc::QUIET:
+      app->Options()->SetIntegerValue("print_level", 0);
+      break;
+  case NumProc::REDUCED:
+      app->Options()->SetIntegerValue("print_level", 5);
+      break;
+  case NumProc::FULL:
+      app->Options()->SetIntegerValue("print_level", 12);
+  }
+
+  // Initialize the IpoptApplication and process the options
+  Ipopt::ApplicationReturnStatus status;
+  status = app->Initialize();
+  if (status != Ipopt::Solve_Succeeded)
+      DUNE_THROW(Dune::Exception, "IPOpt: Error during initialization!");
+
+  // Ask Ipopt to solve the problem
+  status = app->OptimizeTNLP(mynlp);
+
+  // Print helpful text for at least some cases
+  if (status == Ipopt::Invalid_Option)
+    DUNE_THROW(Dune::Exception, "IPOpt returned 'Invalid_Option' error!");
+
+  if (status == Ipopt::Solved_To_Acceptable_Level)
+      std::cout<<"WARNING: Desired tolerance could not be reached, but still acceptable tolerance is reached.\n";
+  else if (status != Ipopt::Solve_Succeeded)
+      DUNE_THROW(Dune::Exception, "IPOpt: Error during optimization!");
+}
+
+#endif

dune/solvers/solvers/linearsolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_LINEARSOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_LINEARSOLVER_HH
+
+#include <dune/solvers/solvers/solver.hh>
+
+namespace Dune {
+
+namespace Solvers {
+
+/** \brief Abstract base class for solvers that solve linear problems
+ *
+ * Linear problems are problems that involve a fixed matrix and a right-hand-side
+ * vector.  Additional constraints are allowed.
+ */
+template <class Matrix, class Vector>
+class LinearSolver : public Solver
+{
+public:
+  /** \brief Constructor taking all relevant data */
+  LinearSolver(VerbosityMode verbosity)
+  : Solver(verbosity)
+  {}
+
+  /** \brief Set up the linear problem to solve */
+  virtual void setProblem(const Matrix& matrix,
+                          Vector& x,
+                          const Vector& rhs) = 0;
+
+};
+
+}   // namespace Solvers
+
+}   // namespace Dune
+
+#endif

dune/solvers/solvers/loopsolver.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+
 #include <cmath>
 #include <limits>
 #include <iostream>
@@ -5,25 +8,13 @@
 #include <dune/solvers/solvers/solver.hh>
 
 template <class VectorType, class BitVectorType>
-void ::LoopSolver<VectorType, BitVectorType>::check() const
-{
-    if (!iterationStep_)
-        DUNE_THROW(SolverError, "You need to supply an iteration step to an iterative solver!");
-
-    iterationStep_->check();
-
-    // check base class
-    IterativeSolver<VectorType,BitVectorType>::check();
-}
-
-template <class VectorType, class BitVectorType>
-void LoopSolver<VectorType, BitVectorType>::preprocess()
+void Dune::Solvers::LoopSolver<VectorType, BitVectorType>::preprocess()
 {
     this->iterationStep_->preprocess();
 }
 
 template <class VectorType, class BitVectorType>
-void LoopSolver<VectorType, BitVectorType>::solve()
+void Dune::Solvers::LoopSolver<VectorType, BitVectorType>::solve()
 {
 
     int i;
@@ -53,7 +44,13 @@ void LoopSolver<VectorType, BitVectorType>::solve()
         std::cout << "     rate";
         std::string header = this->iterationStep_->getOutput();
         std::cout << header;
+
+        for(auto&& c: criteria_)
+            std::cout << c.header();
+
+
         std::cout << std::endl;
+
         std::cout << "-----";
         if (this->useRelativeError_)
             std::cout << "---------------";
@@ -63,6 +60,10 @@ void LoopSolver<VectorType, BitVectorType>::solve()
         std::cout << "---------";
         for(size_t i=0; i<header.size(); ++i)
             std::cout << "-";
+
+        for(auto&& c: criteria_)
+            std::cout << std::string(c.header().size(), '-');
+
         std::cout << std::endl;
     }
 
@@ -77,15 +78,17 @@ void LoopSolver<VectorType, BitVectorType>::solve()
     // Loop until desired tolerance or maximum number of iterations is reached
     for (i=0; i<this->maxIterations_ && (error>this->tolerance_ || std::isnan(error)); i++)
     {
+        iter_ = i;
+
         // Backup of the current solution for the error computation later on
-        VectorType oldSolution = iterationStep_->getSol();
+        VectorType oldSolution = this->iterationStep_->getSol();
 
         // Perform one iteration step
-        iterationStep_->iterate();
+        this->iterationStep_->iterate();
 
         // write iteration to file, if requested
         if (this->historyBuffer_!="")
-            this->writeIterate(iterationStep_->getSol(), i);
+            this->writeIterate(this->iterationStep_->getSol(), i);
 
         // Compute error
         real_type oldNorm = this->errorNorm_->operator()(oldSolution);
@@ -94,64 +97,88 @@ void LoopSolver<VectorType, BitVectorType>::solve()
 
         // Please don't replace this call to 'diff' by computing the norm of the difference.
         // In some nonlinear DD applications the 'diff' method may be nonlinear.
-        real_type normOfCorrection = this->errorNorm_->diff(oldSolution,iterationStep_->getSol());
-        real_type convRate = normOfCorrection / normOfOldCorrection;
+        real_type normOfCorrection = this->errorNorm_->diff(oldSolution,this->iterationStep_->getSol());
+        convRate_ = (normOfOldCorrection > 0)
+            ? normOfCorrection / normOfOldCorrection : 0.0;
         error = normOfCorrection;
         normOfOldCorrection = normOfCorrection;
 
         // If a reference solution has been provided compute the error with respect to it
         if (referenceSolution_)
         {
-            normOfError = this->errorNorm_->diff(iterationStep_->getSol(), *referenceSolution_);
-            convRate = normOfError / normOfOldError;
+            normOfError = this->errorNorm_->diff(this->iterationStep_->getSol(), *referenceSolution_);
+            convRate_ = (normOfOldError > 0) ? normOfError / normOfOldError : 0.0;
             error = normOfError;
             normOfOldError = normOfError;
         }
 
         // Turn the error into the relative error, if requested
-        if (this->useRelativeError_ && !std::isnan(error/oldNorm))
-            error = error / oldNorm;
+        if (this->useRelativeError_ && error != 0)
+            error = (oldNorm == 0) ? std::numeric_limits<real_type>::max()
+                                   : error / oldNorm;
 
-        if (!std::isinf(convRate) && !std::isnan(convRate) && i>0)
+        if (!std::isinf(convRate_) && !std::isnan(convRate_) && i>0)
         {
-            totalConvRate *= convRate;
+            totalConvRate *= convRate_;
             this->maxTotalConvRate_ = std::max(this->maxTotalConvRate_, std::pow(totalConvRate, 1/((real_type)convRateCounter+1)));
             convRateCounter++;
         }
 
         // Output
         if (this->verbosity_ == NumProc::FULL) {
+            std::streamsize const oldPrecision = std::cout.precision();
+            std::ios_base::fmtflags const oldFormatFlags = std::cout.flags();
+
             std::cout << std::setw(5) << i;
 
             if (this->useRelativeError_)
             {
-                std::cout << std::setiosflags(std::ios::scientific);
-                std::cout << std::setw(15) << std::setprecision(7) << error;
-                std::cout << std::resetiosflags(std::ios::scientific);
+                std::cout << std::scientific
+                          << std::setw(15) << std::setprecision(7) << error;
             }
 
             if (referenceSolution_)
             {
-                std::cout << std::setiosflags(std::ios::scientific);
-                std::cout << std::setw(15) << std::setprecision(7) << normOfError;
-                std::cout << std::resetiosflags(std::ios::scientific);
+                std::cout << std::scientific
+                          << std::setw(15) << std::setprecision(7) << normOfError;
             }
 
-            std::cout << std::setiosflags(std::ios::scientific);
-            std::cout << std::setw(15) << std::setprecision(7) << normOfCorrection;
-            std::cout << std::resetiosflags(std::ios::scientific);
+            std::cout << std::scientific
+                      << std::setw(15) << std::setprecision(7) << normOfCorrection;
 
-            std::cout << std::setiosflags(std::ios::fixed);
-            if (i==0)  // We can't estimate the convergence rate at the first iteration
+            if (i == 0)
+                // We can't estimate the convergence rate at the first iteration
                 std::cout << "         ";
             else
-                std::cout << std::setw(9) << std::setprecision(5) << convRate;
-            std::cout << std::resetiosflags(std::ios::fixed);
+                std::cout << std::fixed
+                          << std::setw(9) << std::setprecision(5) << convRate_;
+
+            std::cout << std::setprecision(oldPrecision);
+            std::cout.flags(oldFormatFlags);
 
             std::cout << this->iterationStep_->getOutput();
+        }
+
+        // execute the stop criteria regardless of the verbosity
+        bool stop = false;
+        for(auto&& c: criteria_)
+        {
+            auto r = c();
+            stop = stop or std::get<0>(r);
+            if (this->verbosity_ == NumProc::FULL)
+            {
+              std::cout << std::get<1>(r);
+            }
+        }
+
+        if (this->verbosity_ == NumProc::FULL)
+        {
             std::cout << std::endl;
         }
 
+        if (stop)
+            break;
+
     }
 
 

dune/solvers/solvers/loopsolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
 #ifndef LOOP_SOLVER_HH
 #define LOOP_SOLVER_HH
 
+#include <type_traits>
+
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/iterationsteps/iterationstep.hh>
 #include <dune/solvers/norms/norm.hh>
+#include <dune/solvers/solvers/criterion.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
+
+
+namespace Dune {
+
+namespace Solvers {
 
 /** \brief A solver which consists of a single loop
   *
   *  This class basically implements a loop that calls an iteration procedure
   *  (which is to be supplied by the user).  It also monitors convergence.
   */
-template <class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
+template <class VectorType, class BitVectorType = Dune::Solvers::DefaultBitVector_t<VectorType> >
 class LoopSolver : public IterativeSolver<VectorType, BitVectorType>
 {
     typedef typename VectorType::field_type field_type;
@@ -21,6 +32,7 @@ class LoopSolver : public IterativeSolver<VectorType, BitVectorType>
 public:
 
     /** \brief Constructor taking all relevant data */
+    [[deprecated("Handing over raw pointer in the constructor is deprecated!")]]
     LoopSolver(IterationStep<VectorType, BitVectorType>* iterationStep,
                int maxIterations,
                double tolerance,
@@ -30,17 +42,80 @@ public:
                const VectorType* referenceSolution=0)
         : IterativeSolver<VectorType, BitVectorType>(maxIterations,
                                                      tolerance,
-                                                     errorNorm,
+                                                     *errorNorm,
                                                      verbosity,
                                                      useRelativeError),
-          iterationStep_(iterationStep),
           referenceSolution_(referenceSolution)
-    {}
+    {
+        this->setIterationStep(*iterationStep);
+    }
+
+    /** \brief Constructor taking all relevant data */
+    template <class Step, class RealNorm, class Enable = std::enable_if_t<
+                not std::is_pointer<Step>::value and not std::is_pointer<RealNorm>::value> >
+    LoopSolver(Step&& iterationStep,
+               int maxIterations,
+               double tolerance,
+               RealNorm&& errorNorm,
+               Solver::VerbosityMode verbosity,
+               bool useRelativeError=true,
+               const VectorType* referenceSolution=0)
+        : IterativeSolver<VectorType, BitVectorType>(maxIterations,
+                                                     tolerance,
+                                                     std::forward<RealNorm>(errorNorm),
+                                                     verbosity,
+                                                     useRelativeError),
+          referenceSolution_(referenceSolution)
+    {
+        this->setIterationStep(std::forward<Step>(iterationStep));
+    }
+
+    /**
+     * \brief Add a convergence criterion
+     *
+     * All criteria are checked after each loop. The solver loop will
+     * terminate once any of the supplied criteria evaluates to true.
+     * If you want to terminate only if several criteria are true
+     * you can combine criteria using the overloaded bitwise
+     * logical operations.
+     *
+     * Besides checking the criterion the LoopSolver will
+     * print the diagnostic string returned by the criterion
+     * in a column with the criterions header string on top.
+     *
+     * This method forward all arguments to the constructor of
+     * Criterion. So you can pass a Criterion or any argument
+     * list supported by the Criterion constructors.
+     *
+     * Notice, that the old hard-wired criteria are still active.
+     * If you don't supply any criteria here, the you will get
+     * exactly the old behaviour.
+     *
+     * \attention This is an experimental feature that may still change in the near future.
+     */
+    template<class... Args>
+    void addCriterion(Args&&... args)
+    {
+        criteria_.emplace_back(std::forward<Args>(args)...);
+    }
+
+    /**
+     * \brief Get current iteration count
+     */
+    int iterationCount() const
+    {
+        return iter_;
+    }
+
+
+    /**
+     * \brief Get current convergence rate
+     */
+    auto convergenceRate() const
+    {
+        return convRate_;
+    }
 
-    /**  \brief Checks whether all relevant member variables are set
-      *  \exception SolverError if the iteration step is not set up properly
-      */
-    virtual void check() const;
 
     virtual void preprocess();
 
@@ -49,13 +124,23 @@ public:
       */
     virtual void solve();
 
-    //! The iteration step used by the algorithm
-    IterationStep<VectorType, BitVectorType>* iterationStep_;
-
     const VectorType* referenceSolution_;
+protected:
+    std::vector<Dune::Solvers::Criterion> criteria_;
+
+    int iter_;
 
+    real_type convRate_;
 };
 
+}  // namespace Solvers
+
+}  // namespace Dune
+
+// For backward compatibility: will be removed eventually
+
+using Dune::Solvers::LoopSolver;
+
 #include "loopsolver.cc"
 
 #endif

dune/solvers/solvers/proximalnewtonsolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_PROXIMALNEWTONSOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_PROXIMALNEWTONSOLVER_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
+#include <dune/solvers/common/resize.hh>
+#include <dune/solvers/solvers/criterion.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+
+namespace Dune::Solvers
+{
+  namespace ProximalNewton
+  {
+    /** \brief List of the four stages of the proximal Newton step
+     *
+     *  These are used to select proper regularization rules and are handed over to the
+     *  regularization update method.
+     */
+    enum Stage
+    {
+      minimize,      // first stage: minimizing the second order problem
+      configuration, // second stage: testing the new configuration x + dx
+      descent,       // third stage: testing the descent criteria for the new configuration
+      accepted       // last stage: the new increment was accepted
+    };
+
+    //! A dummy class for g=0 in the ProximalNewtonSolver
+    template<class VectorType>
+    struct ZeroFunctional
+    {
+      double operator()( const VectorType& dx ) const
+      {
+        return 0.0;
+      }
+
+      void updateOffset( const VectorType& x )
+      {
+        // do nothing
+      }
+    };
+
+
+    // A simple regularization updater which doubles in case of failure and halves in case of success
+    struct SimpleRegUpdater
+    {
+      void operator()( double& regWeight, Stage stage) const
+      {
+        if ( stage == Stage::accepted )
+          regWeight *= 0.5;
+        else
+          // make it at least 1.0
+          regWeight = std::max( 1.0, 10.0*regWeight );
+      }
+    };
+  }
+
+
+
+
+
+  /** \brief Generic proximal Newton solver to solve a given minimization problem
+   *
+   *  The proximal Newton solver aims to solve a minimization problem given in the form
+   *      Minimize J(x) := f(x) + g(x)
+   *  where f is a C^2 functional and g is a possibly non-smooth functional.
+   *  During the minimization, a sequence of increments dx as solutions of the second order subproblems
+   *      Minimize 0.5*f''(x)[dx,dx] +  f'(x)[dx] + g(x + dx) + r*||dx||^2
+   *  is computed until the update x := x + dx converges in some sense.
+   *  The user has to provide a suitable regularization strategy to control the regularization weight r,
+   *  and a proper norm ||.|| for the subproblem.
+   */
+  template<class SEA, class NSF, class SOS, class VectorType, class ErrorNorm, class RegUpdater, class BitVectorType = DefaultBitVector_t<VectorType>>
+  class ProximalNewtonSolver : public Solver, public CanIgnore<BitVectorType>
+  {
+  public:
+
+    using SmoothEnergyAssembler = SEA;
+    using NonsmoothFunctional = NSF;
+    using SecondOrderSolver = SOS;
+    using MatrixType = typename SecondOrderSolver::MatrixType;
+
+    void solve() override;
+
+    /** \brief Constructor taking all relevant data
+     *
+     *  \param sea The SmoothEnergyAssembler representing f: It must provide the method
+     *             assembleGradientAndHessian( x, f', f'' ) in order to compute the second order subproblem, and
+     *             the evaluation by computeEnergy( x ) to return f(x)
+     *  \param nsf The NonsmoothFunctional representing g: It must provide the method
+     *             updateOffset( x ) to update the offset in g( x + dx ), and an evaluation operator().
+     *  \param sos The SecondOrderSolver which is able to minimize the second order subproblem. It must provide
+     *             a method minimize( f'', f', g, r, x, ignore) which overwrites the parameter x with the minimizer
+     *             and throws a Dune::Exception in case the minimization failed.
+     *  \param solution This is the solution of the global problem. It is overwritten during the computation and serves
+     *                  also as the initial value.
+     *  \param errorNorm This is the Solvers::EnergyNorm used in the second order problem and also in the descent criteria.
+     *  \param regUpdater The regularization strategy. It must provide a call operator ( r, Stage ) that overwrites r
+     *                    based on the given Stage of the computation
+     *  \param initialRegularizationWeight The initial regularization weight to begin with
+     *  \param maxIterations The maximal number of proximal Newton steps before the Proximal Newton solver aborts the loop
+     *  \param threshold The threshold to stop the iteration once || dx || < threshold
+     *  \param verbosity If the verbosity is set to Solver::FULL the ProximalNewtonSolver will print a table showing
+     *                   the current iterations and some useful information.
+     */
+    ProximalNewtonSolver( const SmoothEnergyAssembler& sea,
+                          NonsmoothFunctional& nsf,
+                          const SecondOrderSolver& sos,
+                          VectorType& solution,
+                          const ErrorNorm& errorNorm,
+                          const RegUpdater& regUpdater,
+                          double& initialRegularizationWeight,
+                          int maxIterations,
+                          double threshold,
+                          Solver::VerbosityMode verbosity)
+    : smoothEnergyAssembler_(&sea)
+    , nonsmoothFunctional_(&nsf)
+    , sos_(&sos)
+    , solution_(&solution)
+    , norm_(&errorNorm)
+    , regUpdater_(regUpdater)
+    , regWeight_(&initialRegularizationWeight)
+    , maxIterations_(maxIterations)
+    , threshold_(threshold)
+    , verbosity_(verbosity)
+    {}
+
+
+
+    /** \brief Add a stop criterion to be executed at the beginning of the loop */
+    template<class... Args>
+    void addStopCriterion(Args&&... args)
+    {
+      stopCriteria_.emplace_back(std::forward<Args>(args)...);
+    }
+
+    /** \brief Add a descent criterion to be executed in the descent stage of the loop */
+    template<class... Args>
+    void addDescentCriterion(Args&&... args)
+    {
+      descentCriteria_.emplace_back(std::forward<Args>(args)...);
+    }
+
+    /** \brief Return the currently computed gradient of smooth energy part */
+    const auto& gradient() const
+    {
+      return *gradientPtr_;
+    }
+
+    /** \brief Return the currently computed hessian of smooth energy part */
+    const auto& hessian() const
+    {
+      return *hessianPtr_;
+    }
+
+    /** \brief Check the existence of a current hessian matrix */
+    bool hasHessian() const
+    {
+      return hessianPtr_ != nullptr;
+    }
+
+    /** \brief Check the existence of a current gradient vector */
+    bool hasGradient() const
+    {
+      return gradientPtr_ != nullptr;
+    }
+
+    /** \brief Set a hessian from the outside */
+    void setHessian( const std::shared_ptr<MatrixType>& hessianPtr )
+    {
+      hessianPtr_ = hessianPtr;
+    }
+
+    /** \brief Set a gradient from the outside */
+    void setGradient( const std::shared_ptr<VectorType>& gradientPtr )
+    {
+      gradientPtr_ = gradientPtr;
+    }
+
+    /** \brief Return the current computed correction in x */
+    const auto& correction() const
+    {
+      return *correction_;
+    }
+
+    /** \brief Access the currently used nonsmooth part (it changes due to the offsets) */
+    const auto& nonsmoothFunctional() const
+    {
+      return *nonsmoothFunctional_;
+    }
+
+    /** \brief direct access to the current regularization weight with read/write possibility */
+    auto& regularizationWeight()
+    {
+      return *regWeight_;
+    }
+
+    /** \brief Get current iteration number */
+    int getIterationCount()
+    {
+      return iter_;
+    }
+
+  private:
+
+    const SmoothEnergyAssembler* smoothEnergyAssembler_;
+    NonsmoothFunctional* nonsmoothFunctional_;
+    const SecondOrderSolver* sos_;
+
+    // current iterate of the solution of the minimization problem
+    VectorType* solution_;
+
+    // increments of the proximal Newton step
+    std::shared_ptr<VectorType> correction_;
+
+    const ErrorNorm* norm_;
+
+    const RegUpdater regUpdater_;
+
+    // current regularization weight
+    double* regWeight_;
+
+    int iter_;
+    int maxIterations_;
+    double threshold_;
+    Solver::VerbosityMode verbosity_;
+
+    // store the different criteria
+    std::vector<Dune::Solvers::Criterion> stopCriteria_;
+    std::vector<Dune::Solvers::Criterion> descentCriteria_;
+
+    // access to the internal data for external criteria
+    std::shared_ptr<MatrixType> hessianPtr_ = nullptr;
+    std::shared_ptr<VectorType> gradientPtr_ = nullptr;
+
+    void printLine( int iter, double usedReg, double normCorrection, double newEnergy, double energyDiff, std::string errorMessage = "") const
+    {
+      std::cout << std::setw( 7) << iter << "  |  "
+      << std::setw(15) << std::setprecision(9) << usedReg << " | "
+      << std::setw(15) << std::setprecision(9) << normCorrection << " | "
+      << std::setw(15) << std::setprecision(9) << newEnergy << " | "
+      << std::setw(15) << std::setprecision(9) << energyDiff << "   "
+      << errorMessage << std::endl;
+    };
+  };
+
+
+
+  template<class SEA, class NSF, class SOS, class V, class EN, class RU, class BV>
+  void ProximalNewtonSolver<SEA,NSF,SOS,V,EN,RU,BV>::solve()
+  {
+    using VectorType = V;
+
+    const bool printOutput = this->verbosity_ == NumProc::FULL;
+
+    auto& regWeight = *regWeight_;
+
+    if ( printOutput )
+    {
+      std::cout << " iterate |   regularization |      correction |          energy | energy difference "<< std::endl;
+      std::cout << "---------+------------------+-----------------+-----------------+-------------------"<< std::endl;
+    }
+
+    iter_ = 0;
+
+    if ( not correction_ )
+      correction_ = std::make_shared<VectorType>();
+
+    // we need a zero vector for computing concrete energy descents later
+    VectorType zeroVector;
+    resizeInitializeZero(*correction_, *solution_);
+    resizeInitializeZero(zeroVector, *solution_);
+
+    using real_type = typename VectorType::field_type;
+    real_type normCorrection = std::numeric_limits<double>::max();
+
+    // start the loop
+    for( iter_ = 0; iter_ < this->maxIterations_; iter_++ )
+    {
+      // check for ||dx|| < threshold
+      if ( (1.0 + regWeight)*normCorrection < threshold_ )
+      {
+        if ( printOutput )
+          std::cout << "ProximalNewtonSolver terminated because of weighted correction is below threshold: " << (1.0 + regWeight)*normCorrection << std::endl;
+        break;
+      }
+
+      // check user added additional stop criteria
+      bool stop = false;
+      for ( auto&& c: stopCriteria_ )
+      {
+        auto r = c();
+        if ( std::get<0>(r) )
+        {
+          if ( printOutput )
+            std::cout << "ProximalNewtonSolver terminated because of a user added stop criterion: " << std::get<1>( r ) << std::endl;
+
+          stop = true;
+          break;
+        }
+      }
+
+      // don't do another loop
+      if ( stop )
+        break;
+
+      // keep a copy
+      auto usedReg = *regWeight_;
+
+
+      // store some information in case the step gets discarded
+      auto oldX = *solution_;
+
+      // assemble the quadratic and linear part if not recycled from previous step
+      if ( not hasGradient() or not hasHessian() )
+      {
+        hessianPtr_ = std::make_shared<MatrixType>();
+        gradientPtr_ = std::make_shared<VectorType>();
+
+        smoothEnergyAssembler_->assembleGradientAndHessian( oldX, *gradientPtr_, *hessianPtr_ );
+      }
+
+      // shift the nonsmoothFunctional by the current x
+      nonsmoothFunctional_->updateOffset( oldX );
+
+      *correction_ = 0.0;
+
+      // remember the old energy: note that nonsmoothFunctional_ is already shifted by oldX
+      auto oldEnergy = smoothEnergyAssembler_->computeEnergy( oldX ) + (*nonsmoothFunctional_)( zeroVector );
+
+      ///////////////////////////////////////////////////////////////////////////////////////////
+      /// Stage I: Try to compute a Proximal Newton step ////////////////////////////////////////
+      ///////////////////////////////////////////////////////////////////////////////////////////
+
+      // compute one Proximal Newton Step with the second order solver
+      try
+      {
+        sos_->minimize( *hessianPtr_, *gradientPtr_, *nonsmoothFunctional_, regWeight, *correction_, this->ignore() );
+      }
+      catch(const MathError& e)
+      {
+        if ( printOutput )
+          printLine( iter_, usedReg, 0, oldEnergy, 0, "The Proximal Newton Step reported an error: " + std::string(e.what()) );
+
+        regUpdater_(regWeight, ProximalNewton::Stage::minimize );
+        continue;
+      }
+
+
+      ///////////////////////////////////////////////////////////////////////////////////////////
+      /// Stage II: Check that new x is a valid configuration ///////////////////////////////////
+      ///////////////////////////////////////////////////////////////////////////////////////////
+
+
+      // if we got here the correction can be used to compute the next x
+      auto newX = oldX;
+      newX += *correction_;
+
+      // compute the newEnergy: check for invalid configuration
+      real_type newEnergy;
+      try
+      {
+        newEnergy = smoothEnergyAssembler_->computeEnergy( newX ) + (*nonsmoothFunctional_)( *correction_ );
+      }
+      catch(const MathError& e)
+      {
+        if ( printOutput )
+          printLine( iter_, usedReg, 0, oldEnergy, 0, "Computing the new energy resulted in an error: " + std::string(e.what()) );
+
+        regUpdater_(regWeight, ProximalNewton::Stage::configuration );
+        continue;
+      }
+
+      // Compute objective function descent
+      auto energyDiff = newEnergy;
+      energyDiff -= oldEnergy;
+
+
+      ///////////////////////////////////////////////////////////////////////////////////////////
+      /// Stage III: Check that the new x fulfills descent criteria /////////////////////////////
+      ///////////////////////////////////////////////////////////////////////////////////////////
+
+      // check user added additional descent criteria
+      bool accepted = true;
+      std::string errorMessage;
+      for ( auto&& c: descentCriteria_ )
+      {
+        auto r = c();
+        if ( not std::get<0>( r ) )
+        {
+          if ( printOutput )
+            errorMessage = std::get<1>( r );
+
+          accepted = false;
+          break;
+        }
+      }
+
+      if ( not accepted )
+      {
+        if ( printOutput )
+          printLine( iter_, usedReg, 0, oldEnergy, 0, "The following descent criterion was not accepted: " + errorMessage );
+
+        regUpdater_(regWeight, ProximalNewton::Stage::descent );
+        continue;
+      }
+
+
+      normCorrection = (*norm_)( *correction_ );
+
+      if ( printOutput )
+      {
+        printLine( iter_, usedReg, normCorrection, newEnergy, energyDiff );
+      }
+
+
+      ///////////////////////////////////////////////////////////////////////////////////////////
+      /// Stage IV: Update the regularization weight for the next step  /////////////////////////
+      ///////////////////////////////////////////////////////////////////////////////////////////
+
+      regUpdater_(regWeight, ProximalNewton::Stage::accepted );
+
+
+      // seems like the step was accepted:
+      *solution_ = newX;
+
+      // reset gradient and hessian since x is updated
+      gradientPtr_.reset();
+      hessianPtr_.reset();
+    }
+  }
+} // namespace Dune::Solvers
+
+
+#endif