Add DILU preconditioner

LICENSE.md

@@ -50,6 +50,7 @@ Copyright holders:
 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

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/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/preconditioners.hh

@@ -22,6 +22,7 @@
 #include "istlexception.hh"
 #include "matrixutils.hh"
 #include "gsetc.hh"
+#include "dilu.hh"
 #include "ildl.hh"
 #include "ilu.hh"
 
@@ -515,7 +516,171 @@ namespace Dune {
   };
   DUNE_REGISTER_PRECONDITIONER("jac", defaultPreconditionerBlockLevelCreator<Dune::SeqJac>());
 
+  /*!
+     \brief Sequential DILU preconditioner.
+
+     Wraps the naked ISTL generic DILU preconditioner into the solver framework.
+
+     The Diagonal Incomplete LU factorization (DILU) is a simplified variant of the ILU(0) factorization,
+     where only the diagonal elements are factorised. This preconditioner can be written as
+
+     M = (D + L_A) D^{-1} (D + U_A),
+
+     where L_A and U_A are the strictly lower and upper parts of A and D is the diagonal matrix containing
+     the generated pivot values. The 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.
+
+     Note that when applying the preconditioner, M is never explicitly created but instead a lower and
+     upper triangualr solve is performed using the values of A and D^-1. When applying the preconditioner,
+     we are working with the residual d = b - Ax and update v = x_{n+1} - x_n, such that:
+
+     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
+
+     This means that the preconditioner only requires an additional storage of the diagonal matrix D^{-1}
+
+    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
+
+     \tparam M The matrix type to operate on
+     \tparam X Type of the update
+     \tparam Y Type of the defect
+     \tparam l Ignored. Just there to have the same number of template arguments
+     as other preconditioners.
+   */
+  template <class M, class X, class Y, int l = 1>
+  class SeqDILU : public Preconditioner<X, Y>
+  {
+  public:
+    //! \brief The matrix type the preconditioner is for.
+    using matrix_type = M;
+    //! block type of matrix
+    using block_type = typename matrix_type::block_type;
+    //! \brief The domain type of the preconditioner.
+    using domain_type = X;
+    //! \brief The range type of the preconditioner.
+    using range_type = Y;
+
+    //! \brief The field type of the preconditioner.
+    using field_type = typename X::field_type;
+
+    //! \brief scalar type underlying the field_type
+    using scalar_field_type = Simd::Scalar<field_type>;
+    //! \brief real scalar type underlying the field_type
+    using real_field_type = typename FieldTraits<scalar_field_type>::real_type;
+
+    /*! \brief Constructor.
+
+       Constructor invoking DILU gets all parameters to operate the prec.
+       \param A The matrix to operate on.
+       \param w The relaxation factor.
+     */
+    SeqDILU(const M &A, real_field_type w)
+        : _A_(A),
+          _w(w),
+          wNotIdentity_([w]
+                        {using std::abs; return abs(w - real_field_type(1)) > 1e-15; }())
+    {
+      Dinv_.resize(_A_.N());
+      CheckIfDiagonalPresent<M, l>::check(_A_);
+      DILU::blockDILUDecomposition(_A_, Dinv_);
+    }
+
+    /*!
+      \brief Constructor.
+
+      \param A The assembled linear operator to use.
+      \param configuration ParameterTree containing preconditioner parameters.
+
+      ParameterTree Key | Meaning
+      ------------------|------------
+      relaxation        | The relaxation factor. default=1.0
+
+      See \ref ISTL_Factory for the ParameterTree layout and examples.
+    */
+    SeqDILU(const std::shared_ptr<const AssembledLinearOperator<M, X, Y>> &A, const ParameterTree &configuration)
+        : SeqDILU(A->getmat(), configuration)
+    {
+    }
 
+    /*!
+       \brief Constructor.
+
+       \param A The matrix to operate on.
+       \param configuration ParameterTree containing preconditioner parameters.
+
+       ParameterTree Key | Meaning
+       ------------------|------------
+      relaxation        | The relaxation factor. default=1.0
+
+       See \ref ISTL_Factory for the ParameterTree layout and examples.
+     */
+    SeqDILU(const M &A, const ParameterTree &config)
+        : SeqDILU(A, config.get<real_field_type>("relaxation", 1.0))
+    {
+    }
+
+    /*!
+       \brief Prepare the preconditioner.
+
+       \copydoc Preconditioner::pre(X&,Y&)
+     */
+    virtual void pre([[maybe_unused]] X &x, [[maybe_unused]] Y &b)
+    {
+    }
+
+    /*!
+       \brief Apply the preconditioner.
+
+       \copydoc Preconditioner::apply(X&,const Y&)
+     */
+    virtual void apply(X &v, const Y &d)
+    {
+
+      DILU::blockDILUBacksolve(_A_, Dinv_, v, d);
+
+      if (wNotIdentity_)
+      {
+        v *= _w;
+      }
+    }
+
+    /*!
+       \brief Clean up.
+
+       \copydoc Preconditioner::post(X&)
+     */
+    virtual void post([[maybe_unused]] X &x)
+    {
+    }
+
+    //! Category of the preconditioner (see SolverCategory::Category)
+    virtual SolverCategory::Category category() const
+    {
+      return SolverCategory::sequential;
+    }
+
+  protected:
+    std::vector<block_type> Dinv_;
+    //! \brief The matrix we operate on.
+    const M &_A_;
+    //! \brief The relaxation factor to use.
+    const real_field_type _w;
+    //! \brief true if w != 1.0
+    const bool wNotIdentity_;
+  };
+  DUNE_REGISTER_PRECONDITIONER("dilu", defaultPreconditionerBlockLevelCreator<Dune::SeqDILU>());
 
   /*!
      \brief Sequential ILU preconditioner.

dune/istl/test/CMakeLists.txt

@@ -17,6 +17,8 @@ dune_add_test(SOURCES bcrsnormtest.cc)
 
 dune_add_test(SOURCES cgconditiontest.cc)
 
+dune_add_test(SOURCES dilutest.cc)
+
 dune_add_test(SOURCES dotproducttest.cc)
 
 dune_add_test(SOURCES complexmatrixtest.cc)

dune/istl/test/dilutest.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
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+/** \file \brief Tests the DILU preconditioner.
+ */
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/test/testsuite.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/operators.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/solvers.hh>
+
+
+
+void seqDILUApplyIsCorrect1()
+{
+    /*
+        Tests that applying the dilu preconditioner matches the expected result
+        for a 2x2 matrix, with block size 2x2.
+                 A               x     =     b
+        | | 3  1| | 1  0| |   | |1| |     | |2| |
+        | | 2  1| | 0  1| |   | |2| |     | |1| |
+        |                 | x |     |  =  |     |
+        | | 0  0| |-1  0| |   | |1| |     | |3| |
+        | | 0  0| | 0 -1| |   | |1| |     | |4| |
+    */
+
+
+    const int N = 2;
+    const int bz = 2;
+    const int nonZeroes = 3;
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+
+    using ft = typename Vector::field_type;
+    const ft myEps((ft)1e-6);
+    Dune::TestSuite t;
+
+    Matrix A(N, N, nonZeroes, Matrix::row_wise);
+    for (auto row = A.createbegin(); row != A.createend(); ++row) {
+        row.insert(row.index());
+        if (row.index() == 0) {
+            row.insert(row.index() + 1);
+        }
+    }
+
+    A[0][0][0][0]=3.0;
+    A[0][0][0][1]=1.0;
+    A[0][0][1][0]=2.0;
+    A[0][0][1][1]=1.0;
+
+    A[0][1][0][0]=1.0;
+    A[0][1][1][1]=1.0;
+
+    A[1][1][0][0]=-1.0;
+    A[1][1][1][1]=-1.0;
+
+
+    Vector x(2);
+    x[0][0] = 1.0;
+    x[0][1] = 2.0;
+    x[1][0] = 1.0;
+    x[1][1] = 1.0;
+
+    Vector b(2);
+    b[0][0] = 2.0;
+    b[0][1] = 1.0;
+    b[1][0] = 3.0;
+    b[1][1] = 4.0;
+
+
+    auto D_00 = A[0][0];
+    auto D_00_inv = D_00;
+    D_00_inv.invert();
+    // D_11 = A_11 - L_10 D_0_inv U_01 = A_11
+    auto D_11 = A[1][1];
+    auto D_11_inv = D_11;
+    D_11_inv.invert();
+
+    // define: z = M^-1(b - Ax)
+    // where: M = (D + L_A) A^-1 (D + U_A)
+    // lower triangular solve: (D + L) y = b - Ax
+    // y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)]
+    Vector y(2);
+
+    auto rhs = b[0];
+    A[0][0].mmv(x[0], rhs);
+    A[0][1].mmv(x[1], rhs);
+    D_00_inv.mv(rhs, y[0]);
+
+    // y_1 = D_11_inv*(b_1 - (A_10*x_0 + A_11*x_1) - A_10*y_0) = D_11_inv*(b_1 - A_11*x_1)
+    rhs = b[1];
+    A[1][1].mmv(x[1], rhs);
+    D_11_inv.mv(rhs, y[1]);
+
+
+    // upper triangular solve: (D + U) z = Dy
+    // z_1 = y_1
+    Vector z(2);
+    z[1] = y[1];
+
+    // z_0 = y_0 - D_00_inv*A_01*z_1
+    z[0] = y[0];
+    auto temp = D_00_inv*A[0][1];
+    temp.mmv(z[1], z[0]);
+
+    // z is now the update x_k+1 - x_k
+
+    Dune::SeqDILU<Matrix,Vector,Vector> seqDILU(A, 1);
+    Vector v(2);
+    Vector d(b);
+    A.mmv(x, d);
+    seqDILU.apply(v, d);
+
+    // compare calculated update z from equations to update v from DILU preconditioner
+    for (int i = 0; i < 2; ++i) {
+        for (int j = 0; j < 2; ++j) {
+            t.check(std::abs(v[i][j] - z[i][j]) <= myEps)
+             << " Error in SeqDILUApplyIsCorrect1, v[" << i <<"][" << j << "]=" << v[i][j] << " != z[" << i <<"][" << j << "]=" << z[i][j];
+        }
+    }
+}
+
+
+void seqDILUApplyIsCorrect2()
+{
+    /*
+        Tests that applying the DILU preconditioner matches the expected result
+        for a 3x3 matrix, with block size 3x3.
+
+                          A                       x    =     b
+        | | 3  1  2| | 0  0  0| | 0  0  0| |   | |1| |    | |2| |
+        | | 2  3  1| | 0  0  0| | 0  0  0| |   | |2| |    | |1| |
+        | | 2  1  0| | 0  0  0| | 0  0  0| |   | |3| |    | |2| |
+        |                                  |   |     |    |     |
+        | | 0  0  0| | 1  0  1| | 1  0  2| |   | |1| |    | |2| |
+        | | 0  0  0| | 4  1  0| | 0  1  1| | x | |3| |  = | |3| |
+        | | 0  0  0| | 3  1  3| | 0  1  3| |   | |2| |    | |2| |
+        |                                  |   |     |    |     |
+        | | 0  0  0| | 1  0  2| | 1  3  2| |   | |1| |    | |0| |
+        | | 0  0  0| | 0  1  4| | 2  1  3| |   | |0| |    | |2| |
+        | | 0  0  0| | 5  1  1| | 3  1  2| |   | |2| |    | |1| |
+    */
+
+    const int N = 3;
+    const int bz = 3;
+    const int nonZeroes = 5;
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+
+    using ft = typename Vector::field_type;
+    const ft myEps((ft)1e-6);
+    Dune::TestSuite t;
+
+    Matrix A(N, N, nonZeroes, Matrix::row_wise);
+    for (auto row = A.createbegin(); row != A.createend(); ++row) {
+        if (row.index() == 0) {
+            row.insert(row.index());
+        }
+        else if (row.index() == 1) {
+            row.insert(row.index());
+            row.insert(row.index() + 1);
+        }
+        else if (row.index() == 2) {
+            row.insert(row.index() - 1);
+            row.insert(row.index());
+        }
+    }
+
+    A[0][0][0][0]=3.0;    A[1][1][0][0]=1.0;    A[1][2][0][0]=1.0;
+    A[0][0][0][1]=1.0;    A[1][1][0][1]=0.0;    A[1][2][0][1]=0.0;
+    A[0][0][0][2]=2.0;    A[1][1][0][2]=1.0;    A[1][2][0][2]=2.0;
+    A[0][0][1][0]=2.0;    A[1][1][1][0]=4.0;    A[1][2][1][0]=0.0;
+    A[0][0][1][1]=3.0;    A[1][1][1][1]=1.0;    A[1][2][1][1]=1.0;
+    A[0][0][1][2]=1.0;    A[1][1][1][2]=0.0;    A[1][2][1][2]=1.0;
+    A[0][0][2][0]=2.0;    A[1][1][2][0]=3.0;    A[1][2][2][0]=0.0;
+    A[0][0][2][1]=1.0;    A[1][1][2][1]=1.0;    A[1][2][2][1]=1.0;
+    A[0][0][2][2]=0.0;    A[1][1][2][2]=3.0;    A[1][2][2][2]=3.0;
+
+    A[2][1][0][0]=1.0;    A[2][2][0][0]=1.0;
+    A[2][1][0][1]=0.0;    A[2][2][0][1]=3.0;
+    A[2][1][0][2]=2.0;    A[2][2][0][2]=2.0;
+    A[2][1][1][0]=0.0;    A[2][2][1][0]=2.0;
+    A[2][1][1][1]=1.0;    A[2][2][1][1]=1.0;
+    A[2][1][1][2]=4.0;    A[2][2][1][2]=3.0;
+    A[2][1][2][0]=5.0;    A[2][2][2][0]=3.0;
+    A[2][1][2][1]=1.0;    A[2][2][2][1]=1.0;
+    A[2][1][2][2]=1.0;    A[2][2][2][2]=2.0;
+
+    Vector x(3);
+    x[0][0] = 1.0;    x[1][0] = 1.0;    x[2][0] = 1.0;
+    x[0][1] = 2.0;    x[1][1] = 3.0;    x[2][1] = 0.0;
+    x[0][2] = 3.0;    x[1][2] = 2.0;    x[2][2] = 2.0;
+
+    Vector b(3);
+    b[0][0] = 2.0;    b[1][0] = 2.0;    b[2][0] = 0.0;
+    b[0][1] = 1.0;    b[1][1] = 3.0;    b[2][1] = 2.0;
+    b[0][2] = 2.0;    b[1][2] = 2.0;    b[2][2] = 1.0;
+
+
+    // D_00 = A_00
+    auto D_00 = A[0][0];
+    auto D_00_inv = D_00;
+    D_00_inv.invert();
+    // D_11 = A_11 - L_10 D_00_inv U_01
+    //      = A_11
+    auto D_11 = A[1][1];
+    auto D_11_inv = D_11;
+    D_11_inv.invert();
+    // D_22 = A_22 - A_20 D_00_inv A_02 - A_21 D_11_inv A_12
+    //      = A_22 - A_21 D_11_inv A_12
+    auto D_22 = A[2][2] - A[2][1]*D_11_inv*A[1][2];
+    auto D_22_inv = D_22;
+    D_22_inv.invert();
+
+    // define: z = M^-1(b - Ax)
+    // where: M = (D + L_A) A^-1 (D + U_A)
+    // lower triangular solve: (D + L) y = b - Ax
+
+    Vector y(3);
+    // y_0 = D_00_inv*[b_0 - (A_00*x_0 + A_01*x_1)]
+    //     = D_00_inv*[b_0 - A_00*x_0]
+    auto rhs = b[0];
+    A[0][0].mmv(x[0], rhs);
+    D_00_inv.mv(rhs, y[0]);
+
+    // y_1 = D_11_inv*(b_1 - (A_10*x_0 + A_11*x_1 + A_12*x_2) - A_10*y_0)
+    //     = D_11_inv*(b_1 - A_11*x_1)
+    rhs = b[1];
+    A[1][1].mmv(x[1], rhs);
+    A[1][2].mmv(x[2], rhs);
+    D_11_inv.mv(rhs, y[1]);
+
+    // y_2 = D_22_inv*(b_2 - (A_20*x_0 + A_21*x_1 + A_22*x_2) - (A_20*y_0 + A_21*y_1))
+    //     = D_22_inv*(b_2 - (A_21*x_1 + A_22*x_2) - (A_21*y_1))
+    rhs = b[2];
+    A[2][1].mmv(x[1], rhs);
+    A[2][2].mmv(x[2], rhs);
+    A[2][1].mmv(y[1], rhs);
+    D_22_inv.mv(rhs, y[2]);
+
+
+    // upper triangular solve: (D + U) z = Dy
+    Vector z(3);
+    // z_2 = y_2
+    z[2] = y[2];
+
+    // z_1 = y_1 - D_11_inv*A_12*z_2
+    z[1] = y[1];
+    auto temp = D_11_inv*A[1][2];
+    temp.mmv(z[2], z[1]);
+
+    // z_0 = y_0 - D_00_inv(A_01*z_1 + A_02*z_2)
+    // z_0 = y_0
+    z[0] = y[0];
+
+    // z is now the update x_k+1 - x_k
+
+    Dune::SeqDILU<Matrix,Vector,Vector> seqDILU(A, 1);
+    Vector v(3);
+    Vector d(b);
+    A.mmv(x, d);
+    seqDILU.apply(v, d);
+
+    for (int i = 0; i < 3; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            t.check(std::abs(v[i][j] - z[i][j]) <= myEps)
+             << " Error in SeqDILUApplyIsCorrect2, v[" << i <<"][" << j << "]=" << v[i][j] << " != z[" << i <<"][" << j << "]=" << z[i][j];
+        }
+    }
+}
+
+
+void seqDILUApplyIsEqualToSeqILUApply()
+{
+    /*
+        Tests that applying the DILU preconditioner is equivalent to applying an ILU preconditioner
+        for a block diagonal matrix.
+
+                 A               x     =     b
+        | | 3  1| | 0  0| |   | |1| |     | |2| |
+        | | 2  1| | 0  0| |   | |2| |     | |1| |
+        |                 | x |     |  =  |     |
+        | | 0  0| |-1  0| |   | |1| |     | |3| |
+        | | 0  0| | 0 -1| |   | |1| |     | |4| |
+    */
+
+    const int N = 2;
+    const int bz = 2;
+    const int nonZeroes = 2;
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, bz, bz>>;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double, bz>>;
+
+    using ft = typename Vector::field_type;
+    const ft myEps((ft)1e-6);
+    Dune::TestSuite t;
+
+    Matrix A(N, N, nonZeroes, Matrix::row_wise);
+    for (auto row = A.createbegin(); row != A.createend(); ++row) {
+        row.insert(row.index());
+    }
+
+    A[0][0][0][0]=3.0;
+    A[0][0][0][1]=1.0;
+    A[0][0][1][0]=2.0;
+    A[0][0][1][1]=1.0;
+
+    A[1][1][0][0]=-1.0;
+    A[1][1][1][1]=-1.0;
+
+    Vector x(2);
+    x[0][0] = 1.0;
+    x[0][1] = 2.0;
+    x[1][0] = 1.0;
+    x[1][1] = 1.0;
+
+    Vector b(2);
+    b[0][0] = 2.0;
+    b[0][1] = 1.0;
+    b[1][0] = 3.0;
+    b[1][1] = 4.0;
+
+    Dune::SeqDILU<Matrix,Vector,Vector> seqDILU(A, 1);
+    Vector v_DILU(2);
+    Vector d_DILU(b);
+    A.mmv(x, d_DILU);
+    seqDILU.apply(v_DILU, d_DILU);
+
+    Dune::SeqILU<Matrix,Vector,Vector> seqILU(A, 1);
+    Vector v_ILU(2);
+    Vector d_ILU(b);
+    A.mmv(x, d_ILU);
+    seqDILU.apply(v_ILU, d_ILU);
+
+    for (int i = 0; i < 2; ++i) {
+        for (int j = 0; j < 2; ++j) {
+            t.check(std::abs(v_DILU[i][j] - v_ILU[i][j]) <= myEps)
+            << " Error in SeqDILUApplyIsEqualToSeqILUApply, v_DILU[" << i <<"][" << j << "]=" << v_DILU[i][j] << " != v_ILU[" << i <<"][" << j << "]=" << v_ILU[i][j];
+        }
+    }
+}
+
+
+
+
+
+
+int main() try
+{
+  seqDILUApplyIsCorrect1();
+  seqDILUApplyIsCorrect2();
+  seqDILUApplyIsEqualToSeqILUApply();
+
+  return 0;
+}
+catch (std::exception& e) {
+  std::cerr << e.what() << std::endl;
+  return 1;
+}

