.gitlab-ci.yml

@@ -22,15 +22,6 @@ variables:
 # than one thread.
   OMP_NUM_THREADS: 1
 
-debian-11-gcc-9-17-with-checking:
-  image: registry.dune-project.org/docker/ci/debian:11
-  script: duneci-standard-test
-  variables:
-    DUNECI_TOOLCHAIN:   gcc-9-17
-    DUNECI_CXXFLAGS: "-DDUNE_ISTL_WITH_CHECKING"
-    DUNECI_CMAKE_FLAGS: "-DDUNE_MAX_TEST_CORES=4 -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE"
-  tags: [duneci]
-
 # Check for spelling mistakes in text
 code-spelling-check:
   stage: .pre

CHANGELOG.md

@@ -3,7 +3,14 @@ SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE
 SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 -->
 
-# Master (will become release 2.10)
+# Master (will become release 2.11)
+
+- `FastAMG` can now be used with non-nested matrices of type `BCRSMatrix<T>`
+  where `T` is a plain number type instead of `FieldMatrix<T,m,n>`.
+
+# Release 2.10
+
+- Improve testing support on Laplacian matrices with an optional diagonal regularization parameter.
 
 - Base the implementation of `VariableBlockVector` on `std::vector` as the storage type. Note that
   this prevents from using `bool` as block type that was possible before.

CMakeLists.txt

