.gitlab-ci.yml

@@ -22,33 +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_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]
-
 # Check for spelling mistakes in text
 code-spelling-check:
   stage: .pre

CHANGELOG.md

@@ -3,8 +3,43 @@ 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.
+
+- 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
 
@@ -21,6 +56,9 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 - 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
 
@@ -48,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*`.
 

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,20 @@ 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)
+
 # 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 +45,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--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,15 +45,17 @@ 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
-2004--2021    Oliver Sander
+2023          Vaishnavi Sanchi
+2004--2024    Oliver Sander
 2015--2019    Linus Seelinger
 2013          Bård Skaflestad
 2018          Mathis Springwald
+2023          Jakob Torben
 2006--2010    Martin Weiser
 2019--2020    Kilian Weishaupt
 2015          Sebastian Westerheide

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

@@ -493,12 +493,12 @@ 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;
 
@@ -775,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)
@@ -1218,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!
      */
@@ -2199,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/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:
 
@@ -410,10 +407,10 @@ namespace Imp {
     typedef typename A::size_type size_type;
 
     //! 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
 
@@ -552,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;
     }
 
@@ -619,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:
 
@@ -638,15 +635,15 @@ namespace Imp {
     typedef typename A::size_type size_type;
 
     //! 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
@@ -682,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;
     }
 
@@ -741,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:
 
@@ -754,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
 
@@ -961,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
@@ -995,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:
 
@@ -1008,22 +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;
+    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
@@ -1062,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 */
 

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

+// 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:
+#ifndef DUNE_ISTL_DILU_HH
+#define DUNE_ISTL_DILU_HH
+
+#include <cmath>
+#include <complex>
+#include <map>
+#include <vector>
+#include <sstream>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/scalarvectorview.hh>
+#include <dune/common/scalarmatrixview.hh>
+
+#include "istlexception.hh"
+
+/** \file
+ * \brief  The diagonal incomplete LU factorization kernels
+ */
+
+namespace Dune
+{
+
+  /** @addtogroup ISTL_Kernel
+          @{
+   */
+
+  namespace DILU
+  {
+
+    /*! compute DILU decomposition of A
+
+        The preconditioner matrix M has the property
+
+        diag(A) = diag(M) = diag((D + L_A) D^{-1} (D + U_A)) = diag(D + L_A D^{-1} U_A)
+
+        such that the diagonal matrix D can be constructed:
+
+        D_11 = A_11
+        D_22 = A22 - L_A_{21} D_{11}^{-1} U_A_{12}
+        and etc
+
+        we store the inverse of D to be used when applying the preconditioner.
+
+        For more details, see: R. Barrett et al., "Templates for the Solution of Linear Systems:
+        Building Blocks for Iterative Methods", 1994. Available from: https://www.netlib.org/templates/templates.pdf
+     */
+    template <class M>
+    void blockDILUDecomposition(M &A, std::vector<typename M::block_type> &Dinv_)
+    {
+      auto endi = A.end();
+      for (auto row = A.begin(); row != endi; ++row)
+      {
+        const auto row_i = row.index();
+        const auto col = row->find(row_i);
+        // initialise Dinv[i] = A[i, i]
+        Dinv_[row_i] = *col;
+      }
+
+      for (auto row = A.begin(); row != endi; ++row)
+      {
+        const auto row_i = row.index();
+        for (auto a_ij = row->begin(); a_ij.index() < row_i; ++a_ij)
+        {
+          const auto col_j = a_ij.index();
+          const auto a_ji = A[col_j].find(row_i);
+          // if A[i, j] != 0 and A[j, i] != 0
+          if (a_ji != A[col_j].end())
+          {
+            // Dinv[i] -= A[i, j] * d[j] * A[j, i]
+            Dinv_[row_i] -= (*a_ij) * Dinv_[col_j] * (*a_ji);
+          }
+        }
+
+        // store the inverse
+        try
+        {
+          Impl::asMatrix(Dinv_[row_i]).invert(); // compute inverse of diagonal block
+        }
+        catch (Dune::FMatrixError &e)
+        {
+          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;
+        }
+      }
+    }
+
+    /*! DILU backsolve
+
+      M = (D + L_A) D^-1 (D + U_A)   (a LU decomposition of M)
+      where L_A and U_A are the strictly lower and upper parts of A.
+
+      Working with residual d = b - Ax and update v = x_{n+1} - x_n
+      solving the product M^-1(d) using upper and lower triangular solve:
+
+      v = M^{-1} d = (D + U_A)^{-1} D (D + L_A)^{-1}  d
+
+      define y = (D + L_A)^{-1} d
+
+      lower triangular solve: (D + L_A) y = d
+      upper triangular solve: (D + U_A) v = D y
+     */
+    template <class M, class X, class Y>
+    void blockDILUBacksolve(const M &A, const std::vector<typename M::block_type> Dinv_, X &v, const Y &d)
+    {
+      using dblock = typename Y::block_type;
+      using vblock = typename X::block_type;
+
+      // lower triangular solve: (D + L_A) y = d
+      auto endi = A.end();
+      for (auto row = A.begin(); row != endi; ++row)
+      {
+        const auto row_i = row.index();
+        dblock rhsValue(d[row_i]);
+        auto &&rhs = Impl::asVector(rhsValue);
+        for (auto a_ij = (*row).begin(); a_ij.index() < row_i; ++a_ij)
+        {
+          // if  A[i][j] != 0
+          // rhs -= A[i][j]* y[j], where v_j stores y_j
+          const auto col_j = a_ij.index();
+          Impl::asMatrix(*a_ij).mmv(v[col_j], rhs);
+        }
+        // y_i = Dinv_i * rhs
+        // storing y_i in v_i
+        auto &&vi = Impl::asVector(v[row_i]);
+        Impl::asMatrix(Dinv_[row_i]).mv(rhs, vi); // (D + L_A)_ii = D_i
+      }
+
+      // upper triangular solve: (D + U_A) v = Dy
+      auto rendi = A.beforeBegin();
+      for (auto row = A.beforeEnd(); row != rendi; --row)
+      {
+        const auto row_i = row.index();
+        // rhs = 0
+        vblock rhs(0.0);
+        for (auto a_ij = (*row).beforeEnd(); a_ij.index() > row_i; --a_ij)
+        {
+          // if A[i][j] != 0
+          // rhs += A[i][j]*v[j]
+          const auto col_j = a_ij.index();
+          Impl::asMatrix(*a_ij).umv(v[col_j], rhs);
+        }
+        // calculate update v = M^-1*d
+        // v_i = y_i - Dinv_i*rhs
+        // before update v_i is y_i
+        auto &&vi = Impl::asVector(v[row_i]);
+        Impl::asMatrix(Dinv_[row_i]).mmv(rhs, vi);
+      }
+    }
+  } // end namespace DILU
+
+  /** @} end documentation */
+
+} // end namespace
+
+#endif

dune/istl/eigenvalue/arpackpp.hh

@@ -204,17 +204,17 @@ namespace Dune
   /**
    * \brief Wrapper to use a range of ARPACK++ eigenvalue solvers
    *
-   *        A class template for performing some eigenvalue algorithms
-   *        provided by the ARPACK++ library which is based on the implicitly
-   *        restarted Arnoldi/Lanczos method (IRAM/IRLM), a synthesis of the
-   *        Arnoldi/Lanczos process with the implicitily shifted QR technique.
-   *        The method is designed to compute eigenvalue-eigenvector pairs of
-   *        large scale sparse nonsymmetric/symmetric matrices. This class
-   *        template uses the algorithms to compute the dominant (i.e. largest
-   *        magnitude) and least dominant (i.e. smallest magnitude) eigenvalue
-   *        as well as the spectral condition number of square, symmetric
-   *        matrices and to compute the largest and smallest singular value as
-   *        well as the spectral condition number of nonsymmetric matrices.
+   * A class template for performing some eigenvalue algorithms
+   * provided by the ARPACK++ library which is based on the implicitly
+   * restarted Arnoldi/Lanczos method (IRAM/IRLM), a synthesis of the
+   * Arnoldi/Lanczos process with the implicitily shifted QR technique.
+   * The method is designed to compute eigenvalue-eigenvector pairs of
+   * large scale sparse nonsymmetric/symmetric matrices. This class
+   * template uses the algorithms to compute the dominant (i.e. largest
+   * magnitude) and least dominant (i.e. smallest magnitude) eigenvalue
+   * as well as the spectral condition number of square, symmetric
+   * matrices and to compute the largest and smallest singular value as
+   * well as the spectral condition number of nonsymmetric matrices.
    *
    * \note For a recent version of the ARPACK++ library working with recent
    *       compiler versions see "http://reuter.mit.edu/software/arpackpatch/"
@@ -235,8 +235,8 @@ namespace Dune
    *                     respectively singular values shall be considered;
    *                     is assumed to have blocklevel 2.
    * \tparam BlockVector Type of the associated vectors; compatible with the
-   *                     rows of a BCRSMatrix object (if #rows >= #ncols) or
-   *                     its columns (if #rows < #ncols).
+   *                     rows of a BCRSMatrix object (if \#rows \>= \#ncols) or
+   *                     its columns (if \#rows \< \#ncols).
    *
    * \author Sebastian Westerheide.
    */
@@ -256,14 +256,14 @@ namespace Dune
      *                            update iterations allowed; depending on the
      *                            algorithm, c*nIterationsMax iterations may
      *                            be performed, where c is a natural number.
-     * \param[in] verbosity_level Verbosity setting;
-     *                            >= 1: algorithms print a preamble and
-     *                                  the final result,
-     *                            >= 2: algorithms print information about
-     *                                  the problem solved using ARPACK++,
-     *                            >= 3: the final result output includes
-     *                                  the approximated eigenvector,
-     *                            >= 4: sets the ARPACK(++) verbosity mode.
+     * \param[in] verbosity_level Verbosity setting;<br>
+     *                            \>= 1: algorithms print a preamble and
+     *                                  the final result,<br>
+     *                            \>= 2: algorithms print information about
+     *                                  the problem solved using ARPACK++,<br>
+     *                            \>= 3: the final result output includes
+     *                                  the approximated eigenvector,<br>
+     *                            \>= 4: sets the ARPACK(++) verbosity mode.
      */
     ArPackPlusPlus_Algorithms (const BCRSMatrix& m,
                                const unsigned int nIterationsMax = 100000,
@@ -603,8 +603,8 @@ namespace Dune
      *                     (0 == machine precision).
      * \param[out] sigma   The approximated largest singular value.
      * \param[out] x       The associated approximated right-singular
-     *                     vector (if #rows >= #ncols) respectively
-     *                     left-singular vector (if #rows < #ncols).
+     *                     vector (if \#rows \>= \#ncols) respectively
+     *                     left-singular vector (if \#rows \< \#ncols).
      */
     inline void computeNonSymMax (const Real& epsilon,
                                   BlockVector& x, Real& sigma) const
@@ -715,8 +715,8 @@ namespace Dune
      *                     (0 == machine precision).
      * \param[out] sigma   The approximated smallest singular value.
      * \param[out] x       The associated approximated right-singular
-     *                     vector (if #rows >= #ncols) respectively
-     *                     left-singular vector (if #rows < #ncols).
+     *                     vector (if \#rows \>= \#ncols) respectively
+     *                     left-singular vector (if \#rows \< \#ncols).
      */
     inline void computeNonSymMin (const Real& epsilon,
                                   BlockVector& x, Real& sigma) const