.gitlab-ci.yml

@@ -31,25 +31,6 @@ debian-11-gcc-9-17-with-checking:
     DUNECI_CMAKE_FLAGS: "-DDUNE_MAX_TEST_CORES=4 -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE"
   tags: [duneci]
 
-# 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

@@ -5,6 +5,16 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 # Master (will become release 2.10)
 
+- `setupLaplacian` has now 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.
@@ -15,6 +25,17 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
   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`,
@@ -30,6 +51,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
 

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
 
-cmake_minimum_required(VERSION 3.13)
+cmake_minimum_required(VERSION 3.16)
 project(dune-istl C CXX)
 
 # guess build tree of dune-common
@@ -39,4 +39,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

@@ -46,10 +46,12 @@ Copyright holders:
 2018--2022    Lukas Renelt
 2021--2022    Santiago Ospina De Los Ríos
 2008--2013    Uli Sack
+2023          Vaishnavi Sanchi
 2004--2021    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

@@ -17,7 +17,7 @@ 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

@@ -19,7 +19,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 +28,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/FindARPACKPP.cmake

@@ -123,12 +123,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()

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/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,7 +5,7 @@
 #ifndef DUNE_ISTL_BASEARRAY_HH
 #define DUNE_ISTL_BASEARRAY_HH
 
-#include "assert.h"
+#include <cassert>
 #include <cmath>
 #include <cstddef>
 #include <memory>
@@ -41,13 +41,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 +53,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&;
@@ -302,7 +295,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 +305,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&;

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/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 <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)
+        {
+          DUNE_THROW(MatrixBlockError, "DILU failed to invert matrix block D[" << row_i << "]"
+                                                                               << e.what();
+                     th__ex.r = row_i;);
+        }
+      }
+    }
+
+    /*! 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

dune/istl/eigenvalue/test/cond2test.cc

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

dune/istl/foreach.hh

@@ -167,8 +167,12 @@ std::pair<std::size_t,std::size_t> flatMatrixForEach(Matrix&& matrix, F&& f, std
         {
           auto&& entry = *colIt;
           auto colIdx = colIt.index();
-          auto [ dummyRows, dummyCols ] = flatMatrixForEach(entry, f, rowOffset + rowIdx*blockRows, colOffset + colIdx*blockCols);
-          assert( dummyRows == blockRows and dummyCols == blockCols and "we need the same size of each block in this matrix type");
+#ifndef NDEBUG
+          // only instantiate return value in debug mode (for the assert)
+          auto [ numRows, numCols ] =
+#endif
+          flatMatrixForEach(entry, f, rowOffset + rowIdx*blockRows, colOffset + colIdx*blockCols);
+          assert( numRows == blockRows and numCols == blockCols and "we need the same size of each block in this matrix type");
         }
       }
 

dune/istl/io.hh

@@ -246,6 +246,36 @@ namespace Dune {
     s.precision(oldprec);
   }
 
+  namespace Impl
+  {
+    template<class B>
+    void printInnerMatrixElement(std::ostream& s,
+                                 const B& innerMatrixElement,
+                                 int innerrow, int innercol)
+    {
+      s<<innerMatrixElement<<" ";
+    }
+
+    template<class B, int n>
+    void printInnerMatrixElement(std::ostream& s,
+                                 const ScaledIdentityMatrix<B,n> innerMatrixElement,
+                                 int innerrow, int innercol)
+    {
+      if (innerrow == innercol)
+        s<<innerMatrixElement.scalar()<<" ";
+      else
+        s<<"-";
+    }
+
+    template<class B, int n, int m>
+    void printInnerMatrixElement(std::ostream& s,
+                                 const FieldMatrix<B,n,m> innerMatrixElement,
+                                 int innerrow, int innercol)
+    {
+      s<<innerMatrixElement[innerrow][innercol]<<" ";
+    }
+  }
+
   /**
    * \brief Prints a BCRSMatrix with fixed sized blocks.
    *
@@ -267,13 +297,13 @@ namespace Dune {
    * @param width The number of nonzero blocks to print in one line.
    * @param precision The precision to use when printing the numbers.
    */
-  template<class B, int n, int m, class A>
+  template<class A, class InnerMatrixType>
   void printSparseMatrix(std::ostream& s,
-                         const BCRSMatrix<FieldMatrix<B,n,m>,A>& mat,
+                         const BCRSMatrix<InnerMatrixType,A>& mat,
                          std::string title, std::string rowtext,
                          int width=3, int precision=2)
   {
-    typedef BCRSMatrix<FieldMatrix<B,n,m>,A> Matrix;
+    typedef BCRSMatrix<InnerMatrixType,A> Matrix;
     // remember old flags
     std::ios_base::fmtflags oldflags = s.flags();
     // set the output format
@@ -290,6 +320,8 @@ namespace Dune {
 
     typedef typename Matrix::ConstRowIterator Row;
 
+    constexpr int n = std::decay_t<decltype(Impl::asMatrix(std::declval<InnerMatrixType>()))>::rows;
+    constexpr int m = std::decay_t<decltype(Impl::asMatrix(std::declval<InnerMatrixType>()))>::cols;
     for(Row row=mat.begin(); row != mat.end(); ++row) {
       int skipcols=0;
       bool reachedEnd=false;
@@ -323,7 +355,7 @@ namespace Dune {
             }
             for(int innercol=0; innercol < m; ++innercol) {
               s.width(9);
-              s<<(*col)[innerrow][innercol]<<" ";
+              Impl::printInnerMatrixElement(s,*col,innerrow,innercol);
             }
 
             s<<"|";
@@ -540,10 +572,10 @@ namespace Dune {
     }
 
     //! Writes (recursively) the svg content for an sparse matrix
-    template <class Mat, class Stream, class SVGMatrixOptions,
+    template <class Stream, class Mat, class SVGMatrixOptions,
               class RowPrefix, class ColPrefix>
     std::pair<std::size_t, size_t>
-    writeSVGMatrix(const Mat &mat, Stream &out, SVGMatrixOptions opts,
+    writeSVGMatrix(Stream &out, const Mat &mat, SVGMatrixOptions opts,
                   RowPrefix row_prefix, ColPrefix col_prefix) {
       // get values to fill the offsets
       const auto& block_size = opts.block_size;
@@ -599,7 +631,7 @@ namespace Dune {
           NullStream dev0;
           // get size of sub-block
           auto sub_size =
-              writeSVGMatrix(val, dev0, opts, row_prefix, col_prefix);
+              writeSVGMatrix(dev0, val, opts, row_prefix, col_prefix);
 
           // if we didn't see col size before
           if (col_offsets[col + 1] == null_offset) // write it in the offset vector
@@ -636,7 +668,7 @@ namespace Dune {
           ss << "<svg x='" << col_offsets[col] << "' y='" << row_offsets[row]
             << "' width='" << width << "' height='" << height << "'>\n";
           // write a nested svg with the contents of the sub-block
-          writeSVGMatrix(val, ss, opts, row_prefix, col_prefix);
+          writeSVGMatrix(ss, val, opts, row_prefix, col_prefix);
           ss << "</svg>\n";
         });
         col_offset = col_offsets.back();
@@ -804,12 +836,34 @@ namespace Dune {
    * @param opts          SVG Options object
    */
   template <class Mat, class SVGOptions = DefaultSVGMatrixOptions>
-  void writeSVGMatrix(const Mat &mat, std::ostream &out,
-                      SVGOptions opts = {}) {
+  void writeSVGMatrix(std::ostream &out, const Mat &mat, SVGOptions opts = {}) {
     // We need a vector that can fit all the multi-indices for rows and columns
     using IndexPrefix = Dune::ReservedVector<std::size_t, blockLevel<Mat>()>;
     // Call overload for Mat type
-    Impl::writeSVGMatrix(mat, out, opts, IndexPrefix{}, IndexPrefix{});
+    Impl::writeSVGMatrix(out, mat, opts, IndexPrefix{}, IndexPrefix{});
+  }
+
+  /**
+   * @brief Writes the visualization of matrix in the SVG format
+   * @details The default visualization writes a rectangle for every block
+   *          bounding box. This is enough to visualize patterns. If you need a
+   *          more advance SVG object in each matrix block (e.g. color scale, or
+   *          write the value in text), just provide a custom SVGOptions that
+   *          fulfills the DefaultSVGMatrixOptions interface.
+   *
+   * @tparam Mat          Matrix type to write
+   * @tparam SVGOptions   Options object type (see DefaultSVGMatrixOptions)
+   * @param mat           The matrix to write
+   * @param out           A output stream to write SVG to
+   * @param opts          SVG Options object
+   *
+   * \deprecated Use signature where std::stream is the first argument. This
+   *             method will be removed after Dune 2.10.
+   */
+  template <class Mat, class SVGOptions = DefaultSVGMatrixOptions>
+  [[deprecated("Use signature where std::stream is the first argument. This will be removed after Dune 2.10.")]]
+  void writeSVGMatrix(const Mat &mat, std::ostream &out, SVGOptions opts = {}) {
+    writeSVGMatrix(out, mat, opts);
   }
 
   /** @} end documentation */

dune/istl/matrix.hh

@@ -37,7 +37,7 @@ namespace MatrixImp
        elegant solution.
    */
   template<class B, class A=std::allocator<B> >
-  class DenseMatrixBase : public Imp::block_vector_unmanaged<B,A>
+  class DenseMatrixBase : public Imp::block_vector_unmanaged<B,typename A::size_type>
                               // this derivation gives us all the blas level 1 and norms
                               // on the large array. However, access operators have to be
                               // overwritten.
@@ -79,7 +79,7 @@ namespace MatrixImp
     /** constructor without arguments makes empty vector,
             object cannot be used yet
      */
-    DenseMatrixBase () : Imp::block_vector_unmanaged<B,A>()
+    DenseMatrixBase () : Imp::block_vector_unmanaged<B,size_type>()
     {
       // nothing is known ...
       rows_ = 0;
@@ -89,10 +89,10 @@ namespace MatrixImp
     /** make vector with given number of blocks each having a constant size,
             object is fully usable then.
 
-            \param _nblocks Number of blocks
-            \param m Number of elements in each block
+     * \param rows    Number of rows
+     * \param columns Number of columns
      */
-    DenseMatrixBase (size_type rows, size_type columns) : Imp::block_vector_unmanaged<B,A>()
+    DenseMatrixBase (size_type rows, size_type columns) : Imp::block_vector_unmanaged<B,size_type>()
     {
       // and we can allocate the big array in the base class
       this->n = rows*columns;
@@ -228,7 +228,7 @@ namespace MatrixImp
     //! assign from scalar
     DenseMatrixBase& 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;
     }
 

dune/istl/matrixindexset.hh

@@ -142,6 +142,22 @@ namespace Dune {
     /** \brief Return the number of rows */
     size_type rows() const {return rows_;}
 
+    /** \brief Return the number of columns */
+    size_type cols() const {return cols_;}
+
+    /**
+     * \brief Return column indices of entries in given row
+     *
+     * This returns a range of all column indices
+     * that have been added for the given column.
+     * Since there are different internal implementations
+     * of this range, the result is stored in a std::variant<...>
+     * which has to be accessed using `std::visit`.
+     */
+    const auto& columnIndices(size_type row) const {
+      return indices_[row];
+    }
+
     /** \brief Return the number of entries in a given row */
     size_type rowsize(size_type row) const {
       return std::visit([&](const auto& rowIndices) {
@@ -193,8 +209,7 @@ namespace Dune {
 
       for (size_type row=0; row<rows_; row++) {
         std::visit([&](const auto& rowIndices) {
-          for(size_type col : rowIndices)
-             matrix.addindex(row, col);
+          matrix.setIndicesNoSort(row, rowIndices.begin(), rowIndices.end());
         }, indices_[row]);
       }