dune/solvers/iterationsteps/mmgstep.hh

@@ -54,7 +54,7 @@ public:
     virtual void nestedIteration();
 
     //! Set the hasObstacle bitfield
-    DUNE_DEPRECATED_MSG("Setting by raw pointer is deprecated. Use l-value or r-value or a shared_ptr")
+    [[deprecated("Setting by raw pointer is deprecated. Use l-value or r-value or a shared_ptr")]]
     void setHasObstacles(Dune::BitSetVector<dim>* hasObstacle) {
       hasObstacle_ = Dune::stackobject_to_shared_ptr(*hasObstacle);
     }
@@ -66,7 +66,7 @@ public:
     }
 
     //! Set the obstacle field
-    DUNE_DEPRECATED_MSG("Setting by raw pointer is deprecated. Use l-value or r-value or a shared_ptr")
+    [[deprecated("Setting by raw pointer is deprecated. Use l-value or r-value or a shared_ptr")]]
     void setObstacles(ObstacleVectorType* obstacles) {
       obstacles_ = Dune::stackobject_to_shared_ptr(*obstacles);
     }

dune/solvers/iterationsteps/multigridstep.cc

@@ -73,17 +73,14 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
 
     xHierarchy_.resize(numLevels());
     rhsHierarchy_.resize(numLevels());
-    matrixHierarchy_.resize(numLevels());
     ignoreNodesHierarchy_.resize(numLevels());
 
-    for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+    for (int i=0; i<int(xHierarchy_.size())-1; i++)
     {
-        matrixHierarchy_[i].reset();
         xHierarchy_[i].reset();
         ignoreNodesHierarchy_[i] = NULL;
     }
 
-    matrixHierarchy_.back() = this->mat_;
     xHierarchy_.back()      = Dune::stackobject_to_shared_ptr(*(this->x_));
     rhsHierarchy_.back()    = *(this->rhs_);
 
@@ -117,22 +114,35 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
     // /////////////////////////////////////////////////////
     //   Assemble the complete hierarchy of matrices
     // /////////////////////////////////////////////////////
-    for (int i=this->numLevels()-2; i>=0; i--)
+    if (not matrixHierarchyManuallySet_)
     {
-        this->matrixHierarchy_[i] = std::make_shared<MatrixType>();
-        this->xHierarchy_[i] = std::make_shared<VectorType>();
+        matrixHierarchy_.resize(numLevels());
+        for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+            matrixHierarchy_[i].reset();
+        matrixHierarchy_.back() = this->mat_;
+        for (int i=this->numLevels()-2; i>=0; i--)
+        {
+            this->matrixHierarchy_[i] = std::make_shared<MatrixType>();
 
-        // Compute which entries are present in the (sparse) coarse grid stiffness
-        // matrices.
-        this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+            // Compute which entries are present in the (sparse) coarse grid stiffness
+            // matrices.
+            this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+
+            // setup matrix
+            this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+        }
+        matrixHierarchyManuallySet_ = false;
+    }
 
-        // setup matrix
-        this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+    for (int i=this->numLevels()-2; i>=0; i--)
+    {
+        this->xHierarchy_[i] = std::make_shared<VectorType>();
 
         // Set solution vector sizes for the lower levels
         MatrixVector::resize(*(this->xHierarchy_[i]),*this->matrixHierarchy_[i]);
     }
 
+
     // /////////////////////////////////////////////////////
     //   Setup dirichlet bitfields
     // /////////////////////////////////////////////////////

dune/solvers/iterationsteps/multigridstep.hh