@@ -20,6 +20,12 @@ list(APPEND CMAKE_MODULE_PATH ${dune-common_MODULE_PATH}
 #include the dune macros
 include(DuneMacros)
 
+# deactivate global include-directories
+dune_policy(SET DP_DEFAULT_INCLUDE_DIRS NEW)
+
+# deactivate global calls to add_dune_all_flags in tests
+dune_policy(SET DP_TEST_ADD_ALL_FLAGS NEW)
+
 # start a dune project with information from dune.module
 dune_project()
 
@@ -28,6 +34,9 @@ dune_add_library(duneistl INTERFACE
   EXPORT_NAME ISTL
   LINK_LIBRARIES Dune::Common)
 
+# set include directories to target
+dune_default_include_directories(duneistl INTERFACE)
+
 add_subdirectory(cmake/modules)
 add_subdirectory(dune)
 add_subdirectory(doc)

LICENSE.md

@@ -3,37 +3,41 @@ Copyright holders:
 2014--2016    Marco Agnese
 2003--2013    Peter Bastian
 2020          Jean Benezech
-2004--2022    Markus Blatt
-2014--2022    Ansgar Burchardt
+2004--2024    Markus Blatt
+2022          Eduardo Bueno
+2022          Samuel Burbulla
+2014--2023    Ansgar Burchardt
 2004--2005    Adrian Burri
 2017--2018    Matthew Collins
-2008--2021    Andreas Dedner
+2008--2024    Andreas Dedner
 2018--2022    Nils-Arne Dreier
-2003--2022    Christian Engwer
+2003--2023    Christian Engwer
 2005--2019    Jorrit Fahlke
 2008--2017    Bernd Flemisch
 2017--2020    Janick Gerstenberger
 2015--2020    Felix Gruber
-2005--2021    Carsten Gräser
-2012--2021    Christoph Grüninger
+2005--2024    Carsten Gräser
+2012--2024    Christoph Grüninger
 2016--2020    René Heß
 2018          Claus-Justus Heine
 2016--2019    Stephan Hilb
-2018          Lasse Hinrichsen
+2018--2024    Lasse Hinrichsen
 2012--2013    Olaf Ippisch
-2018--2021    Patrick Jaap
+2018--2024    Patrick Jaap
 2013--2019    Dominic Kempf
 2015          Emmanouil Kiagias
-2004--2022    Robert Klöfkorn
-2016--2022    Timo Koch
+2004--2023    Robert Klöfkorn
+2016--2023    Timo Koch
 2007          Sreejith Pulloor Kuttanikkad
 2013--2016    Arne Morten Kvarving
 2009--2013    Andreas Lauser
 2012--2017    Tobias Malkmus
 2007--2009    Sven Marnach
 2013          René Milk
+2024          Alex Müller
 2013--2019    Steffen Müthing
 2016          Maikel Nadolski
+2023--2024    Lisa Julia Nebel
 2003--2005    Thimo Neubauer
 2010--2012    Rebecca Neumann
 2008--2018    Martin Nolte
@@ -41,13 +45,13 @@ Copyright holders:
 2013--2014    Marian Piatkowski
 2011--2016    Elias Pipping
 2013          Jurgis Pods
-2018--2022    Simon Praetorius
+2018--2024    Simon Praetorius
 2009          Atgeirr Rasmussen
 2018--2022    Lukas Renelt
-2021--2022    Santiago Ospina De Los Ríos
+2021--2024    Santiago Ospina De Los Ríos
 2008--2013    Uli Sack
 2023          Vaishnavi Sanchi
-2004--2021    Oliver Sander
+2004--2024    Oliver Sander
 2015--2019    Linus Seelinger
 2013          Bård Skaflestad
 2018          Mathis Springwald

cmake/modules/AddARPACKPPFlags.cmake

@@ -12,12 +12,19 @@
 #
 #       A list of targets to use ARPACKPP with.
 #
+include_guard(GLOBAL)
+
+set_package_properties("ARPACK" PROPERTIES
+  PURPOSE "Solve large scale eigenvalue problems")
+
+set_package_properties("ARPACKPP" PROPERTIES
+  PURPOSE "C++ interface for ARPACK")
 
 function(add_dune_arpackpp_flags _targets)
   if(ARPACKPP_FOUND)
     foreach(_target ${_targets})
       target_link_libraries(${_target} PUBLIC ${ARPACKPP_DUNE_LIBRARIES})
-      target_compile_definitions(${_target} PUBLIC ENABLE_ARPACKPP=1)
+      target_compile_definitions(${_target} PUBLIC HAVE_ARPACKPP=1)
       target_compile_options(${_target} PUBLIC ${ARPACKPP_DUNE_COMPILE_FLAGS})
     endforeach()
   endif()

cmake/modules/AddSuperLUFlags.cmake

@@ -12,6 +12,10 @@
 #
 #       A list of targets to use SuperLU with.
 #
+include_guard(GLOBAL)
+
+set_package_properties("SuperLU" PROPERTIES
+  PURPOSE "Direct solver for linear system, based on LU decomposition")
 
 # set HAVE_SUPERLU for config.h
 set(HAVE_SUPERLU ${SuperLU_FOUND})
@@ -19,7 +23,7 @@ set(HAVE_SUPERLU ${SuperLU_FOUND})
 # register all SuperLU related flags
 if(SuperLU_FOUND)
   dune_register_package_flags(
-    COMPILE_DEFINITIONS "ENABLE_SUPERLU=1"
+    COMPILE_DEFINITIONS "HAVE_SUPERLU=1"
     LIBRARIES SuperLU::SuperLU)
 endif()
 
@@ -28,7 +32,7 @@ function(add_dune_superlu_flags _targets)
   if(SuperLU_FOUND)
     foreach(_target ${_targets})
       target_link_libraries(${_target} PUBLIC SuperLU::SuperLU)
-      target_compile_definitions(${_target} PUBLIC ENABLE_SUPERLU=1)
+      target_compile_definitions(${_target} PUBLIC HAVE_SUPERLU=1)
     endforeach()
   endif()
 endfunction(add_dune_superlu_flags)

cmake/modules/FindARPACK.cmake

@@ -25,6 +25,11 @@
 #    system paths.
 #
 
+# text for feature summary
+set_package_properties("ARPACK" PROPERTIES
+  URL "https://www.arpack.org"
+  DESCRIPTION "ARnoldi PACKage")
+
 # look for library, only at positions given by the user
 find_library(ARPACK_LIBRARY
   NAMES "arpack"
@@ -81,8 +86,3 @@ else()
     "Determining location of ARPACK failed:\n"
     "Libraries to link against: ${ARPACK_LIBRARIES}\n\n")
 endif()
-
-# text for feature summary
-set_package_properties("ARPACK" PROPERTIES
-  DESCRIPTION "ARnoldi PACKage"
-  PURPOSE "Solve large scale eigenvalue problems")

cmake/modules/FindARPACKPP.cmake

@@ -28,6 +28,11 @@
 #    system paths.
 #
 
+# text for feature summary
+set_package_properties("ARPACKPP" PROPERTIES
+  URL "https://github.com/m-reuter/arpackpp"
+  DESCRIPTION "ARPACK++")
+
 # find ARPACK which is required by ARPACK++
 find_package(ARPACK)
 
@@ -123,17 +128,12 @@ else()
     "Libraries to link against: ${ARPACKPP_LIBRARIES}\n\n")
 endif()
 
