.gitlab-ci.yml

 ---
 # Install external dependencies
-before_script:
-  - duneci-install-module https://git.imp.fu-berlin.de/agnumpde/dune-matrix-vector.git
 
-dune:git  clang  C++17:
-  image: registry.dune-project.org/docker/ci/dune:git-debian-10-clang-7-libcpp-17
+# For release-based jobs, we use corresponding images,
+# that already contain all core and staging dependencies.
+.release_based_job:
+  before_script:
+    - duneci-install-module https://gitlab.dune-project.org/fufem/dune-matrix-vector.git
+
+# For branch-based-based jobs, we need to install all dependencies manually.
+.branch_base_job:
+  before_script:
+    - . /duneci/bin/duneci-init-job
+    - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git
+    - duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git
+    - duneci-install-module https://gitlab.dune-project.org/core/dune-localfunctions.git
+    - duneci-install-module https://gitlab.dune-project.org/staging/dune-uggrid.git
+    - duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git
+    - duneci-install-module https://gitlab.dune-project.org/core/dune-istl.git
+    - duneci-install-module https://gitlab.dune-project.org/fufem/dune-matrix-vector.git
+
+# The 2.9 release is the one in Debian 12 ('bookworm', 'stable' at the time of writing)
+dune:2.9 debian-11 gcc-10 C++20:
+  extends: .release_based_job
+  variables:
+    DUNECI_BRANCH: releases/2.9
+  image: registry.dune-project.org/docker/ci/dune:2.9-debian-11-gcc-10-20
+  script: duneci-standard-test
+
+# The 2.10 release is the one in Debian 13 ('trixie', 'testing' at the time of writing)
+# To test against 2.10 we need an image with the corresponding version of the core modules.
+# Unfortunately, there is no 2.10 image with a recent enough compiler.
+dune:2.10 debian-11 gcc-10  C++20:
+  extends: .release_based_job
+  variables:
+    DUNECI_BRANCH: releases/2.10
+  image: registry.dune-project.org/docker/ci/dune:2.10-debian-11-gcc-10-20
   script: duneci-standard-test
 
-dune:git  gcc-9  C++20:
-  image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-9-20
+# Also test with the current development branch
+dune:git debian-11 clang-13  C++20:
+  extends: .branch_base_job
+#  image: registry.dune-project.org/docker/ci/dune:git-debian-11-clang-13-20
+  image: registry.dune-project.org/docker/ci/debian:11
+  variables:
+    DUNECI_TOOLCHAIN: clang-13-20
   script: duneci-standard-test
 
-dune:git  gcc-8  C++17:
-  image: registry.dune-project.org/docker/ci/dune:git-debian-10-gcc-8-17
+dune:git debian-11 gcc-10  C++20:
+  extends: .branch_base_job
+#  image: registry.dune-project.org/docker/ci/dune:git-debian-11-gcc-10-20
+  image: registry.dune-project.org/docker/ci/debian:11
+  variables:
+    DUNECI_TOOLCHAIN: gcc-10-20
   script: duneci-standard-test
 
 # Check for spelling mistakes in text

CHANGELOG.md

-# Master (will become release 2.9)
+# Master (will become release 2.11)
+
+- ...
+
+# Release 2.10
+
+- Deprecate the file `tuplevector.hh`.  An equivalent file exists in `dune-common`
+  since 2016. Please use that from now on.
+
+- `UMFPackSolver` accepts each correctly formed blocked matrix. The user has to make sure that the vector types of `x` and `rhs` are compatible to the matrix.
+  The main advantage is that it is now possible to use `MultiTypeBlockMatrix`.
+
+- A new solver `ProximalNewtonSolver` is added which solves non-smooth minimization problems.
+
+# Release 2.9
 
 - The internal matrix of the`EnergyNorm` can now be accessed by `getMatrix()`.
 
+- The default `BitVectorType` of the class `IterationStep` is now
+  `Solvers::DefaultBitVector_t<VectorType>` rather than `Dune::BitSetVector`.
+  This should do the right thing in more situations, while being fully
+  backward-compatible.
+
 - `codespell` spell checker is now active for automated spell checking in the Gitlab CI. 
   To exclude false positives add the words to the `--ignore-words-list` in `.gitlab-ci.yml`.
 

CMakeLists.txt

