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>;
 
@@ -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/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/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/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/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/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();

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>;

dune/solvers/test/obstacletnnmgtest.cc

@@ -76,7 +76,7 @@ void solveObstacleProblemByTNNMG(const GridType& grid, const MatrixType& mat, Ve
     tnnmgStep.preprocess();
     tnnmgStep.nestedIteration();
 
-    // cretae loop solver
+    // create loop solver
     Solver solver(tnnmgStep, maxIterations, tolerance, norm, Solver::FULL);
 
     // solve problem

dune/solvers/test/proximalnewtonsolvertest.cc

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+
+#include <config.h>
+
+#include <cmath>
+#include <iostream>
+#include <sstream>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/norms/twonorm.hh>
+#include <dune/solvers/solvers/cholmodsolver.hh>
+#include <dune/solvers/solvers/proximalnewtonsolver.hh>
+
+#include "common.hh"
+
+using namespace Dune;
+
+// grid dimension
+const int dim = 3;
+
+// f(x) = 1/4 * ||x||^4
+template<class Matrix, class Vector>
+struct F
+{
+  void assembleGradientAndHessian( const Vector& x, Vector& grad, Matrix& hess ) const
+  {
+    auto norm2 = x.two_norm2();
+
+    grad = x;
+    grad *= norm2;
+
+    for(size_t i=0; i<x.size(); i++)
+      for(size_t j=0; j<x.size(); j++)
+        hess[i][j] = (i==j)*norm2 + 2*x[i]*x[j];
+  }
+
+  double computeEnergy( const Vector& x ) const
+  {
+    return  0.25* x.two_norm2() * x.two_norm2();
+  }
+};
+
+
+// g(x) = ( -1, -1, -1, ... , -1)^T * x
+template<class Vector>
+struct G
+{
+  double operator()( const Vector& x ) const
+  {
+    auto realX = x;
+    realX += offset_;
+
+    double sum = 0.0;
+    for ( size_t i=0; i<realX.size(); i++ )
+      sum -= realX[i];
+
+    return sum;
+  }
+
+  void updateOffset( const Vector& offset )
+  {
+    offset_ = offset;
+  }
+
+  Vector offset_;
+};
+
+template<class Matrix, class Vector>
+struct SecondOrdersolver
+{
+  using MatrixType = Matrix;
+
+  template<class BitVector>
+  void minimize( const Matrix& hess, const Vector& grad, const G<Vector>& g, double reg, Vector& dx, const BitVector& /*ignore*/) const
+  {
+    // add reg
+    auto HH = hess;
+    for ( size_t i=0; i<dx.size(); i++ )
+      HH[i][i] += reg;
+
+    // add g to grad
+    auto gg = grad;
+    gg *= -1.0;
+    for ( size_t i=0; i<dx.size(); i++ )
+      gg[i] += 1.0;
+
+
+    Solvers::CholmodSolver cholmod( HH, dx, gg );
+    cholmod.solve();
+  }
+};
+
+
+
+int main (int argc, char *argv[])
+{
+  // initialize MPI, finalize is done automatically on exit
+  [[maybe_unused]] MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+  const int size = 10;
+
+  using Vector = FieldVector<double,size>;
+  using Matrix = FieldMatrix<double,size>;
+  using BitVector = std::bitset<size>;
+
+  // this is the analytical solution of the problem
+  Vector exactSol, zeroVector;
+  exactSol = std::pow( size, -1.0/3.0 );
+
+  zeroVector = 0.0;
+
+  F<Matrix,Vector> f;
+  G<Vector> g;
+  g.updateOffset(zeroVector);
+  SecondOrdersolver<Matrix,Vector> sos;
+
+
+  // choose some random initial vector
+  Vector sol;
+  for ( size_t i=0; i<sol.size(); i++ )
+    sol[i] = i;
+
+  TwoNorm<Vector> norm;
+
+  // create the solver
+  Solvers::ProximalNewton::SimpleRegUpdater regUpdater;
+  double reg = 42.0;
+  Solvers::ProximalNewtonSolver proximalNewtonSolver( f, g, sos, sol, norm, regUpdater, reg, 100, 1e-14, Solver::FULL);
+
+  // set some empty ignore field
+  BitVector ignore;
+  proximalNewtonSolver.setIgnore( ignore );
+
+  // go!
+  proximalNewtonSolver.solve();
+
+  // compute diff the exact solution
+  sol -= exactSol;
+
+  std::cout << "Difference to exact solution = " << sol.two_norm() << std::endl;
+
+  bool passed = sol.two_norm() < 1e-15;
+
+  return passed ? 0 : 1;
+}

dune/solvers/test/umfpacksolvertest.cc

@@ -9,6 +9,8 @@
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/tuplevector.hh>
+#include <dune/common/version.hh>
 
 #include <dune/istl/bcrsmatrix.hh>
 
@@ -62,6 +64,99 @@ struct UMFPackSolverTestSuite
 
 
 
+/**  \brief test for the UMFPackSolver class with multitype matrix */
+template <size_t blocksize>
+struct UMFPackMultiTypeSolverTestSuite
+{
+  template <class GridType>
+  bool check(const GridType& grid)
+  {
+    double tol = 1e-11;
+
+    bool passed = true;
+
+    using Problem =
+    SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>;
+    Problem p(grid.leafGridView());
+
+    // create a multitype problem
+    using BlockRow = MultiTypeBlockVector<typename Problem::Matrix,typename Problem::Matrix>;
+    using Matrix = MultiTypeBlockMatrix<BlockRow,BlockRow>;
+
+    using Vector = MultiTypeBlockVector<typename Problem::Vector,typename Problem::Vector>;
+
+    Matrix A;
+    Vector u, rhs;
+
+    using namespace Dune::Indices;
+    // initialize each block
+    A[_0][_0] = p.A;
+    A[_0][_1] = p.A;
+    A[_1][_0] = p.A;
+    A[_1][_1] = p.A;
+
+    u[_0] = p.u;
+    u[_1] = p.u;
+    rhs[_0] = p.rhs;
+    rhs[_1] = p.rhs;
+
+    using IgnoreBlock = decltype(p.ignore);
+    using Ignore = TupleVector<IgnoreBlock,IgnoreBlock>;
+
+    Ignore ignore;
+    ignore[_0] = p.ignore;
+    ignore[_1] = p.ignore;
+
+    typedef Solvers::UMFPackSolver<Matrix,Vector> Solver;
+    Solver solver(A,u,rhs);
+
+    solver.setIgnore(ignore);
+
+    // solve blockdiagonal
+    A[_0][_1] *= 0.0;
+    A[_1][_0] *= 0.0;
+    solver.solve();
+
+    // the solution is ( u_ex, u_ex )
+    if (p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) > tol)
+    {
+      std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(u[_0],p.u_ex) + p.energyNorm.diff(u[_1],p.u_ex) << std::endl;
+      passed = false;
+    }
+
+    // solve a problem more generic problem without ignore field
+
+    A[_0][_1] = p.A;
+    A[_1][_0] = p.A;
+
+    // make A regular
+    A[_0][_0] *= 2.0;
+
+    u[_0] = p.u;
+    u[_1] = p.u;
+
+    // solve problem
+    Solver solver2(A,u,rhs);
+    solver2.preprocess();
+    solver2.solve();
+
+    // check residual
+    auto resisual = rhs;
+    A.mmv(u,resisual);
+
+    if (p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) > tol)
+    {
+      std::cout << "The UMFPackSolver did not produce a satisfactory result. ||residual||=" << p.energyNorm(resisual[_0]) + p.energyNorm(resisual[_1]) << std::endl;
+      passed = false;
+    }
+
+
+    return passed;
+  }
+};
+
+
+
 int main(int argc, char** argv)
 {
     Dune::MPIHelper::instance(argc, argv);
@@ -71,9 +166,15 @@ int main(int argc, char** argv)
     UMFPackSolverTestSuite<0> testsuite0;
     UMFPackSolverTestSuite<1> testsuite1;
     UMFPackSolverTestSuite<2> testsuite2;
+
+    UMFPackMultiTypeSolverTestSuite<3> testsuite3;
+
     passed = checkWithStandardGrids(testsuite0);
     passed = checkWithStandardGrids(testsuite1);
     passed = checkWithStandardGrids(testsuite2);
+#if DUNE_VERSION_GTE(DUNE_ISTL, 2, 10)
+    passed = checkWithStandardGrids(testsuite3);
+#endif
 
     return passed ? 0 : 1;
 }

dune/solvers/transferoperators/genericmultigridtransfer.hh

@@ -10,7 +10,7 @@
 
 #include <dune/geometry/referenceelements.hh>
 
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
 #include <dune/matrix-vector/axpy.hh>
 #include <dune/matrix-vector/transformmatrix.hh>
@@ -81,7 +81,7 @@ public:
         int cols = grid.size(cL, dim);
 
         // A factory for the shape functions
-        typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
+        typedef typename Dune::LagrangeLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
         typedef typename P1FECache::FiniteElementType FEType;
         P1FECache p1FECache;
 
@@ -228,7 +228,7 @@ public:
         int cols = gvCoarse.size(dim);
 
         // A factory for the shape functions
-        typedef typename Dune::PQkLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
+        typedef typename Dune::LagrangeLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
         typedef typename P1FECache::FiniteElementType FEType;
         P1FECache p1FECache;