-# set HAVE_ARPACKPP for config.h
+# set HAVE_ARPACKPP for conditional tests
 set(HAVE_ARPACKPP ${ARPACKPP_FOUND})
 
 # register all ARPACK++ related flags
 if(ARPACKPP_FOUND)
-  dune_register_package_flags(COMPILE_DEFINITIONS "ENABLE_ARPACKPP=1"
+  dune_register_package_flags(COMPILE_DEFINITIONS "HAVE_ARPACKPP=1"
                               LIBRARIES "${ARPACKPP_LIBRARIES}"
                               COMPILE_OPTIONS "${ARPACKPP_DUNE_COMPILE_FLAGS}")
 endif()
-
-# text for feature summary
-set_package_properties("ARPACKPP" PROPERTIES
-  DESCRIPTION "ARPACK++"
-  PURPOSE "C++ interface for ARPACK")

cmake/modules/FindSuperLU.cmake

@@ -25,8 +25,8 @@
 # text for feature summary
 include(FeatureSummary)
 set_package_properties("SuperLU" PROPERTIES
-  DESCRIPTION "Supernodal LU"
-  PURPOSE "Direct solver for linear system, based on LU decomposition")
+  URL "https://portal.nersc.gov/project/sparse/superlu/"
+  DESCRIPTION "Supernodal LU")
 
 set(SUPERLU_INT_TYPE "int" CACHE STRING
   "The integer version that SuperLU was compiled for (Default is int.

config.h.cmake

@@ -30,16 +30,10 @@
 
 /* end private */
 
-/* Define to ENABLE_SUPERLU if the SuperLU library is available */
-#cmakedefine HAVE_SUPERLU ENABLE_SUPERLU
-
 /* Define to the integer type that SuperLU was compiled for
    See e.g. what int_t is defined to in slu_sdefs.h */
 #cmakedefine SUPERLU_INT_TYPE @SUPERLU_INT_TYPE@
 
-/* Define to ENABLE_ARPACKPP if the ARPACK++ library is available */
-#cmakedefine HAVE_ARPACKPP ENABLE_ARPACKPP
-
 /* Define to the version of dune-istl */
 #define DUNE_ISTL_VERSION "${DUNE_ISTL_VERSION}"
 

dune.module

@@ -2,11 +2,11 @@
 # SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 Module: dune-istl
-Version: 2.10-git
+Version: 2.11-git
 Author: The Dune Core developers
 Maintainer: dune-devel@lists.dune-project.org
 Description: Iterative solver template library which provides generic sparse matrix/vector classes and a variety of solvers based on these classes
 URL: https://gitlab.dune-project.org/core/dune-istl
 Python-Requires:
-Depends: dune-common (>= 2.10)
+Depends: dune-common (>= 2.11)
 Whitespace-Hook: Yes

dune/istl/basearray.hh

@@ -5,11 +5,12 @@
 #ifndef DUNE_ISTL_BASEARRAY_HH
 #define DUNE_ISTL_BASEARRAY_HH
 
-#include "assert.h"
+#include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <memory>
 #include <algorithm>
+#include <type_traits>
 
 #include "istlexception.hh"
 #include <dune/common/iteratorfacades.hh>
@@ -97,17 +98,27 @@ namespace Imp {
       friend class RealIterator<ValueType>;
 
       //! constructor
-      RealIterator ()
-        : p(0), i(0)
-      {}
+      RealIterator () = default;
 
-      RealIterator (const B* _p, B* _i) : p(_p), i(_i)
-      {   }
+      RealIterator (const B* _p, B* _i)
+        : p(_p), i(_i)
+      {}
 
-      RealIterator(const RealIterator<ValueType>& it)
-        : p(it.p), i(it.i)
+      template <class T_,
+        std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
+      RealIterator (const RealIterator<T_>& other)
+        : p(other.p), i(other.i)
       {}
 
+      template <class T_,
+        std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
+      RealIterator& operator= (const RealIterator<T_>& other)
+      {
+        p = other.p;
+        i = other.i;
+        return *this;
+      }
+
       //! return index
       size_type index () const
       {
@@ -163,8 +174,8 @@ namespace Imp {
         i+=d;
       }
 
-      const B* p;
-      B* i;
+      const B* p = nullptr;
+      B* i = nullptr;
     };
 
     //! iterator type for sequential access
@@ -349,22 +360,29 @@ namespace Imp {
       friend class RealIterator<ValueType>;
 
       //! constructor
-      RealIterator ()
-        : p(0), j(0), i(0)
-      {}
+      RealIterator () = default;
 
       //! constructor
       RealIterator (B* _p, size_type* _j, size_type _i)
         : p(_p), j(_j), i(_i)
-      {       }
+      {}
 
-      /**
-       * @brief Copy constructor from mutable iterator
-       */
-      RealIterator(const RealIterator<ValueType>& it)
-        : p(it.p), j(it.j), i(it.i)
+      template <class T_,
+        std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
+      RealIterator (const RealIterator<T_>& other)
+        : p(other.p), j(other.j), i(other.i)
       {}
 
+      template <class T_,
+        std::enable_if_t<std::is_same_v<std::remove_const_t<T>, std::remove_const_t<T_>>, int> = 0>
+      RealIterator& operator= (const RealIterator<T_>& other)
+      {
+        p = other.p;
+        j = other.j;
+        i = other.i;
+        return *this;
+      }
+
 
       //! equality
       bool equals (const RealIterator<ValueType>& it) const
@@ -424,9 +442,9 @@ namespace Imp {
         return p[i];
       }
 
-      B* p;
-      size_type* j;
-      size_type i;
+      B* p = nullptr;
+      size_type* j = nullptr;
+      size_type i = 0;
     };
 
     /** @brief The iterator type. */

dune/istl/cholmod.hh

@@ -11,7 +11,7 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include<dune/istl/solver.hh>
-#include <dune/istl/solverfactory.hh>
+#include <dune/istl/solverregistry.hh>
 #include <dune/istl/foreach.hh>
 
 #include <vector>
@@ -250,7 +250,7 @@ public:
 
   /** @brief simple forward to apply(X&, Y&, InverseOperatorResult&)
     */
-  void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
+  void apply (Vector& x, Vector& b, [[maybe_unused]] double reduction, InverseOperatorResult& res) override
   {
     apply(x,b,res);
   }
@@ -260,7 +260,7 @@ public:
    * The method assumes that setMatrix() was called before
    * In the case of a given ignore field the corresponding entries of both in x and b will stay untouched in this method.
   */
-  void apply(Vector& x, Vector& b, InverseOperatorResult& res)
+  void apply(Vector& x, Vector& b, InverseOperatorResult& res) override
   {
     // do nothing if N=0
     if ( nIsZero_ )
@@ -461,6 +461,10 @@ public:
 
     });
 
+    // Remove old factor that may be left over from a previous run
+    if (L_)
+      CholmodMethod::free_factor(&L_, &c_);
+
     // Now analyse the pattern and optimal row order
     L_ = CholmodMethod::analyze(M.get(), &c_);
 
@@ -468,7 +472,7 @@ public:
     CholmodMethod::factorize(M.get(), L_, &c_);
   }
 