dune/istl/test/preconditionerstest.cc

@@ -30,6 +30,7 @@ namespace Dune
   template class SeqSOR<Mat1, Vec1, Vec1>;
   template class SeqSSOR<Mat1, Vec1, Vec1>;
   template class Richardson<Vec1, Vec1>;
+  template class SeqDILU<Mat1, Vec1, Vec1>;
   template class SeqILU<Mat1, Vec1, Vec1>;
   template class SeqILDL<Mat1, Vec1, Vec1>;
 
@@ -37,6 +38,7 @@ namespace Dune
   template class SeqSOR<Mat2, Vec2, Vec2>;
   template class SeqSSOR<Mat2, Vec2, Vec2>;
   template class Richardson<Vec2, Vec2>;
+  template class SeqDILU<Mat2, Vec2, Vec2>;
   template class SeqILU<Mat2, Vec2, Vec2>;
   template class SeqILDL<Mat2, Vec2, Vec2>;
 
@@ -94,6 +96,10 @@ void testAllPreconditioners(const Matrix& matrix, const Vector& b)
   SeqJac<Matrix,Vector,Vector> seqJac(matrix, 1,1.0);
   testPreconditioner(matrix, b, x, seqJac);
 
+  x = 0;
+  SeqDILU<Matrix,Vector,Vector> seqDILU(matrix, 1.1);
+  testPreconditioner(matrix, b, x, seqDILU);
+
   x = 0;
   SeqILU<Matrix,Vector,Vector> seqILU(matrix, 3, 1.2, true);
   testPreconditioner(matrix, b, x, seqILU);

dune/python/istl/preconditioners.hh

@@ -164,6 +164,20 @@ namespace Dune
               The Jacobi iteration can only be applied if the matrix diagonal entries are all non-zeros.
         )doc" );
 
+      module.def( "SeqDILU", [] ( const M &A, field_type w ) {
+          return static_cast< Preconditioner * >( new SeqDILU< M, X, Y >( A, w ) );
+        }, "matrix"_a, "relaxation"_a = field_type( 1 ),
+        R"doc(Diagonal incomplete LU factorization preconditioner
+
+          Args:
+              matrix:      matrix to precondition
+              relaxation:  factor to relax the iterations (default: 1)
+
+          Returns:
+              ISTL Sequential diagonal incomplete LU factorization preconditioner
+
+        )doc" );
+
       module.def( "SeqILU", [] ( const M &A, int n, field_type w ) {
           return static_cast< Preconditioner * >( new SeqILU< M, X, Y >( A, n, w ) );
         }, "matrix"_a, "iterations"_a = 1, "relaxation"_a = field_type( 1 ),