-cmake_minimum_required(VERSION 2.8.6)
+if(dune-common_VERSION VERSION_GREATER_EQUAL 2.10.0)
+    cmake_minimum_required(VERSION 3.16)
+else()
+    cmake_minimum_required(VERSION 3.13)
+endif()
 project(dune-solvers CXX)
 
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
@@ -22,9 +26,17 @@ dune_project()
 find_package(SuiteSparse OPTIONAL_COMPONENTS UMFPACK)
 include(AddSuiteSparseFlags)
 
+# Create library target and export it as Dune::Solvers
+dune_add_library(dunesolvers EXPORT_NAME Solvers LINK_LIBRARIES ${DUNE_LIBS})
+
+dune_register_package_flags(LIBRARIES dunesolvers)
+
 add_subdirectory("dune")
 add_subdirectory("doc")
 add_subdirectory("cmake/modules")
 
-# finalize the dune project, e.g. generating config.h etc.
-finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
+if(dune-common_VERSION VERSION_GREATER_EQUAL 2.10.0)
+    finalize_dune_project()
+else()
+    finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
+endif()

cmake/modules/FindIPOpt.cmake

@@ -9,12 +9,12 @@ find_library(DL_LIBRARY dl)
 find_path(IPOPT_INCLUDE_DIR
   NAMES "IpNLP.hpp"
   PATHS ${IPOPT_ROOT}
-  PATH_SUFFIXES "include" "include/coin"
+  PATH_SUFFIXES "include" "include/coin" "include/coin-or"
   NO_DEFAULT_PATH
 )
 find_path(IPOPT_INCLUDE_DIR
   NAMES "IpNLP.hpp"
-  PATH_SUFFIXES "include" "include/coin"
+  PATH_SUFFIXES "include" "include/coin" "include/coin-or"
 )
 
 find_library(IPOPT_LIBRARY
@@ -40,7 +40,7 @@ find_library(HSL_LIBRARY
 
 find_package_handle_standard_args(hsl DEFAULT_MSG HSL_LIBRARY)
 find_package_handle_standard_args(dl DEFAULT_MSG DL_LIBRARY)
-find_package_handle_standard_args(Ipopt DEFAULT_MSG IPOPT_INCLUDE_DIR IPOPT_LIBRARY)
+find_package_handle_standard_args(IPOpt DEFAULT_MSG IPOPT_INCLUDE_DIR IPOPT_LIBRARY)
 
 if(IPOPT_FOUND)
     set(HAVE_IPOPT ENABLE_IPOPT)

dune.module

@@ -4,7 +4,7 @@
 
 #Name of the module
 Module: dune-solvers
-Version: 2.9-git
+Version: 2.11-git
 Maintainer: oliver.sander@tu-dresden.de
 #depending on
 Depends: dune-common dune-grid dune-istl dune-localfunctions dune-matrix-vector

dune/solvers/CMakeLists.txt

-dune_add_library("dunesolvers"
-    iterationsteps/blockgssteps.cc
-    solvers/criterion.cc)
-
-dune_register_package_flags(LIBRARIES dunesolvers)
-
 add_subdirectory("common")
 add_subdirectory("iterationsteps")
 add_subdirectory("norms")
@@ -12,6 +6,10 @@ add_subdirectory("solvers")
 add_subdirectory("test")
 add_subdirectory("transferoperators")
 
+target_sources(dunesolvers PRIVATE
+    iterationsteps/blockgssteps.cc
+    solvers/criterion.cc)
+
 install(FILES
     computeenergy.hh
     DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/solvers)

dune/solvers/common/CMakeLists.txt

 install(FILES
-    arithmetic.hh
     boxconstraint.hh
     canignore.hh
     copyorreference.hh

dune/solvers/common/arithmetic.hh

-// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=8 sw=4 sts=4:
-#ifndef ARITHMETIC_HH
-#define ARITHMETIC_HH
-
-#warning arithmetic.hh is deprecated, please use dune-matrix-vector instead!
-
-#include <dune/matrix-vector/axpy.hh>
-#include <dune/matrix-vector/axy.hh>
-#include <dune/matrix-vector/crossproduct.hh>
-#include <dune/matrix-vector/matrixtraits.hh>
-#include <dune/matrix-vector/promote.hh>
-#include <dune/matrix-vector/scalartraits.hh>
-#include <dune/matrix-vector/transpose.hh>
-
-namespace Arithmetic
-{
-    using Dune::MatrixVector::Promote;
-    using Dune::MatrixVector::ScalarTraits;
-    using Dune::MatrixVector::MatrixTraits;
-    using Dune::MatrixVector::Transposed;
-    using Dune::MatrixVector::addProduct;
-    using Dune::MatrixVector::subtractProduct;
-    using Dune::MatrixVector::crossProduct;
-    using Dune::MatrixVector::transpose;
-    using Dune::MatrixVector::Axy;
-    using Dune::MatrixVector::bmAxy;
-}
-
-#endif

dune/solvers/common/defaultbitvector.hh

@@ -10,7 +10,7 @@
 #include <dune/istl/bvector.hh>
 #include <dune/istl/multitypeblockvector.hh>
 
-#include <dune/solvers/common/tuplevector.hh>
+#include <dune/common/tuplevector.hh>
 
 namespace Dune {
 namespace Solvers {

dune/solvers/common/staticmatrixtools.hh

@@ -16,8 +16,6 @@
 #include <dune/matrix-vector/scalartraits.hh>
 #include <dune/matrix-vector/transformmatrix.hh>
 
-#include "arithmetic.hh"
-
 namespace StaticMatrix
 {
         // type promotion (smallest matrix that can hold the sum of two matrices *******************

dune/solvers/common/tuplevector.hh

@@ -3,9 +3,9 @@
 #ifndef DUNE_SOLVERS_COMMON_TUPLEVECTOR_HH
 #define DUNE_SOLVERS_COMMON_TUPLEVECTOR_HH
 
-#include <tuple>
+#warning This file is deprecated!  Use tuplevector.hh from the dune-common module instead!
 
-#include <dune/common/indices.hh>
+#include <dune/common/tuplevector.hh>
 
 namespace Dune
 {
@@ -16,38 +16,7 @@ namespace Solvers
    * \brief A std::tuple that allows access to its element via operator[]
    */
   template<class... T>
-  class TupleVector : public std::tuple<T...>
-  {
-    using Base = std::tuple<T...>;
-
-  public:
-
-    // inherit all constructors from std::tuple
-    using Base::Base;
-
-    /** \brief Const access to the tuple elements */
-    template<std::size_t i>
-    auto operator[](const Dune::index_constant<i>&) const
-      ->decltype(std::get<i>(*this))
-    {
-      return std::get<i>(*this);
-    }
-
-    /** \brief Non-const access to the tuple elements */
-    template<std::size_t i>
-    auto operator[](const Dune::index_constant<i>&)
-      ->decltype(std::get<i>(*this))
-    {
-      return std::get<i>(*this);
-    }
-
-    /** \brief Number of elements of the tuple */
-    static constexpr std::size_t size()
-    {
-      return std::tuple_size<Base>::value;
-    }
-
-  };
+  using TupleVector = Dune::TupleVector<T...>;
 
 }  // namespace Functions
 

dune/solvers/iterationsteps/blockgssteps.hh

@@ -3,6 +3,7 @@
 
 #include <cassert>
 #include <functional>
+#include <type_traits>
 
 #include <dune/common/parametertree.hh>
 
@@ -181,7 +182,7 @@ template <class LinearSolver>
 auto truncateSymmetrically(LinearSolver&& linearSolver) {
   return [linearSolver = std::move(linearSolver) ](
       const auto& m, const auto& b, const auto& ignore) {
-    using Return = typename std::result_of<LinearSolver(decltype(m), decltype(b))>::type;
+    using Return = std::invoke_result_t<LinearSolver, decltype(m), decltype(b)>;
     if (ignore.all())
       return Return(0);
 

dune/solvers/iterationsteps/iterationstep.hh

@@ -6,8 +6,8 @@
 #include <vector>
 #include <string>
 
-#include <dune/common/bitsetvector.hh>
 #include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
 #include <dune/solvers/common/numproc.hh>
 
 namespace Dune {
@@ -15,7 +15,7 @@ namespace Dune {
 namespace Solvers {
 
     //! Base class for iteration steps being called by an iterative solver
-template<class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
+template<class VectorType, class BitVectorType = Solvers::DefaultBitVector_t<VectorType> >
 class IterationStep : virtual public NumProc, public CanIgnore<BitVectorType>
     {
     public:

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

@@ -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/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 ] ];
+        }
+      });
     }
   }