-  virtual SolverCategory::Category category() const
+  SolverCategory::Category category() const override
   {
     return SolverCategory::Category::sequential;
   }
@@ -524,34 +528,31 @@ private:
   std::vector<std::size_t> subIndices_;
 };
 
-  struct CholmodCreator{
-    template<class F> struct isValidBlock : std::false_type{};
-    template<int k> struct isValidBlock<FieldVector<double,k>> : std::true_type{};
-    template<int k> struct isValidBlock<FieldVector<float,k>> : std::true_type{};
-
-    template<class TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator()(TL /*tl*/, const M& mat, const Dune::ParameterTree& /*config*/,
-               std::enable_if_t<isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      using D = typename Dune::TypeListElement<1, TL>::type;
-      auto solver = std::make_shared<Dune::Cholmod<D>>();
-      solver->setMatrix(mat);
-      return solver;
-    }
-
-    // second version with SFINAE to validate the template parameters of Cholmod
-    template<typename TL, typename M>
-    std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
-                                          typename Dune::TypeListElement<2, TL>::type>>
-    operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
-                std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
-    {
-      DUNE_THROW(UnsupportedType, "Unsupported Type in Cholmod");
-    }
-  };
-  DUNE_REGISTER_DIRECT_SOLVER("cholmod", Dune::CholmodCreator());
+  DUNE_REGISTER_SOLVER("cholmod",
+                       [](auto opTraits, const auto& op, const Dune::ParameterTree& config)
+                       -> std::shared_ptr<typename decltype(opTraits)::solver_type>
+                       {
+                         using OpTraits = decltype(opTraits);
+                         using M = typename OpTraits::matrix_type;
+                         using D = typename OpTraits::domain_type;
+                         // works only for sequential operators
+                         if constexpr (OpTraits::isParallel){
+                           if(opTraits.getCommOrThrow(op).communicator().size() > 1)
+                             DUNE_THROW(Dune::InvalidStateException, "CholMod works only for sequential operators.");
+                         }
+                         if constexpr (OpTraits::isAssembled &&
+                                       // check whether the Matrix field_type is double or float
+                                       (std::is_same_v<typename FieldTraits<D>::field_type, double> ||
+                                       std::is_same_v<typename FieldTraits<D>::field_type, float>)){
+                           const auto& A = opTraits.getAssembledOpOrThrow(op);
+                           const M& mat = A->getmat();
+                           auto solver = std::make_shared<Dune::Cholmod<D>>();
+                           solver->setMatrix(mat);
+                           return solver;
+                         }
+                         DUNE_THROW(UnsupportedType,
+                                    "Unsupported Type in Cholmod (only double and float supported)");
+                       });
 
 } /* namespace Dune */
 

