dune/istl/common/registry.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
+#ifndef DUNE_ISTL_COMMON_REGISTRY_HH
+#define DUNE_ISTL_COMMON_REGISTRY_HH
+
+#include <cstddef>
+#include <iostream>
+#include <memory>
+#include <string>
+#include <utility>
+
+#include "counter.hh"
+
+#include <dune/common/typelist.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/parameterizedobject.hh>
+
+#define DUNE_REGISTRY_PUT(Tag, id, ...)               \
+  namespace {                                   \
+    template<>                                  \
+    struct Registry<Tag, DUNE_GET_COUNTER(Tag)>       \
+    {                                           \
+      static auto getCreator()                  \
+      {                                         \
+        return __VA_ARGS__;                     \
+      }                                         \
+      static std::string name() { return id; }  \
+    };                                          \
+  }                                             \
+  DUNE_INC_COUNTER(Tag)
+
+
+namespace Dune {
+  namespace {
+    template<class Tag, std::size_t index>
+    struct Registry;
+  }
+
+  namespace {
+    template<template<class> class Base, class V, class Tag, typename... Args>
+    auto registryGet(Tag , std::string name, Args... args)
+    {
+      constexpr auto count = DUNE_GET_COUNTER(Tag);
+      std::shared_ptr<Base<V> > result;
+      Dune::Hybrid::forEach(std::make_index_sequence<count>{},
+                            [&](auto index) {
+                              using Reg = Registry<Tag, index>;
+                              if(!result && Reg::name() == name) {
+                                result = Reg::getCreator()(Dune::MetaType<V>{}, args...);
+                              }
+                            });
+      return result;
+    }
+
+    /*
+      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.
+    */
+    template<class V, class Type, class Tag, class... Args>
+    int addRegistryToFactory(Dune::ParameterizedObjectFactory<Type(Args...), std::string>& factory,
+                              Tag){
+      constexpr auto count = DUNE_GET_COUNTER(Tag);
+      Dune::Hybrid::forEach(std::make_index_sequence<count>{},
+                            [&](auto index) {
+                              // we first get the generic lambda
+                              // and later specialize it with given parameters.
+                              // doing all at once leads to an ICE with g++-6
+                              using Reg = Registry<Tag, index>;
+                              auto genericcreator = Reg::getCreator();
+                              factory.define(Reg::name(), [genericcreator](Args... args){
+                                                            return genericcreator(V{}, args...);
+                                                          });
+                            });
+      return count;
+    }
+  } // end anonymous namespace
+} // end namespace Dune
+
+#endif // DUNE_ISTL_COMMON_REGISTRY_HH

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/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
+
 add_subdirectory(test)
 
 #install headers

dune/istl/eigenvalue/arpackpp.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_EIGENVALUE_ARPACKPP_HH
@@ -13,6 +15,7 @@
 #include <dune/common/fvector.hh>     // provides Dune::FieldVector
 #include <dune/common/exceptions.hh>  // provides DUNE_THROW(...)
 
+#include <dune/istl/blocklevel.hh>     // provides Dune::blockLevel
 #include <dune/istl/bvector.hh>        // provides Dune::BlockVector
 #include <dune/istl/istlexception.hh>  // provides Dune::ISTLError
 #include <dune/istl/io.hh>             // provides Dune::printvector(...)
@@ -64,7 +67,7 @@ namespace Dune
       {
         // assert that BCRSMatrix type has blocklevel 2
         static_assert
-          (BCRSMatrix::blocklevel == 2,
+          (blockLevel<BCRSMatrix>() == 2,
             "Only BCRSMatrices with blocklevel 2 are supported.");
 
         // allocate memory for auxiliary block vector objects
@@ -201,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/"
@@ -232,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.
    */
@@ -253,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,
@@ -593,15 +596,15 @@ namespace Dune
     /**
      * \brief Assume the matrix to be nonsymmetric and perform IRLM
      *        to compute an approximation sigma of its largest
-     *        singlar value and the corresponding approximation x of
+     *        singular value and the corresponding approximation x of
      *        an associated singular vector.
      *
      * \param[in]  epsilon The target relative accuracy of Ritz values
      *                     (0 == machine precision).
-     * \param[out] sigma   The approximated largest singlar value.
+     * \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
@@ -705,15 +708,15 @@ namespace Dune
     /**
      * \brief Assume the matrix to be nonsymmetric and perform IRLM
      *        to compute an approximation sigma of its smallest
-     *        singlar value and the corresponding approximation x of
+     *        singular value and the corresponding approximation x of
      *        an associated singular vector.
      *
      * \param[in]  epsilon The target relative accuracy of Ritz values
      *                     (0 == machine precision).
-     * \param[out] sigma   The approximated smallest singlar value.
+     * \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/poweriteration.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_EIGENVALUE_POWERITERATION_HH
@@ -16,6 +18,7 @@
 
 #include <dune/common/exceptions.hh>  // provides DUNE_THROW(...)
 
+#include <dune/istl/blocklevel.hh>      // provides Dune::blockLevel
 #include <dune/istl/operators.hh>       // provides Dune::LinearOperator
 #include <dune/istl/solvercategory.hh>  // provides Dune::SolverCategory::sequential
 #include <dune/istl/solvertype.hh>      // provides Dune::IsDirectSolver
@@ -53,13 +56,13 @@ namespace Dune
           mutable_scaling_(mutable_scaling)
       {}
 
-      virtual void apply (const X& x, Y& y) const
+      void apply (const X& x, Y& y) const override
       {
         y = x;
         y *= immutable_scaling_*mutable_scaling_;
       }
 
-      virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+      void applyscaleadd (field_type alpha, const X& x, Y& y) const override
       {
         X temp(x);
         temp *= immutable_scaling_*mutable_scaling_;
@@ -67,7 +70,7 @@ namespace Dune
       }
 
       //! Category of the linear operator (see SolverCategory::Category)
-      virtual SolverCategory::Category category() const
+      SolverCategory::Category category() const override
       {
         return SolverCategory::sequential;
       }
@@ -105,14 +108,14 @@ namespace Dune
           "Range type of both operators doesn't match!");
       }
 
-      virtual void apply (const domain_type& x, range_type& y) const
+      void apply (const domain_type& x, range_type& y) const override
       {
         op1_.apply(x,y);
         op2_.applyscaleadd(1.0,x,y);
       }
 
-      virtual void applyscaleadd (field_type alpha,
-        const domain_type& x, range_type& y) const
+      void applyscaleadd (field_type alpha,
+        const domain_type& x, range_type& y) const override
       {
         range_type temp(y);
         op1_.apply(x,temp);
@@ -121,7 +124,7 @@ namespace Dune
       }
 
       //! Category of the linear operator (see SolverCategory::Category)
-      virtual SolverCategory::Category category() const
+      SolverCategory::Category category() const override
       {
         return SolverCategory::sequential;
       }
@@ -215,7 +218,7 @@ namespace Dune
     {
       // assert that BCRSMatrix type has blocklevel 2
       static_assert
-        (BCRSMatrix::blocklevel == 2,
+        (blockLevel<BCRSMatrix>() == 2,
          "Only BCRSMatrices with blocklevel 2 are supported.");
 
       // assert that BCRSMatrix type has square blocks