@@ -24,8 +24,6 @@ namespace Dune {
     class MultigridStep : public LinearIterationStep<MatrixType, VectorType, BitVectorType>
     {
 
-        static const int blocksize = VectorType::block_type::dimension;
-
     public:
         using LinearStepType = LinearIterationStep<MatrixType, VectorType, BitVectorType>;
 
@@ -95,7 +93,7 @@ namespace Dune {
         }
 
         template <class DerivedTransfer>
-        DUNE_DEPRECATED_MSG("Consider setting the transfer operators via smart pointers instead.")
+        [[deprecated("Consider setting the transfer operators via smart pointers instead.")]]
         void setTransferOperators(const std::vector<DerivedTransfer*>& transfer)
         {
             mgTransfer_.resize(transfer.size());
@@ -153,7 +151,7 @@ namespace Dune {
         virtual void setMGType(int mu, int nu1, int nu2);
 
         /** \brief Set the smoother iteration step */
-        DUNE_DEPRECATED_MSG("Consider setting the smoother via smart pointer, reference or temporaries instead.")
+        [[deprecated("Consider setting the smoother via smart pointer, reference or temporaries instead.")]]
         void setSmoother(LinearStepType* smoother)
         {
             presmootherDefault_ = postsmootherDefault_ = Dune::stackobject_to_shared_ptr(*smoother);
@@ -171,7 +169,7 @@ namespace Dune {
         }
 
         /** \brief Set pre- and post smoothers individually */
-        DUNE_DEPRECATED_MSG("Consider setting the smoother via smart pointer, reference or temporaries instead.")
+        [[deprecated("Consider setting the smoother via smart pointer, reference or temporaries instead.")]]
         void setSmoother(LinearStepType* preSmoother,
                          LinearStepType* postSmoother)
         {
@@ -190,7 +188,7 @@ namespace Dune {
         }
 
         /** \brief Set the smoother iteration step for a particular level */
-        DUNE_DEPRECATED_MSG("Consider setting the smoother via smart pointer, reference or temporaries instead.")
+        [[deprecated("Consider setting the smoother via smart pointer, reference or temporaries instead.")]]
         void setSmoother(LinearStepType* smoother, std::size_t level)
         {
             levelWiseSmoothers_[level] = Dune::stackobject_to_shared_ptr(*smoother);
@@ -210,6 +208,24 @@ namespace Dune {
             basesolver_ = wrap_own_share<Solver>(std::forward<BaseSolver>(baseSolver));
         }
 
+        /**
+         * \brief Set pre-coarsened matrix hierarchy
+         *
+         * This allows to set the coarsened matrix hierarchy manually.
+         * It is assumed, that the finest matrix coincides with
+         * the one passed to setProblem(). The latter is ignored.
+         * This has to be called before preprocess().
+         */
+        template<class M>
+        void setMatrixHirarchy(const std::vector<std::shared_ptr<M>>& matrixHierarchy)
+        {
+            matrixHierarchy_.resize(matrixHierarchy.size());
+            for(auto i : Dune::range(matrixHierarchy.size()))
+                matrixHierarchy_[i] = matrixHierarchy[i];
+            matrixHierarchyManuallySet_ = true;
+        }
+
+
     protected:
         /** \brief The presmoothers, one for each level */
         std::vector<std::shared_ptr<LinearStepType> > presmoother_;
@@ -228,6 +244,8 @@ namespace Dune {
 
         //! Number of coarse corrections
         int mu_;
+
+        bool matrixHierarchyManuallySet_= false;
     public:
         //! Variable used to store the current level
         int level_;

dune/solvers/iterationsteps/obstacletnnmgstep.hh

@@ -40,7 +40,7 @@
  * details please refer to
  *
  *   'Multigrid Methods for Obstacle Problems', C. Gräser and R. Kornhuber
- *   Journal of Computational Mathematics 27 (1), 2009, pp. 1-44 
+ *   Journal of Computational Mathematics 27 (1), 2009, pp. 1-44
  *
  * A much more general and flexible implementation that can also be used for
  * energy functionals with nonquadratic (anisotropic) smooth part
@@ -241,7 +241,7 @@ class ObstacleTNNMGStep
                         truncatedResidual_[i][ii] = 0;
             }
 
-            // apply linear multigrid to approximatively solve the truncated linear system
+            // apply linear multigrid to approximately solve the truncated linear system
             linearMGStep_.setProblem(truncatedMat, coarseCorrection_, truncatedResidual_);
             linearMGStep_.ignoreNodes_ = ignoreNodes_;
             linearMGStep_.preprocess();
@@ -345,13 +345,13 @@ class ObstacleTNNMGStep
          * is ignored if it is associated to an ignored fine dof
          * in the sense that the corresponding transfer operators entry is 1.
          *
-         * On each level a fixed number of v-cycles is performed.
+         * On each level a fixed number of V-cycles is performed.
          *
          * This method can be called before the preprocess() method.
          * It does not change the ObstacleTNNMGStep state despite
          * updating the solution vector.
          *
-         * \param coarseIterationSteps Number of v-cycle performed on the corse levels.
+         * \param coarseIterationSteps Number of V-cycles performed on the coarse levels.
          */
         void nestedIteration(int coarseIterationSteps=2)
         {

dune/solvers/iterationsteps/projectedlinegsstep.cc

@@ -178,7 +178,7 @@ void ProjectedLineGSStep<MatrixType, VectorType, BitVectorType>::iterate()
 
         const int current_block_size = this->blockStructure_[b_num].size();
 
-        //! compute and save the residuals for the curent block:
+        //! compute and save the residuals for the current block:
         // determine the (permuted) residuals r[p(i)],..., r[p(i+current_block_size-1)]
         // this means that we determine the residuals for the current block
         VectorType permuted_r_i(current_block_size);

dune/solvers/norms/energynorm.hh

@@ -60,6 +60,11 @@ namespace Solvers {
             matrixProvider_ = [=]() -> const MatrixType& { return *matrix; };
         }
 
+        //! \brief get the energy norm matrix
+        const MatrixType& getMatrix() const {
+            return matrixProvider_();
+        }
+
         //! \brief sets to use the current problem matrix of the linear iteration step
         template<class BV>
         void setIterationStep(LinearIterationStep<MatrixType, VectorType, BV>* step) {
@@ -89,7 +94,7 @@ namespace Solvers {
         }
 
         // \brief Compute the squared norm for a given vector and matrix
-        DUNE_DEPRECATED static field_type normSquared(const VectorType& u,
+        [[deprecated]] static field_type normSquared(const VectorType& u,
                                                       const MatrixType& A,
                                                       const field_type tol=1e-10)
         {

dune/solvers/solvers/CMakeLists.txt

@@ -6,6 +6,7 @@ install(FILES
     linearsolver.hh
     loopsolver.cc
     loopsolver.hh
+    proximalnewtonsolver.hh
     quadraticipopt.hh
     solver.hh
     tcgsolver.cc

dune/solvers/solvers/cholmodsolver.hh

@@ -46,7 +46,14 @@ 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,
@@ -57,56 +64,117 @@ public:
     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
    *
-   * Warning! 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.
+   * 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;
-    x_ = &x;
+    isFactorized_ = false;
+  }
+
+  /** \brief Set the rhs vector for the linear problem
+   */
+  void setRhs(const VectorType& rhs)
+  {
     rhs_ = &rhs;
   }
 
-  virtual void solve()
+  /** \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;
-    Dune::Cholmod<VectorType> solver;
 
-    // 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
-    solver.cholmodCommonObject().metis_dswitch = 0.0;
+    // 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_;
 
-      // setMatrix will decompose the matrix
-      solver.setMatrix(*matrix_);
-      // get the error code from Cholmod in case something happened
-      errorCode_ = solver.cholmodCommonObject().status;
-
-      solver.apply(*x_, mutableRhs, statistics);
+      // The back-substitution
+      istlCholmodSolver_.apply(*x_, mutableRhs, statistics);
     }
     else
     {
       const auto& ignore = this->ignore();
 
-      // The setMatrix method stores only the selected entries of the matrix (determined by the ignore field)
-      solver.setMatrix(*matrix_,&ignore);
-      // get the error code from Cholmod in case something happened
-      errorCode_ = solver.cholmodCommonObject().status;
-
       // create flat vectors
       using T = typename VectorType::field_type;
       std::size_t flatSize = flatVectorForEach(ignore, [](...){});
@@ -148,10 +216,15 @@ public:
       });
 
       // Solve the modified system
-      solver.apply(*x_, modifiedRhs, statistics);
+      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
@@ -165,6 +238,9 @@ public:
     return errorCode_;
   }
 
+  //! dune-istl Cholmod solver object
+  Dune::Cholmod<VectorType> istlCholmodSolver_;
+
   //! Pointer to the system matrix
   const MatrixType* matrix_;
 
@@ -176,6 +252,13 @@ public:
 
   //! 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

dune/solvers/solvers/criterion.hh

@@ -52,7 +52,7 @@ namespace Dune {
        * of the header string coincide to get a readable log.
        */
       template<class F,
-        std::enable_if_t<std::is_convertible<std::result_of_t<F()>, Result>::value, int> = 0>
+        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)
@@ -76,7 +76,7 @@ namespace Dune {
        * of the header string coincide to get a readable log.
        */
       template<class F,
-        std::enable_if_t<std::is_convertible<std::result_of_t<F()>, std::string>::value, int> = 0>
+        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)

dune/solvers/solvers/iterativesolver.hh

@@ -40,7 +40,7 @@ namespace Dune {
         {}
 
         /** \brief Constructor taking all relevant data */
-        DUNE_DEPRECATED_MSG("Handing over raw pointer in the constructor is deprecated!")
+        [[deprecated("Handing over raw pointer in the constructor is deprecated!")]]
         IterativeSolver(int maxIterations,
                         double tolerance,
                         const Norm<VectorType>* errorNorm,

dune/solvers/solvers/linearipopt.hh

@@ -584,7 +584,7 @@ void LinearIPOptSolver<VectorType,JacobianType>::solve()
       app->Options()->SetIntegerValue("print_level", 12);
   }
 
-  // Intialize the IpoptApplication and process the options
+  // Initialize the IpoptApplication and process the options
   Ipopt::ApplicationReturnStatus status;
   status = app->Initialize();
   if (status != Ipopt::Solve_Succeeded)
@@ -598,7 +598,7 @@ void LinearIPOptSolver<VectorType,JacobianType>::solve()
     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 accetable tolerance is reached.\n";
+      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!");
 }

dune/solvers/solvers/loopsolver.hh

@@ -32,7 +32,7 @@ class LoopSolver : public IterativeSolver<VectorType, BitVectorType>
 public:
 
     /** \brief Constructor taking all relevant data */
-    DUNE_DEPRECATED_MSG("Handing over raw pointer in the constructor is deprecated!")
+    [[deprecated("Handing over raw pointer in the constructor is deprecated!")]]
     LoopSolver(IterationStep<VectorType, BitVectorType>* iterationStep,
                int maxIterations,
                double tolerance,

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

dune/solvers/solvers/quadraticipopt.hh

@@ -797,7 +797,7 @@ void QuadraticIPOptSolver<MatrixType,VectorType,JacobianType>::solve()
       app->Options()->SetIntegerValue("print_level", 12);
   }
 
-  // Intialize the IpoptApplication and process the options
+  // Initialize the IpoptApplication and process the options
   Ipopt::ApplicationReturnStatus status;
   status = app->Initialize();
   if (status != Ipopt::Solve_Succeeded)
@@ -813,10 +813,24 @@ void QuadraticIPOptSolver<MatrixType,VectorType,JacobianType>::solve()
   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::Search_Direction_Becomes_Too_Small) {
-      std::array<Ipopt::Number,4> inf;
-      app->Statistics()->Infeasibilities(inf[0],inf[1],inf[2],inf[3]);
-      if (inf[3]>std::max(1e-10,this->tolerance_))
-              DUNE_THROW(Dune::Exception,Dune::formatString("Problem could not be solved to acceptable accuracy %d",inf[3]));
+      Ipopt::Number dual_inf;        // dual infeasibility (Gradient of Lagrangian not zero)
+      Ipopt::Number constr_viol;     // violation of constraints
+#if IPOPT_VERSION_MAJOR>=3 && IPOPT_VERSION_MINOR>=14
+      Ipopt::Number varbounds_viol;  // violation of variable bounds
+#endif
+      Ipopt::Number complementarity; // violation of complementarity
+      Ipopt::Number kkt_error;       // KKT error
+
+      app->Statistics()->Infeasibilities(dual_inf,
+                                         constr_viol,
+#if IPOPT_VERSION_MAJOR>=3 && IPOPT_VERSION_MINOR>=14
+                                         varbounds_viol,
+#endif
+                                         complementarity,
+                                         kkt_error);
+
+      if (kkt_error>std::max(1e-10,this->tolerance_))
+              DUNE_THROW(Dune::Exception,Dune::formatString("Problem could not be solved to acceptable accuracy %d", kkt_error));
 
   } else if (status != Ipopt::Solve_Succeeded)
       DUNE_THROW(Dune::Exception, "IPOpt: Error during optimization!");

dune/solvers/solvers/trustregionsolver.hh

@@ -128,10 +128,10 @@ protected:
     /** \brief The maximal trust-region radius in the maximum-norm */
     field_type maxRadius_;
 
-    /** \brief If the iteration is very successfull we enlarge the radius by this factor.*/
+    /** \brief If the iteration is very successful we enlarge the radius by this factor.*/
     field_type enlargeRadius_;
 
-    /** \brief If the iteration is unsuccessfull we shrink the radius by this factor.*/
+    /** \brief If the iteration is unsuccessful we shrink the radius by this factor.*/
     field_type shrinkRadius_;
 
     /** \brief If the ratio between predicted and achieved decrease is above this number, the iteration is very successful.*/

dune/solvers/solvers/umfpacksolver.hh

@@ -12,7 +12,9 @@
 #include <dune/common/exceptions.hh>
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/version.hh>
 
+#include <dune/istl/foreach.hh>
 #include <dune/istl/solver.hh>
 #include <dune/istl/umfpack.hh>
 #include <dune/istl/io.hh>
@@ -61,8 +63,6 @@ public:
 
   void solve() override
   {
-    // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
-    VectorType mutableRhs = *rhs_;
 
     if (not this->hasIgnore())
     {
@@ -72,11 +72,15 @@ public:
       Dune::InverseOperatorResult statistics;
       Dune::UMFPack<MatrixType> solver(*matrix_);
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
+
+      // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
+      VectorType mutableRhs = *rhs_;
       solver.apply(*x_, mutableRhs, statistics);
 
     }
     else
     {
+#if DUNE_VERSION_LTE(DUNE_ISTL, 2, 9)
       ///////////////////////////////////////////////////////////////////////////////////////////
       //  Extract the set of matrix rows that do not correspond to ignored degrees of freedom.
       //  Unfortunately, not all cases are handled by the ISTL UMFPack solver.  Currently,
@@ -99,6 +103,7 @@ public:
             DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
         }
       }
+#endif
 
       // Construct the solver
       Dune::InverseOperatorResult statistics;
@@ -106,49 +111,75 @@ public:
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
       // We eliminate all rows and columns(!) from the matrix that correspond to ignored degrees of freedom.
       // Here is where the sparse LU decomposition is happenening.
+#if DUNE_VERSION_LTE(DUNE_ISTL, 2, 9)
       solver.setSubMatrix(*matrix_,nonIgnoreRows);
+#else
+      solver.setMatrix(*matrix_,this->ignore());
+#endif
+
+      // total number of dofs
+      auto N = flatVectorForEach(*rhs_, [](auto&&, auto&&){});
 
       // We need to modify the rhs vector by static condensation.
-      VectorType shortRhs(nonIgnoreRows.size());
-      int shortRowCount=0;
-      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
-        shortRhs[shortRowCount] = mutableRhs[*it];
+      std::vector<bool> flatIgnore(N);
+      size_t numberOfIgnoredDofs = 0;
 
-      shortRowCount = 0;
-      for (size_t i=0; i<matrix_->N(); i++)
+      flatVectorForEach(this->ignore(), [&](auto&& entry, auto&& offset)
       {
-        if constexpr (std::is_convertible<decltype(this->ignore()[i]), const bool>::value) {
-          if (this->ignore()[i])
-            continue;
-        } else {
-          if (this->ignore()[i].all())
-            continue;
+        flatIgnore[offset] = entry;
+        if ( entry )
+        {
+          numberOfIgnoredDofs++;
         }
+      });
 
-        auto cIt    = (*matrix_)[i].begin();
-        auto cEndIt = (*matrix_)[i].end();
 
-        for (; cIt!=cEndIt; ++cIt)
-          if constexpr (std::is_convertible<decltype(this->ignore()[cIt.index()]), const bool>::value) {
-            if (this->ignore()[cIt.index()])
-              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
-          } else {
-            if (this->ignore()[cIt.index()].all())
-              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
-          }
+      using field_type = typename MatrixType::field_type;
+      std::vector<field_type> shortRhs(N-numberOfIgnoredDofs);
 
-        shortRowCount++;
-      }
+      // mapping of long indices to short indices
+      std::vector<size_t> subIndices(N,std::numeric_limits<size_t>::max());
+
+      int shortRhsCount=0;
+      flatVectorForEach(*rhs_, [&](auto&& entry, auto&& offset)
+      {
+        if ( not flatIgnore[offset] )
+        {
+          shortRhs[shortRhsCount] = entry;
+          subIndices[offset] = shortRhsCount;
+          shortRhsCount++;
+        }
+      });
+
+      std::vector<field_type> flatX(N);
+      flatVectorForEach(*x_, [&](auto&& entry, auto&& offset)
+      {
+        flatX[offset] = entry;
+      });
+
+
+      flatMatrixForEach(*matrix_, [&](auto&& entry, auto&& row, auto&& col)
+      {
+        if ( flatIgnore[col] and not flatIgnore[row] )
+        {
+          shortRhs[ subIndices[ row ] ] -= entry * flatX[ col ];
+        }
+      });
 
       // Solve the reduced system
-      VectorType shortX(nonIgnoreRows.size());
-      solver.apply(shortX, shortRhs, statistics);
+      std::vector<field_type> shortX(N-numberOfIgnoredDofs);
 
-      // Blow up the solution vector
-      shortRowCount=0;
-      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
-        (*x_)[*it] = shortX[shortRowCount];
+      // call the raw-pointer variant of the ISTL UMPACK solver
+      solver.apply(shortX.data(), shortRhs.data());
 
+      // Blow up the solution vector
+      flatVectorForEach(*x_, [&](auto&& entry, auto&& offset)
+      {
+        if ( not flatIgnore[ offset ] )
+        {
+          entry = shortX[ subIndices[ offset ] ];
+        }
+      });
     }
   }
 

dune/solvers/test/CMakeLists.txt

@@ -26,4 +26,6 @@ endif()
 if(SuiteSparse_CHOLMOD_FOUND)
   dune_add_test(SOURCES cholmodsolvertest.cc)
   add_dune_suitesparse_flags(cholmodsolvertest)
+  dune_add_test(SOURCES proximalnewtonsolvertest.cc)
+  add_dune_suitesparse_flags(proximalnewtonsolvertest)
 endif()

dune/solvers/test/cholmodsolvertest.cc

@@ -59,7 +59,80 @@ struct CHOLMODSolverTestSuite
   }
 };
 
+// Check whether we can reuse a factorization for several right hand sides
+bool checkMultipleSolves()
+{
+  bool passed = true;
+
+  using Matrix = BCRSMatrix<double>;
+  using Vector = BlockVector<double>;
+
+  // Construct an example matrix
+  Matrix matrix(4, 4,
+                3,    // Expected average number of nonzeros
+                      // per row
+                0.2,  // Size of the compression buffer,
+                      // as a fraction of the expected total number of nonzeros
+                BCRSMatrix<double>::implicit);
+
+  matrix.entry(0,0) = 2;
+  matrix.entry(0,1) = 1;
+  matrix.entry(1,0) = 1;
+  matrix.entry(1,1) = 2;
+  matrix.entry(1,2) = 1;
+  matrix.entry(2,1) = 1;
+  matrix.entry(2,2) = 2;
+  matrix.entry(2,3) = 1;
+  matrix.entry(3,2) = 1;
+  matrix.entry(3,3) = 2;
+
+  matrix.compress();
+
+  // Set up the solver
+  Solvers::CholmodSolver<Matrix,Vector> cholmodSolver;
+
+  // Where to write the solution
+  Vector x(4);
+  cholmodSolver.setSolutionVector(x);
+
+  // Factorize the matrix
+  cholmodSolver.setMatrix(matrix);
+  cholmodSolver.factorize();
+
+  // Test with a particular right-hand side
+  Vector rhs = {3,4,4,3};
+
+  // The corresponding solution
+  Vector reference = {1,1,1,1};
+
+  cholmodSolver.setRhs(rhs);
+  cholmodSolver.solve();
+
+  // Check whether result is correct
+  auto diff = reference;
+  diff -= x;
+  if (diff.infinity_norm() > 1e-8)
+  {
+    std::cerr << "CholmodSolver didn't produce the correct result!" << std::endl;
+    passed = false;
+  }
 
+  // Test with a different right-hand side
+  // This reuses the original factorization.
+  rhs = {4,6,6,5};
+  reference = {1,2,1,2};
+  cholmodSolver.solve();
+
+  diff = reference;
+  diff -= x;
+  if (diff.infinity_norm() > 1e-8)
+  {
+    std::cerr << "CholmodSolver didn't produce the correct result!" << std::endl;
+    passed = false;
+  }
+
+  return passed;
+}
 
 int main(int argc, char** argv)
 {
@@ -69,8 +142,10 @@ int main(int argc, char** argv)
 
     CHOLMODSolverTestSuite<1> testsuite1;
     CHOLMODSolverTestSuite<2> testsuite2;
-    passed = checkWithStandardGrids(testsuite1);
-    passed = checkWithStandardGrids(testsuite2);
+    passed &= checkWithStandardGrids(testsuite1);
+    passed &= checkWithStandardGrids(testsuite2);
+
+    passed &= checkMultipleSolves();
 
     return passed ? 0 : 1;
 }

dune/solvers/test/common.hh

@@ -12,7 +12,7 @@
 
 #include <dune/geometry/quadraturerules.hh>
 
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
 #include <dune/grid/yaspgrid.hh>
 
@@ -38,7 +38,7 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
 {
     static const int dim = GridView::Grid::dimension;
 
-    typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+    typedef typename Dune::LagrangeLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
     typedef typename FiniteElementCache::FiniteElementType FiniteElement;
 
     const auto& indexSet = gridView.indexSet();
@@ -73,7 +73,7 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
     static constexpr int dim = GridView::Grid::dimension;
     static constexpr int dimworld = GridView::Grid::dimensionworld;
 
-    typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+    typedef typename Dune::LagrangeLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
     typedef typename FiniteElementCache::FiniteElementType FiniteElement;
 
     typedef typename Dune::FieldVector<double, dimworld> GlobalCoordinate;
@@ -89,7 +89,7 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
 
         int localSize = fe.size();
 
-        // get quadrature rule of appropiate order (P1/Q1)
+        // get quadrature rule of appropriate order (P1/Q1)
         int order = (element.type().isSimplex())
                     ? 2*(1-1)
                     : 2*(dim-1);
@@ -143,7 +143,7 @@ void assemblePQ1Mass(const GridView& gridView, Matrix& matrix)
 {
     static const int dim = GridView::Grid::dimension;
 
-    typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+    typedef typename Dune::LagrangeLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
     typedef typename FiniteElementCache::FiniteElementType FiniteElement;
 
     typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
@@ -205,7 +205,7 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
 {
     static const int dim = GridView::Grid::dimension;
 
-    typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+    typedef typename Dune::LagrangeLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
     typedef typename FiniteElementCache::FiniteElementType FiniteElement;
 
     typedef typename FiniteElement::Traits::LocalBasisType::Traits::RangeType RangeType;
@@ -220,7 +220,7 @@ void assemblePQ1RHS(const GridView& gridView, Vector& r, const Function& f)
 
         int localSize = fe.size();
 
-        // get quadrature rule of appropiate order (P1/Q1)
+        // get quadrature rule of appropriate order (P1/Q1)
         int order = (element.type().isSimplex())
                     ? 2*1
                     : 2*dim;
@@ -283,7 +283,7 @@ void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
     static const int dim = GridView::Grid::dimension;
 
     typedef typename GridView::IndexSet IndexSet;
-    typedef typename Dune::PQkLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
+    typedef typename Dune::LagrangeLocalFiniteElementCache<double, double, dim, 1> FiniteElementCache;
     typedef typename FiniteElementCache::FiniteElementType FiniteElement;
 
     const IndexSet& indexSet = gridView.indexSet();
@@ -437,7 +437,7 @@ bool checkWithStandardGrids(TestSuite& suite)
     passed = passed and checkWithGrid<Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double,2> > >(suite, 6);
     passed = passed and checkWithGrid<Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double,3> > >(suite, 4);
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
     passed = passed and checkWithGrid<Dune::UGGrid<2> >(suite, 5);
     passed = passed and checkWithGrid<Dune::UGGrid<3> >(suite, 3);
 #endif

dune/solvers/test/energynormtest.cc

@@ -145,6 +145,12 @@ struct CustomMultiTypeBlockVector : public Dune::MultiTypeBlockVector<Args...> {
 template<class... Args>
 struct CustomMultiTypeBlockMatrix : public Dune::MultiTypeBlockMatrix<Args...> {
   constexpr static size_t blocklevel = 1; // fake value needed for BCRSMatrix nesting
+
+  // HACK for istl compatibility
+  static constexpr size_t size(){
+    return sizeof...(Args);
+  }
+
   template<class K, typename = std::enable_if_t<Dune::IsNumber<K>::value>>
   void operator*=(K v) { Dune::Hybrid::forEach(*this, [v](auto& b) { b *= v; }); } // allow multiplication by scalar
 };
@@ -228,8 +234,8 @@ int main() try
   /// test with multi blocked vectors
   // types
   using MVector = CustomMultiTypeBlockVector<BVector, BVector>;
-  using MMatrix0 = CustomMultiTypeBlockMatrix<BMatrix, BMatrix>;
-  using MMatrix1 = CustomMultiTypeBlockMatrix<BMatrix, BMatrix>;
+  using MMatrix0 = CustomMultiTypeBlockVector<BMatrix, BMatrix>;
+  using MMatrix1 = CustomMultiTypeBlockVector<BMatrix, BMatrix>;
   using MMatrix = CustomMultiTypeBlockMatrix<MMatrix0, MMatrix1>;
   // instance setup
   using namespace Dune::Hybrid;
@@ -243,8 +249,8 @@ int main() try
   /// test with blocked multitype vectors
   // types
   using NVector = CustomMultiTypeBlockVector<FVector, FVector>;
-  using NMatrix0 = CustomMultiTypeBlockMatrix<FMatrix, FMatrix>;
-  using NMatrix1 = CustomMultiTypeBlockMatrix<FMatrix, FMatrix>;
+  using NMatrix0 = CustomMultiTypeBlockVector<FMatrix, FMatrix>;
+  using NMatrix1 = CustomMultiTypeBlockVector<FMatrix, FMatrix>;
   using NMatrix = CustomMultiTypeBlockMatrix<NMatrix0, NMatrix1>;
   constexpr size_t bnSize = 4;
   using BNVector = BlockVector<NVector>;