dune/istl/common/registry.hh

@@ -53,7 +53,7 @@ namespace Dune {
     }
 
     /*
-      Register all creators from the registry in the Parameterizedobjectfactory An
+      Register all creators from the registry in the Parameterizedobjectfactory. An
       object of V is passed in the creator and should be used to determine the
       template arguments.
     */

dune/istl/dilu.hh

@@ -9,6 +9,7 @@
 #include <complex>
 #include <map>
 #include <vector>
+#include <sstream>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/scalarvectorview.hh>
@@ -81,9 +82,13 @@ namespace Dune
         }
         catch (Dune::FMatrixError &e)
         {
-          DUNE_THROW(MatrixBlockError, "DILU failed to invert matrix block D[" << row_i << "]"
-                                                                               << e.what();
-                     th__ex.r = row_i;);
+          std::ostringstream sstream;
+          sstream << THROWSPEC(MatrixBlockError)
+            << "DILU failed to invert matrix block D[" << row_i << "]" << e.what();
+          MatrixBlockError ex;
+          ex.message(sstream.str());
+          ex.r = row_i;
+          throw ex;
         }
       }
     }

dune/istl/eigenvalue/poweriteration.hh

@@ -56,13 +56,13 @@ namespace Dune
           mutable_scaling_(mutable_scaling)
       {}
 
-      virtual void apply (const X& x, Y& y) const
+      void apply (const X& x, Y& y) const override
       {
         y = x;
         y *= immutable_scaling_*mutable_scaling_;
       }
 
