.gitlab-ci.yml

@@ -22,32 +22,16 @@ variables:
 # than one thread.
   OMP_NUM_THREADS: 1
 
-debian-11-gcc-9-17-with-checking:
+# Check for spelling mistakes in text
+code-spelling-check:
+  stage: .pre
+  # Avoid the global 'before_script'
+  before_script: ""
   image: registry.dune-project.org/docker/ci/debian:11
-  script: duneci-standard-test
-  variables:
-    DUNECI_TOOLCHAIN:   gcc-9-17
-    DUNECI_CMAKE_FLAGS:   "CC=gcc-9 CXX=g++-9 -DCMAKE_CXX_FLAGS='-std=c++17 -O2 -g -Wall -fdiagnostics-color=always -DDUNE_ISTL_WITH_CHECKING' -DDUNE_MAX_TEST_CORES=4 -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE"
-  tags: [duneci]
-
-# no numpy version available for this image on the gitlab registry
-# disable python bindings for now
-ubuntu:18.04 clang-6-17:
-  image: registry.dune-project.org/docker/ci/ubuntu:18.04
-  script: duneci-standard-test
-  stage: test
-  variables:
-    DUNECI_TOOLCHAIN: clang-6-17
-    DUNECI_CMAKE_FLAGS: "-DDUNE_ENABLE_PYTHONBINDINGS=OFF"
-  tags: [duneci]
-ubuntu:18.04 clang-5-17:
-  image: registry.dune-project.org/docker/ci/ubuntu:18.04
-  script: duneci-standard-test
-  stage: test
-  variables:
-    DUNECI_TOOLCHAIN: clang-5-17
-    DUNECI_CMAKE_FLAGS: "-DDUNE_ENABLE_PYTHONBINDINGS=OFF"
   tags: [duneci]
+  script:
+  - codespell
+    --ignore-words-list ba,eiter,equil,fo,parms
 
 reuse:
   stage: .pre

CHANGELOG.md

@@ -3,7 +3,64 @@ 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.9)
+# 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.
+
+- A method `BCRSMatrix::setIndicesNoSort()` was added. Similar
+  to `BCRSMatrix::setIndices()` this allows to insert all indices
+  of a row at once, but - incontrast to the latter - does not sort them.
+  This allows to speed up insertion if indices are already sorted.
+
+- `UMFPACK` is extended to arbitrary blocked matrix structures. This includes `MultiTypeBlockMatrix`. The external interface is unchanged.
+  However, internally the `BCCSMatrixInitializer` is replaced by direct calls of `flatMatrixForEach` similar to `Cholmod`. This requires a
+  compatible vector field of "ignored" degrees of freedom. The method `setSubMatrix` with a top-level index set is preserved for compatibility.
+
+- The internal storage in `MatrixIndexSet` was changed from `std::set` to a sorted `std::vector`
+  (with a `std::set` fallback for very dense rows) to improve performance. The stored index
+  type was changed from `std::size_t` to `uint32_t` to reduce memory consumption and improve
+  performance. Hence, `MatrixIndexSet` can no longer be used for very large matrices with more
+  than 2^32 columns.
+
+- Added flag 'useFixedOrder' to the coarsen method of AMGs ParallelIndicesCoarsener.
+  If set to true, during the creation of the coarser matrix (by accumulation and restriction
+  to fewer processes) the coarse matrix rows are ordered in a fixed manner, making parallel runs
+  reproducible; the runtime is possibly not ideal though.
+  If set to false (which is the default), the row order depends on the order of messages received from the
+  processes responsible for the respective parts of the finer grid. Then the indices on the coarser grid
+  may differ from run to run.
+
+- Define `field_type` and `real_type` in `MultiTypeBlock[Matrix|Vector]` only if a common type of these
+  types exist over all blocks in the container. Otherwise it is defined as `Std::nonesuch`.
+
+## Deprecations and removals
+
+- The deprecated CMake variables `SUPERLU_INCLUDE_DIRS`, `SUPERLU_LIBRARIES`,
+  `SUPERLU_DUNE_COMPILE_FLAGS`, and `SUPERLU_DUNE_LIBRARIES` are removed, they are
+  replaced by the target `SuperLU::SuperLU`.
+
+- Multiple containers no longer provide the deprecated member function `blocklevel`. Use the free function
+  `blockLevel()`.
+
+- Removed deprecated function `initSolverFactories()`, use `initSolverFactories<O>` instead.
+
+- Removed deprecated `MultyTypeBlockMatrix::size()`, use `N()` instead.
+
+- Removed deprecated `MultyTypeBlockVector::count()`, use `N()` instead.
+
+- Deprecated `writeSVGMatrix` with `std::stream` as the second argument. Use the method
+  with `std::stream` as the first argument.
+
+
+# Release 2.9
 
 - Add `const` qualifier to `LinearOperator` and `ScalarProduct` in
   `IterativeSolver`. In particular, the constructors of iterative solvers have