-      virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+      void applyscaleadd (field_type alpha, const X& x, Y& y) const override
       {
         X temp(x);
         temp *= immutable_scaling_*mutable_scaling_;
@@ -70,7 +70,7 @@ namespace Dune
       }
 
       //! Category of the linear operator (see SolverCategory::Category)
-      virtual SolverCategory::Category category() const
+      SolverCategory::Category category() const override
       {
         return SolverCategory::sequential;
       }
@@ -108,14 +108,14 @@ namespace Dune
           "Range type of both operators doesn't match!");
       }
 
-      virtual void apply (const domain_type& x, range_type& y) const
+      void apply (const domain_type& x, range_type& y) const override
       {
         op1_.apply(x,y);
         op2_.applyscaleadd(1.0,x,y);
       }
 
-      virtual void applyscaleadd (field_type alpha,
-        const domain_type& x, range_type& y) const
+      void applyscaleadd (field_type alpha,
+        const domain_type& x, range_type& y) const override
       {
         range_type temp(y);
         op1_.apply(x,temp);
@@ -124,7 +124,7 @@ namespace Dune
       }
 
       //! Category of the linear operator (see SolverCategory::Category)
-      virtual SolverCategory::Category category() const
+      SolverCategory::Category category() const override
       {
         return SolverCategory::sequential;
       }

dune/istl/eigenvalue/test/CMakeLists.txt

 # SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
 # SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
+# Link all test targets in this directory against Dune::ISTL
+link_libraries(Dune::ISTL)
+
 dune_add_test(NAME poweriterationtest SOURCES cond2test.cc
   CMD_ARGS 40)
 
@@ -17,4 +20,6 @@ endif()
 
 if(ARPACKPP_FOUND AND SuperLU_FOUND)
   dune_add_test(NAME arpackppsuperlutest SOURCES cond2test.cc)
+  add_dune_superlu_flags(arpackppsuperlutest)
+  add_dune_arpackpp_flags(arpackppsuperlutest)
 endif()

dune/istl/eigenvalue/test/cond2test.cc

 // SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
-#include <config.h>
 
 #include <iostream>
 

dune/istl/foreach.hh

@@ -159,11 +159,11 @@ std::pair<std::size_t,std::size_t> flatMatrixForEach(Matrix&& matrix, F&& f, std
       // check whether we have valid sized blocks
       assert( ( blockRows!=0 or blockCols!=0 ) and "the block size can't be zero");
 
-      for ( auto rowIt = matrix.begin(); rowIt != matrix.end(); rowIt++ )
+      for (auto rowIt = matrix.begin(); rowIt != matrix.end(); ++rowIt)
       {
         auto&& row = *rowIt;
         auto rowIdx = rowIt.index();
-        for ( auto colIt = row.begin(); colIt != row.end(); colIt++ )
+        for (auto colIt = row.begin(); colIt != row.end(); colIt++)
         {
           auto&& entry = *colIt;
           auto colIdx = colIt.index();

dune/istl/ildl.hh

@@ -3,6 +3,8 @@
 #ifndef DUNE_ISTL_ILDL_HH
 #define DUNE_ISTL_ILDL_HH
 
+#include <sstream>
+
 #include <dune/common/scalarvectorview.hh>
 #include <dune/common/scalarmatrixview.hh>
 #include "ilu.hh"
@@ -135,7 +137,14 @@ namespace Dune
       }
       catch( const Dune::FMatrixError &e )
       {
-        DUNE_THROW( MatrixBlockError, "ILDL failed to invert matrix block A[" << i.index() << "][" << ij.index() << "]" << e.what(); th__ex.r = i.index(); th__ex.c = ij.index() );
+        std::ostringstream sstream;
+        sstream << THROWSPEC(MatrixBlockError)
+          << "ILDL failed to invert matrix block A[" << i.index() << "][" << ij.index() << "]" << e.what();
+        MatrixBlockError ex;
+        ex.message(sstream.str());
+          ex.r = i.index();
+          ex.c = ij.index();
+        throw ex;
       }
     }
   }