@@ -29,6 +86,11 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
   used by `CHOLMOD` itself to store the Cholesky factor.  This can be used to
   use the more advanced features of `CHOLMOD`.
 
+- The `Cholmod` class now has an additional template parameter `Index`, which specifies the
+  type used by `CHOLMOD` for integers.  The types `int` (the default) and `SuiteSparse_long`
+  (which is usually `long int`) are supported.  Use `SuiteSparse_long` if your problems
+  are too big for `int`.
+
 - You can now multiply objects of type `ScaledIdentityMatrix` by scalars
   using `operator*`.
 
@@ -47,6 +109,7 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 - Remove deprecated `ImplicitModeOverflowExhausted`, use
   `ImplicitModeCompressionBufferExhausted` instead.
 
+
 # Release 2.8
 
 - Extended the MatrixMarket IO functions for reading and writing vectors with
@@ -190,4 +253,4 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 - The solver infrastructure was updated to support SIMD data types (see
   current changes in `dune-common`). This allows to solve multiple systems
-  simultaniously using vectorization.
+  simultaneously using vectorization.

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
 
-# set up project
+cmake_minimum_required(VERSION 3.16)
 project(dune-istl C CXX)
 
-# general stuff
-cmake_minimum_required(VERSION 3.13)
-
 # guess build tree of dune-common
 if(NOT (dune-common_DIR OR dune-common_ROOT OR
         "${CMAKE_PREFIX_PATH}" MATCHES ".*dune-common.*"))
@@ -23,9 +20,23 @@ 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()
 
+# create library target and export it as Dune::ISTL
+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)
@@ -37,4 +48,4 @@ if( DUNE_ENABLE_PYTHONBINDINGS )
 endif()
 
 # finalize the dune project, e.g. generating config.h etc.
-finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
+finalize_dune_project()

LICENSE.md

@@ -3,37 +3,41 @@ Copyright holders:
 2014--2016    Marco Agnese
 2003--2013    Peter Bastian
 2020          Jean Benezech
-2004--2020    Markus Blatt
-2014--2020    Ansgar Burchardt
+2004--2024    Markus Blatt
+2022          Eduardo Bueno
+2022          Samuel Burbulla
+2014--2023    Ansgar Burchardt
 2004--2005    Adrian Burri
-2008--2021    Andreas Dedner
-2018--2021    Nils-Arne Dreier
-2003--2021    Christian Engwer
+2017--2018    Matthew Collins
+2008--2024    Andreas Dedner
+2018--2022    Nils-Arne Dreier
+2003--2023    Christian Engwer
 2005--2019    Jorrit Fahlke
 2008--2017    Bernd Flemisch
 2017--2020    Janick Gerstenberger
-2005--2020    Carsten Gräser
 2015--2020    Felix Gruber
-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
-2018          Lukas Renelt
-2004--2021    Robert Klöfkorn
-2016--2021    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,15 +45,19 @@ Copyright holders:
 2013--2014    Marian Piatkowski
 2011--2016    Elias Pipping
 2013          Jurgis Pods
-2018--2021    Simon Praetorius
+2018--2024    Simon Praetorius
 2009          Atgeirr Rasmussen
+2018--2022    Lukas Renelt
+2021--2024    Santiago Ospina De Los Ríos
 2008--2013    Uli Sack
-2004--2021    Oliver Sander
-2017--2019    Linus Seelinger
+2023          Vaishnavi Sanchi
+2004--2024    Oliver Sander
+2015--2019    Linus Seelinger
 2013          Bård Skaflestad
 2018          Mathis Springwald
-2019          Kilian Weishaupt
+2023          Jakob Torben
 2006--2010    Martin Weiser
+2019--2020    Kilian Weishaupt
 2015          Sebastian Westerheide
 2011--2012    Matthias Wohlmuth
 2014--2016    Jonathan Youett

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,17 +23,8 @@ 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)
-  # for backwards compatibility only,
-  # they are deprecated and will be removed after Dune 2.8
-  set(SUPERLU_INCLUDE_DIRS ${SUPERLU_INCLUDE_DIR})
-  set(SUPERLU_LIBRARIES    ${SUPERLU_LIBRARY})
-  set(SUPERLU_FOUND ${SuperLU_FOUND})
-  set(SUPERLU_DUNE_COMPILE_FLAGS "-I${SUPERLU_INCLUDE_DIRS}"
-    CACHE STRING "Compile flags used by DUNE when compiling SuperLU programs")
-  set(SUPERLU_DUNE_LIBRARIES ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES}
-    CACHE STRING "Libraries used by DUNE when linking SuperLU programs")
 endif()
 
 # Provide function to set target properties for linking to SuperLU
@@ -37,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/DuneIstlMacros.cmake

@@ -13,8 +13,6 @@ find_package(SuperLU 5.0)
 include(AddSuperLUFlags)
 find_package(ARPACKPP)
 include(AddARPACKPPFlags)
-find_package(SuiteSparse OPTIONAL_COMPONENTS CHOLMOD LDL SPQR UMFPACK)
-include(AddSuiteSparseFlags)
 
 # enable / disable backwards compatibility w.r.t. category
 set(DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE 1

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"
@@ -71,18 +76,13 @@ if(ARPACK_FOUND)
   set(ARPACK_LIBRARIES ${ARPACK_LIBRARY})
   # log result
   file(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeOutput.log
-    "Determing location of ARPACK succeeded:\n"
+    "Determining location of ARPACK succeeded:\n"
     "Libraries to link against: ${ARPACK_LIBRARIES}\n\n")
   set(ARPACK_DUNE_LIBRARIES ${ARPACK_LIBRARIES}
     CACHE STRING "Libraries used by DUNE when linking ARPACK programs")
 else()
-  # log errornous result
+  # log erroneous result
   file(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log
-    "Determing location of ARPACK failed:\n"
+    "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)
 
@@ -100,7 +105,7 @@ if(ARPACKPP_FOUND)
   endif(ARPACKPP_LIBRARY)
   # log result
   file(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeOutput.log
-    "Determing location of ARPACK++ succeeded:\n"
+    "Determining location of ARPACK++ succeeded:\n"
     "Include directory: ${ARPACKPP_INCLUDE_DIRS}\n"
     "Libraries to link against: ${ARPACKPP_LIBRARIES}\n\n")
 
@@ -116,24 +121,19 @@ if(ARPACKPP_FOUND)
   set(ARPACKPP_DUNE_LIBRARIES ${ARPACKPP_LIBRARIES}
     CACHE STRING "Libraries used by DUNE when linking ARPACK++ programs")
 else()
-  # log errornous result
+  # log erroneous result
   file(APPEND ${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeError.log
-    "Determing location of ARPACK++ failed:\n"
+    "Determining location of ARPACK++ failed:\n"
     "Include directory: ${ARPACKPP_INCLUDE_DIRS}\n"
     "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

@@ -10,17 +10,6 @@
 #    :code:`SuperLU_FOUND`
 #       True if SuperLU available and usable.
 #
-#    :code:`SUPERLU_INCLUDE_DIRS`
-#       Path to the SuperLU include dirs.
-#       .. deprecated:: 2.8
-#          Use target SuperLU::SuperLU instead.
-#
-#    :code:`SUPERLU_LIBRARIES`
-#       Name to the SuperLU library.
-#       .. deprecated:: 2.8
-#          Use target SuperLU::SuperLU instead.
-#
-#
 #    This module provides the following imported targets, if found:
 #
 #    :code:`SuperLU:SuperLU`
@@ -36,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}"
 

doc/istl.tex

@@ -649,7 +649,7 @@ contains the computed update to
 the current guess, i.~e. $\vec v = M^{-1} \vec
 d$ where $M$ is the approximate inverse of the operator $A$
 characterizing the preconditioner. The method \lstinline!void post(X& x)!
-should be called after all computations to give the precondtioner the
+should be called after all computations to give the preconditioner the
 chance to clean allocated resources. 
 
 See Table \ref{tab:precond} for a list of available

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.9)
+Depends: dune-common (>= 2.11)
 Whitespace-Hook: Yes

dune/istl/CMakeLists.txt

@@ -18,6 +18,7 @@ install(FILES
    btdmatrix.hh
    bvector.hh
    cholmod.hh
+   dilu.hh
    foreach.hh
    gsetc.hh
    ildl.hh

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>
@@ -32,7 +33,7 @@ namespace Imp {
        -  random access
 
            This container has *NO* memory management at all,
-           copy constuctor, assignment and destructor are all default.
+           copy constructor, assignment and destructor are all default.
 
            The constructor is made protected to emphasize that objects
        are only usable in derived classes.
@@ -41,13 +42,9 @@ namespace Imp {
            Setting the compile time switch DUNE_ISTL_WITH_CHECKING
            enables error checking.
 
-   \todo There shouldn't be an allocator argument here, because the array is 'unmanaged'.
-         And indeed, of the allocator, only its size_type is used.  Hence, the signature
-         of this class should be changed to <class B, int stype>
-
    \internal This class is an implementation detail, and should not be used outside of dune-istl.
    */
-  template<class B, class A=std::allocator<B> >
+  template<class B, class ST=std::size_t >
   class base_array_unmanaged
   {
   public:
@@ -57,11 +54,8 @@ namespace Imp {
     //! export the type representing the components
     typedef B member_type;
 
-    //! export the allocator type
-    typedef A allocator_type;
-
     //! the type for the index access
-    typedef typename A::size_type size_type;
+    typedef ST size_type;
 
     //! the type used for references
     using reference = B&;
@@ -104,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
       {
@@ -170,8 +174,8 @@ namespace Imp {
         i+=d;
       }
 
-      const B* p;
-      B* i;
+      const B* p = nullptr;
+      B* i = nullptr;
     };
 
     //! iterator type for sequential access
@@ -291,7 +295,7 @@ namespace Imp {
            -  find returning iterator
 
            This container has *NO* memory management at all,
-           copy constuctor, assignment and destructor are all default.
+           copy constructor, assignment and destructor are all default.
 
            The constructor is made protected to emphasize that objects
        are only usably in derived classes.
@@ -302,7 +306,7 @@ namespace Imp {
 
     \internal This class is an implementation detail, and should not be used outside of dune-istl.
    */
-  template<class B, class A=std::allocator<B> >
+  template<class B, class ST=std::size_t >
   class compressed_base_array_unmanaged
   {
   public:
@@ -312,11 +316,8 @@ namespace Imp {
     //! export the type representing the components
     typedef B member_type;
 
-    //! export the allocator type
-    typedef A allocator_type;
-
     //! The type used for the index access
-    typedef typename A::size_type size_type;
+    typedef ST size_type;
 
     //! the type used for references
     using reference = B&;
@@ -359,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
@@ -434,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/bcrsmatrix.hh

@@ -408,8 +408,8 @@ namespace Dune {
      The compress() method returns a value of type Dune::CompressionStatistics, which
      you can inspect to tune the construction parameters `_avg` and `compressionBufferSize`.
 
-     Use of copy constructor, assignment operator and matrix vector arithmetics
-     are not supported until the matrix is fully built.
+     Use of copy constructor, assignment operator and matrix vector arithmetic
+     is not supported until the matrix is fully built.
 
      The following sample code constructs a \f$ 10 \times 10\f$ matrix, with an expected
      number of two entries per matrix row.  The compression buffer size is set to 0.4.
@@ -493,19 +493,15 @@ namespace Dune {
     //! export the allocator type
     typedef A allocator_type;
 
-    //! implement row_type with compressed vector
-    typedef Imp::CompressedBlockVectorWindow<B,A> row_type;
-
     //! The type for the index access and the size
     typedef typename A::size_type size_type;
 
+    //! implement row_type with compressed vector
+    typedef Imp::CompressedBlockVectorWindow<B,size_type> row_type;
+
     //! The type for the statistics object returned by compress()
     typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
 
-    //! increment block level counter
-    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
-    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
-
     //! we support two modes
     enum BuildMode {
       /**
@@ -779,6 +775,7 @@ namespace Dune {
      * @param _avg expected average number of entries per row
      * @param compressionBufferSize fraction of _n*_avg which is expected to be
      *   needed for elements that exceed _avg entries per row.
+     * \param bm Matrix build mode
      *
      */
     BCRSMatrix (size_type _n, size_type _m, size_type _avg, double compressionBufferSize, BuildMode bm)
@@ -1222,11 +1219,35 @@ namespace Dune {
       *pos = col;
     }
 
+    //! Set all column indices for row from the given iterator range.
+    /**
+     * The iterator range has to be of the same length as the previously set row size.
+     * The entries in the iterator range must be sorted and must not contain duplicate values.
+     * This method will insert the indices in the given order.
+     *
+     * Calling this method overwrites any previously set column indices!
+     *
+     */
+    template<typename It>
+    void setIndicesNoSort(size_type row, It begin, It end)
+    {
+      size_type row_size = r[row].size();
+      size_type* col_begin = r[row].getindexptr();
+      size_type* col_end;
+      // consistency check between allocated row size and number of passed column indices
+      if ((col_end = std::copy(begin,end,r[row].getindexptr())) != col_begin + row_size)
+        DUNE_THROW(BCRSMatrixError,"Given size of row " << row
+                   << " (" << row_size
+                   << ") does not match number of passed entries (" << (col_end - col_begin) << ")");
+    }
+
+
     //! Set all column indices for row from the given iterator range.
     /**
      * The iterator range has to be of the same length as the previously set row size.
      * The entries in the iterator range do not have to be in any particular order, but
      * must not contain duplicate values.
+     * This method will insert the indices and sort them afterwards.
      *
      * Calling this method overwrites any previously set column indices!
      */
@@ -1467,10 +1488,18 @@ namespace Dune {
       overflow.clear();
 
       //determine average number of entries and memory usage
-      std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
-      nnz_ = diff;
-      stats.avg = (double) (nnz_) / (double) n;
-      stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
+      if ( n == 0)
+      {
+        stats.avg = 0;
+        stats.mem_ratio = 1;
+      }
+      else
+      {
+        std::ptrdiff_t diff = (r[n-1].getindexptr() + r[n-1].getsize() - j_.get());
+        nnz_ = diff;
+        stats.avg = (double) (nnz_) / (double) n;
+        stats.mem_ratio = (double) (nnz_) / (double) allocationSize_;
+      }
 
       //matrix is now built
       ready = built;
@@ -2195,12 +2224,13 @@ namespace Dune {
      * setup matrix) results in all the structure and data being lost. Please
      * call deallocate() before calling allocate in this case.
      *
-     * @param row The number of rows the matrix should contain.
+     * @param rows The number of rows the matrix should contain.
      * @param columns the number of columns the matrix should contain.
      * @param allocationSize The number of nonzero entries the matrix should hold (if omitted
      * defaults to 0).
-     * @param allocateRow Whether we have to allocate the row pointers, too. (Defaults to
+     * @param allocateRows Whether we have to allocate the row pointers, too. (Defaults to
      * true)
+     * \param allocate_data Whether data array needs to be allocated
      */
     void allocate(size_type rows, size_type columns, size_type allocationSize, bool allocateRows, bool allocate_data)
     {

dune/istl/bdmatrix.hh

@@ -50,10 +50,6 @@ namespace Dune {
     //! The type for the index access and the size
     typedef typename A::size_type size_type;
 
-    //! increment block level counter
-    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
-    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
-
     /** \brief Default constructor */
     BDMatrix() : BCRSMatrix<B,A>() {}
 

dune/istl/btdmatrix.hh

@@ -48,10 +48,6 @@ namespace Dune {
     //! The type for the index access and the size
     typedef typename A::size_type size_type;
 
-    //! increment block level counter
-    [[deprecated("Use free blockLevel function. Will be removed after 2.8.")]]
-    static constexpr auto blocklevel = blockLevel<B>()+1;
-
     /** \brief Default constructor */
     BTDMatrix() : BCRSMatrix<B,A>() {}
 

dune/istl/bvector.hh

@@ -80,8 +80,8 @@ namespace Imp {
 
    \internal This class is an implementation detail, and should not be used outside of dune-istl.
    */
-  template<class B, class A=std::allocator<B> >
-  class block_vector_unmanaged : public base_array_unmanaged<B,A>
+  template<class B, class ST=std::size_t >
+  class block_vector_unmanaged : public base_array_unmanaged<B,ST>
   {
   public:
 
@@ -91,17 +91,14 @@ namespace Imp {
     //! export the type representing the components
     typedef B block_type;
 
-    //! export the allocator type
-    typedef A allocator_type;
-
     //! The size type for the index access
-    typedef typename A::size_type size_type;
+    typedef ST size_type;
 
     //! make iterators available as types
-    typedef typename base_array_unmanaged<B,A>::iterator Iterator;
+    typedef typename base_array_unmanaged<B,ST>::iterator Iterator;
 
     //! make iterators available as types
-    typedef typename base_array_unmanaged<B,A>::const_iterator ConstIterator;
+    typedef typename base_array_unmanaged<B,ST>::const_iterator ConstIterator;
 
     //! for STL compatibility
     typedef B value_type;
@@ -177,8 +174,8 @@ namespace Imp {
      * @param y other (compatible)  vector
      * @return
      */
-    template<class OtherB, class OtherA>
-    auto operator* (const block_vector_unmanaged<OtherB,OtherA>& y) const
+    template<class OtherB, class OtherST>
+    auto operator* (const block_vector_unmanaged<OtherB,OtherST>& y) const
     {
       typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
       PromotedType sum(0);
@@ -198,8 +195,8 @@ namespace Imp {
      * @param y other (compatible) vector
      * @return
      */
-    template<class OtherB, class OtherA>
-    auto dot(const block_vector_unmanaged<OtherB,OtherA>& y) const
+    template<class OtherB, class OtherST>
+    auto dot(const block_vector_unmanaged<OtherB,OtherST>& y) const
     {
       typedef typename PromotionTraits<field_type,typename BlockTraits<OtherB>::field_type>::PromotedType PromotedType;
       PromotedType sum(0);
@@ -338,7 +335,7 @@ namespace Imp {
 
   protected:
     //! make constructor protected, so only derived classes can be instantiated
-    block_vector_unmanaged () : base_array_unmanaged<B,A>()
+    block_vector_unmanaged () : base_array_unmanaged<B,ST>()
     {       }
   };
 
@@ -391,7 +388,7 @@ namespace Imp {
           enables error checking.
    */
   template<class B, class A=std::allocator<B> >
-  class BlockVector : public Imp::block_vector_unmanaged<B,A>
+  class BlockVector : public Imp::block_vector_unmanaged<B,typename A::size_type>
   {
   public:
 
@@ -409,15 +406,11 @@ namespace Imp {
     //! The type for the index access
     typedef typename A::size_type size_type;
 
-    //! increment block level counter
-    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
-    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
-
     //! make iterators available as types
-    typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
+    typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
 
     //! make iterators available as types
-    typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+    typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
 
     //===== constructors and such
 
@@ -556,7 +549,7 @@ namespace Imp {
     BlockVector& operator= (const field_type& k)
     {
       // forward to operator= in base class
-      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
+      (static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
       return *this;
     }
 
@@ -623,7 +616,7 @@ namespace Imp {
 #else
   template<class B, class A=std::allocator<B> >
 #endif
-  class BlockVectorWindow : public Imp::block_vector_unmanaged<B,A>
+  class BlockVectorWindow : public Imp::block_vector_unmanaged<B,typename A::size_type>
   {
   public:
 
@@ -641,20 +634,16 @@ namespace Imp {
     //! The type for the index access
     typedef typename A::size_type size_type;
 
-    //! increment block level counter
-    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
-    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
-
     //! make iterators available as types
-    typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
+    typedef typename Imp::block_vector_unmanaged<B,size_type>::Iterator Iterator;
 
     //! make iterators available as types
-    typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+    typedef typename Imp::block_vector_unmanaged<B,size_type>::ConstIterator ConstIterator;
 
 
     //===== constructors and such
     //! makes empty array
-    BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
+    BlockVectorWindow () : Imp::block_vector_unmanaged<B,size_type>()
     {       }
 
     //! make array from given pointer and size
@@ -690,7 +679,7 @@ namespace Imp {
     //! assign from scalar
     BlockVectorWindow& operator= (const field_type& k)
     {
-      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
+      (static_cast<Imp::block_vector_unmanaged<B,size_type>&>(*this)) = k;
       return *this;
     }
 
@@ -749,8 +738,8 @@ namespace Imp {
 
    \internal This class is an implementation detail, and should not be used outside of dune-istl.
    */
-  template<class B, class A=std::allocator<B> >
-  class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,A>
+  template<class B, class ST=std::size_t >
+  class compressed_block_vector_unmanaged : public compressed_base_array_unmanaged<B,ST>
   {
   public:
 
@@ -762,17 +751,14 @@ namespace Imp {
     //! export the type representing the components
     typedef B block_type;
 
-    //! export the allocator type
-    typedef A allocator_type;
-
     //! make iterators available as types
-    typedef typename compressed_base_array_unmanaged<B,A>::iterator Iterator;
+    typedef typename compressed_base_array_unmanaged<B,ST>::iterator Iterator;
 
     //! make iterators available as types
-    typedef typename compressed_base_array_unmanaged<B,A>::const_iterator ConstIterator;
+    typedef typename compressed_base_array_unmanaged<B,ST>::const_iterator ConstIterator;
 
     //! The type for the index access
-    typedef typename A::size_type size_type;
+    typedef ST size_type;
 
     //===== assignment from scalar
 
@@ -969,7 +955,7 @@ namespace Imp {
 
   protected:
     //! make constructor protected, so only derived classes can be instantiated
-    compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,A>()
+    compressed_block_vector_unmanaged () : compressed_base_array_unmanaged<B,ST>()
     {       }
 
     //! return true if index sets coincide
@@ -1003,8 +989,8 @@ namespace Imp {
 
    \internal This class is an implementation detail, and should not be used outside of dune-istl.
    */
-  template<class B, class A=std::allocator<B> >
-  class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,A>
+  template<class B, class ST=std::size_t >
+  class CompressedBlockVectorWindow : public compressed_block_vector_unmanaged<B,ST>
   {
   public:
 
@@ -1016,26 +1002,19 @@ namespace Imp {
     //! export the type representing the components
     typedef B block_type;
 
-    //! export the allocator type
-    typedef A allocator_type;
-
     //! The type for the index access
-    typedef typename A::size_type size_type;
-
-    //! increment block level counter
-    [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
-    static constexpr unsigned int blocklevel = blockLevel<B>()+1;
+    typedef ST size_type;
 
     //! make iterators available as types
-    typedef typename compressed_block_vector_unmanaged<B,A>::Iterator Iterator;
+    typedef typename compressed_block_vector_unmanaged<B,ST>::Iterator Iterator;
 
     //! make iterators available as types
-    typedef typename compressed_block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+    typedef typename compressed_block_vector_unmanaged<B,ST>::ConstIterator ConstIterator;
 
 
     //===== constructors and such
     //! makes empty array
-    CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,A>()
+    CompressedBlockVectorWindow () : compressed_block_vector_unmanaged<B,ST>()
     {       }
 
     //! make array from given pointers and size
@@ -1074,7 +1053,7 @@ namespace Imp {
     //! assign from scalar
     CompressedBlockVectorWindow& operator= (const field_type& k)
     {
-      (static_cast<compressed_block_vector_unmanaged<B,A>&>(*this)) = k;
+      (static_cast<compressed_block_vector_unmanaged<B,ST>&>(*this)) = k;
       return *this;
     }
 

dune/istl/cholmod.hh

 // SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
 #pragma once
 
-#if HAVE_SUITESPARSE_CHOLMOD
+#if HAVE_SUITESPARSE_CHOLMOD || defined DOXYGEN
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 #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>
@@ -66,16 +68,157 @@ namespace Impl{
     });
   }
 
+  // wrapper class for C function calls to CHOLMOD itself.
+  // The CHOLMOD API has different functions for different index types.
+  template <class Index>
+  struct CholmodMethodChooser;
+
+  // specialization using 'int' to store indices
+  template <>
+  struct CholmodMethodChooser<int>
+  {
+    [[nodiscard]]
+    static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_allocate_dense(nrow,ncol,d,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_sparse* allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, int sorted, int packed, int stype, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
+    {
+      return ::cholmod_analyze(A,c);
+    }
+
+    static int defaults(cholmod_common *c)
+    {
+      return ::cholmod_defaults(c);
+    }
+
+    static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
+    {
+      return ::cholmod_factorize(A,L,c);
+    }
+
+    static int finish(cholmod_common *c)
+    {
+      return ::cholmod_finish(c);
+    }
+
+    static int free_dense (cholmod_dense **X, cholmod_common *c)
+    {
+      return ::cholmod_free_dense(X,c);
+    }
+
+    static int free_factor(cholmod_factor **L, cholmod_common *c)
+    {
+      return ::cholmod_free_factor(L,c);
+    }
+
+    static int free_sparse(cholmod_sparse **A, cholmod_common *c)
+    {
+      return ::cholmod_free_sparse(A,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
+    {
+      return ::cholmod_solve(sys,L,B,c);
+    }
+
+    static int start(cholmod_common *c)
+    {
+      return ::cholmod_start(c);
+    }
+  };
+
+  // specialization using 'SuiteSparse_long' to store indices
+  template <>
+  struct CholmodMethodChooser<SuiteSparse_long>
+  {
+    [[nodiscard]]
+    static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_l_allocate_dense(nrow,ncol,d,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_sparse* allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, int sorted, int packed, int stype, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_l_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
+    {
+      return ::cholmod_l_analyze(A,c);
+    }
+
+    static int defaults(cholmod_common *c)
+    {
+      return ::cholmod_l_defaults(c);
+    }
+
+    static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
+    {
+      return ::cholmod_l_factorize(A,L,c);
+    }
+
+    static int finish(cholmod_common *c)
+    {
+      return ::cholmod_l_finish(c);
+    }
+
+    static int free_dense (cholmod_dense **X, cholmod_common *c)
+    {
+      return ::cholmod_l_free_dense(X,c);
+    }
+
+    static int free_factor (cholmod_factor **L, cholmod_common *c)
+    {
+      return ::cholmod_l_free_factor(L,c);
+    }
+
+    static int free_sparse(cholmod_sparse **A, cholmod_common *c)
+    {
+      return ::cholmod_l_free_sparse(A,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
+    {
+      return ::cholmod_l_solve(sys,L,B,c);
+    }
+
+    static int start(cholmod_common *c)
+    {
+      return ::cholmod_l_start(c);
+    }
+  };
 
 } //namespace Impl
 
 /** @brief Dune wrapper for SuiteSparse/CHOLMOD solver
   *
   * This class implements an InverseOperator between Vector types
+  *
+  * \tparam Vector Data type for solution and right-hand-side vectors
+  * \tparam Index Type used by CHOLMOD for indices.  Must be either 'int' or 'SuiteSparse_long'
+  *   (which is usually `long int`).
   */
-template<class Vector>
+template<class Vector, class Index=int>
 class Cholmod : public InverseOperator<Vector, Vector>
 {
+  static_assert(std::is_same_v<Index,int> || std::is_same_v<Index,SuiteSparse_long>,
+                "Index type must be either 'int' or 'SuiteSparse_long'!");
+
+  using CholmodMethod = Impl::CholmodMethodChooser<Index>;
+
 public:
 
   /** @brief Default constructor
@@ -85,7 +228,7 @@ public:
    */
   Cholmod()
   {
-    cholmod_start(&c_);
+    CholmodMethod::start(&c_);
   }
 
   /** @brief Destructor
@@ -96,8 +239,8 @@ public:
   ~Cholmod()
   {
     if (L_)
-      cholmod_free_factor(&L_, &c_);
-    cholmod_finish(&c_);
+      CholmodMethod::free_factor(&L_, &c_);
+    CholmodMethod::finish(&c_);
   }
 
   // forbid copying to avoid freeing memory twice
@@ -107,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);
   }
@@ -117,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_ )
@@ -144,13 +287,14 @@ public:
     });
 
       // create a cholmod dense object
-    auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
+    auto b3 = make_cholmod_dense(CholmodMethod::allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
+
     // cast because void-ptr
     auto b4 = static_cast<double*>(b3->x);
     std::copy(b2.get(), b2.get() + L_->n, b4);
 
     // solve for a cholmod x object
-    auto x3 = make_cholmod_dense(cholmod_solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
+    auto x3 = make_cholmod_dense(CholmodMethod::solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
     // cast because void-ptr
     auto xp = static_cast<double*>(x3->x);
 
@@ -199,8 +343,8 @@ public:
   void setMatrix(const Matrix& matrix, const Ignore* ignore)
   {
     // count the number of entries and diagonal entries
-    int nonZeros = 0;
-    int numberOfIgnoredDofs = 0;
+    size_t nonZeros = 0;
+    size_t numberOfIgnoredDofs = 0;
 
 
     auto [flatRows,flatCols] = flatMatrixForEach( matrix, [&](auto&& /*entry*/, auto&& flatRowIndex, auto&& flatColIndex){
@@ -216,38 +360,38 @@ public:
       numberOfIgnoredDofs = std::count(flatIgnore.begin(),flatIgnore.end(),true);
     }
 
-    // Total number of rows
-    int N = flatRows - numberOfIgnoredDofs;
-
-    nIsZero_ = (N <= 0);
+    nIsZero_ = (size_t(flatRows) <= numberOfIgnoredDofs);
 
     if ( nIsZero_ )
     {
         return;
     }
 
+    // Total number of rows
+    size_t N = flatRows - numberOfIgnoredDofs;
+
     /*
     * CHOLMOD uses compressed-column sparse matrices, but for symmetric
     * matrices this is the same as the compressed-row sparse matrix used
     * by DUNE.  So we can just store Mᵀ instead of M (as M = Mᵀ).
     */
     const auto deleter = [c = &this->c_](auto* p) {
-      cholmod_free_sparse(&p, c);
+      CholmodMethod::free_sparse(&p, c);
     };
     auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
-      cholmod_allocate_sparse(N,             // # rows
-                              N,             // # cols
-                              nonZeros,      // # of nonzeroes
-                              1,             // indices are sorted ( 1 = true)
-                              1,             // matrix is "packed" ( 1 = true)
-                              -1,            // stype of matrix ( -1 = consider the lower part only )
-                              CHOLMOD_REAL,  // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
-                              &c_            // cholmod_common ptr
-                             ), deleter);
+      CholmodMethod::allocate_sparse(N,             // # rows
+                                     N,             // # cols
+                                     nonZeros,      // # of nonzeroes
+                                     1,             // indices are sorted ( 1 = true)
+                                     1,             // matrix is "packed" ( 1 = true)
+                                     -1,            // stype of matrix ( -1 = consider the lower part only )
+                                     CHOLMOD_REAL,  // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
+                                     &c_            // cholmod_common ptr
+      ), deleter);
 
     // copy the data of BCRS matrix to Cholmod Sparse matrix
-    int* Ap = static_cast<int*>(M->p);
-    int* Ai = static_cast<int*>(M->i);
+    Index* Ap = static_cast<Index*>(M->p);
+    Index* Ai = static_cast<Index*>(M->i);
     double* Ax = static_cast<double*>(M->x);
 
 
@@ -287,7 +431,7 @@ public:
 
     // now accumulate
     Ap[0] = 0;
-    for ( int i=0; i<N; i++ )
+    for ( size_t i=0; i<N; i++ )
     {
       Ap[i+1] += Ap[i];
     }
@@ -317,14 +461,18 @@ 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_ = cholmod_analyze(M.get(), &c_);
+    L_ = CholmodMethod::analyze(M.get(), &c_);
 
     // Do the factorization (this may take some time)
-    cholmod_factorize(M.get(), L_, &c_);
+    CholmodMethod::factorize(M.get(), L_, &c_);
   }
 
-  virtual SolverCategory::Category category() const
+  SolverCategory::Category category() const override
   {
     return SolverCategory::Category::sequential;
   }
@@ -332,7 +480,7 @@ public:
   /** \brief return a reference to the CHOLMOD common object for advanced option settings
    *
    *  The CHOLMOD common object stores all parameters and options for the solver to run
-   *  and can be modified in several ways, see CHOLMOD Userguide for further information
+   *  and can be modified in several ways, see CHOLMOD Userguide for further information.
    */
   cholmod_common& cholmodCommonObject()
   {
@@ -365,7 +513,7 @@ private:
   auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
   {
     const auto deleter = [c](auto* p) {
-      cholmod_free_dense(&p, c);
+      CholmodMethod::free_dense(&p, c);
     };
     return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
   }
@@ -380,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 */