]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Extand PreconditionRelaxation 13159/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 31 Dec 2021 15:30:21 +0000 (16:30 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Sat, 8 Jan 2022 15:07:30 +0000 (16:07 +0100)
doc/news/changes/minor/20220101Munch [new file with mode: 0644]
include/deal.II/lac/precondition.h
include/deal.II/matrix_free/operators.h
tests/lac/precondition_relaxation_01.cc [new file with mode: 0644]
tests/lac/precondition_relaxation_01.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20220101Munch b/doc/news/changes/minor/20220101Munch
new file mode 100644 (file)
index 0000000..a10abb7
--- /dev/null
@@ -0,0 +1,8 @@
+Improved: The class PreconditionRelaxation has been refactored.
+Similarly to PreconditionChebyshev, it can be now used
+stand-alone and users can attach their own preconditioners.
+Furthermore, users can specify multiple iteration steps, 
+which is particularly useful if PreconditionRelaxation is 
+used as smoother in a multigrid algorithm.
+<br>
+(Peter Munch, 2022/01/01)
index 4b0377e8a381f6c60ca12963e25a318d9934aaf6..cc2a1ac3fb693242ab2a54c95ee6c04573b7d14b 100644 (file)
@@ -29,6 +29,7 @@
 #include <deal.II/base/utilities.h>
 
 #include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/block_vector_base.h>
 #include <deal.II/lac/diagonal_matrix.h>
 #include <deal.II/lac/identity_matrix.h>
 #include <deal.II/lac/solver_cg.h>
@@ -418,13 +419,20 @@ public:
     /**
      * Constructor.
      */
-    AdditionalData(const double relaxation = 1.);
+    AdditionalData(const double       relaxation   = 1.,
+                   const unsigned int n_iterations = 1);
 
     /**
      * Relaxation parameter.
      */
     double relaxation;
 
+    /**
+     * Number of smoothing steps to be performed.
+     */
+    unsigned int n_iterations;
+
+
     /*
      * Preconditioner.
      */
@@ -500,6 +508,11 @@ protected:
    */
   double relaxation;
 
+  /**
+   * Number of smoothing steps to be performed.
+   */
+  unsigned int n_iterations;
+
   /*
    * Preconditioner.
    */
@@ -514,6 +527,114 @@ namespace internal
 {
   namespace PreconditionRelaxation
   {
+    template <typename T, typename VectorType>
+    struct has_step
+    {
+    private:
+      static bool
+      detect(...);
+
+      template <typename U>
+      static decltype(
+        std::declval<U const>().step(std::declval<VectorType &>(),
+                                     std::declval<const VectorType &>()))
+      detect(const U &);
+
+    public:
+      static const bool value =
+        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+    };
+
+    template <typename T, typename VectorType>
+    const bool has_step<T, VectorType>::value;
+
+    template <typename T, typename VectorType>
+    struct has_Tstep
+    {
+    private:
+      static bool
+      detect(...);
+
+      template <typename U>
+      static decltype(
+        std::declval<U const>().Tstep(std::declval<VectorType &>(),
+                                      std::declval<const VectorType &>()))
+      detect(const U &);
+
+    public:
+      static const bool value =
+        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+    };
+
+    template <typename T, typename VectorType>
+    const bool has_Tstep<T, VectorType>::value;
+
+    template <typename T, typename VectorType>
+    struct has_jacobi_step
+    {
+    private:
+      static bool
+      detect(...);
+
+      template <typename U>
+      static decltype(
+        std::declval<U const>().Jacobi_step(std::declval<VectorType &>(),
+                                            std::declval<const VectorType &>(),
+                                            1.0))
+      detect(const U &);
+
+    public:
+      static const bool value =
+        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+    };
+
+    template <typename T, typename VectorType>
+    const bool has_jacobi_step<T, VectorType>::value;
+
+    template <typename T, typename VectorType>
+    struct has_SOR_step
+    {
+    private:
+      static bool
+      detect(...);
+
+      template <typename U>
+      static decltype(
+        std::declval<U const>().SOR_step(std::declval<VectorType &>(),
+                                         std::declval<const VectorType &>(),
+                                         1.0))
+      detect(const U &);
+
+    public:
+      static const bool value =
+        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+    };
+
+    template <typename T, typename VectorType>
+    const bool has_SOR_step<T, VectorType>::value;
+
+    template <typename T, typename VectorType>
+    struct has_SSOR_step
+    {
+    private:
+      static bool
+      detect(...);
+
+      template <typename U>
+      static decltype(
+        std::declval<U const>().SSOR_step(std::declval<VectorType &>(),
+                                          std::declval<const VectorType &>(),
+                                          1.0))
+      detect(const U &);
+
+    public:
+      static const bool value =
+        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+    };
+
+    template <typename T, typename VectorType>
+    const bool has_SSOR_step<T, VectorType>::value;
+
     template <typename MatrixType>
     class PreconditionJacobiImpl
     {
@@ -538,13 +659,28 @@ namespace internal
         this->vmult(dst, src);
       }
 
-      template <typename VectorType>
+      template <
+        typename VectorType,
+        typename std::enable_if<has_jacobi_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
       void
       step(VectorType &dst, const VectorType &src) const
       {
         this->A->Jacobi_step(dst, src, this->relaxation);
       }
 
+      template <
+        typename VectorType,
+        typename std::enable_if<!has_jacobi_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
+      void
+      step(VectorType &, const VectorType &) const
+      {
+        Assert(false,
+               ExcMessage(
+                 "Matrix A does not provide a Jacobi_step() function!"));
+      }
+
       template <typename VectorType>
       void
       Tstep(VectorType &dst, const VectorType &src) const
@@ -581,20 +717,48 @@ namespace internal
         this->A->precondition_TSOR(dst, src, this->relaxation);
       }
 
-      template <typename VectorType>
+      template <
+        typename VectorType,
+        typename std::enable_if<has_SOR_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
       void
       step(VectorType &dst, const VectorType &src) const
       {
         this->A->SOR_step(dst, src, this->relaxation);
       }
 
-      template <typename VectorType>
+      template <
+        typename VectorType,
+        typename std::enable_if<!has_SOR_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
+      void
+      step(VectorType &, const VectorType &) const
+      {
+        Assert(false,
+               ExcMessage("Matrix A does not provide a SOR_step() function!"));
+      }
+
+      template <
+        typename VectorType,
+        typename std::enable_if<has_SOR_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
       void
       Tstep(VectorType &dst, const VectorType &src) const
       {
         this->A->TSOR_step(dst, src, this->relaxation);
       }
 
+      template <
+        typename VectorType,
+        typename std::enable_if<!has_SOR_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
+      void
+      Tstep(VectorType &, const VectorType &) const
+      {
+        Assert(false,
+               ExcMessage("Matrix A does not provide a TSOR_step() function!"));
+      }
+
     private:
       const SmartPointer<const MatrixType> A;
       const double                         relaxation;
@@ -658,13 +822,27 @@ namespace internal
                                    pos_right_of_diagonal);
       }
 
-      template <typename VectorType>
+      template <
+        typename VectorType,
+        typename std::enable_if<has_SSOR_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
       void
       step(VectorType &dst, const VectorType &src) const
       {
         this->A->SSOR_step(dst, src, this->relaxation);
       }
 
+      template <
+        typename VectorType,
+        typename std::enable_if<!has_SSOR_step<MatrixType, VectorType>::value,
+                                MatrixType>::type * = nullptr>
+      void
+      step(VectorType &, const VectorType &) const
+      {
+        Assert(false,
+               ExcMessage("Matrix A does not provide a SSOR_step() function!"));
+      }
+
       template <typename VectorType>
       void
       Tstep(VectorType &dst, const VectorType &src) const
@@ -724,104 +902,175 @@ namespace internal
       const std::vector<size_type> &inverse_permutation;
     };
 
-    template <typename T, typename VectorType>
-    struct has_step
-    {
-    private:
-      static bool
-      detect(...);
-
-      template <typename U>
-      static decltype(
-        std::declval<U const>().step(std::declval<VectorType &>(),
-                                     std::declval<const VectorType &>()))
-      detect(const U &);
-
-    public:
-      static const bool value =
-        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
-    };
-
-    template <typename T, typename VectorType>
-    const bool has_step<T, VectorType>::value;
-
-    template <typename T, typename VectorType>
-    struct has_Tstep
-    {
-    private:
-      static bool
-      detect(...);
-
-      template <typename U>
-      static decltype(
-        std::declval<U const>().Tstep(std::declval<VectorType &>(),
-                                      std::declval<const VectorType &>()))
-      detect(const U &);
-
-    public:
-      static const bool value =
-        !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
-    };
-
-    template <typename T, typename VectorType>
-    const bool has_Tstep<T, VectorType>::value;
-
     template <
+      typename MatrixType,
       typename PreconditionerType,
       typename VectorType,
       typename std::enable_if<has_step<PreconditionerType, VectorType>::value,
                               PreconditionerType>::type * = nullptr>
     void
-    step(const PreconditionerType &preconditioner,
+    step(const MatrixType &,
+         const PreconditionerType &preconditioner,
          VectorType &              dst,
-         const VectorType &        src)
+         const VectorType &        src,
+         const double              relaxation,
+         VectorType &,
+         VectorType &)
     {
+      Assert(relaxation == 1.0, ExcInternalError());
+
+      (void)relaxation;
+
       preconditioner.step(dst, src);
     }
 
     template <
+      typename MatrixType,
       typename PreconditionerType,
       typename VectorType,
       typename std::enable_if<!has_step<PreconditionerType, VectorType>::value,
                               PreconditionerType>::type * = nullptr>
     void
-    step(const PreconditionerType &preconditioner,
+    step(const MatrixType &        A,
+         const PreconditionerType &preconditioner,
          VectorType &              dst,
-         const VectorType &        src)
+         const VectorType &        src,
+         const double              relaxation,
+         VectorType &              residual,
+         VectorType &              tmp)
     {
-      Assert(false, ExcNotImplemented());
-      (void)preconditioner;
-      (void)dst;
-      (void)src;
+      residual.reinit(dst, true);
+      tmp.reinit(dst, true);
+
+      A.vmult(residual, dst);
+      residual.sadd(-1.0, 1.0, src);
+
+      preconditioner.vmult(tmp, residual);
+      dst.add(relaxation, tmp);
     }
 
     template <
+      typename MatrixType,
       typename PreconditionerType,
       typename VectorType,
       typename std::enable_if<has_Tstep<PreconditionerType, VectorType>::value,
                               PreconditionerType>::type * = nullptr>
     void
-    Tstep(const PreconditionerType &preconditioner,
+    Tstep(const MatrixType &,
+          const PreconditionerType &preconditioner,
           VectorType &              dst,
-          const VectorType &        src)
+          const VectorType &        src,
+          const double              relaxation,
+          VectorType &,
+          VectorType &)
     {
+      Assert(relaxation == 1.0, ExcInternalError());
+
+      (void)relaxation;
+
       preconditioner.Tstep(dst, src);
     }
 
     template <
+      typename MatrixType,
       typename PreconditionerType,
       typename VectorType,
       typename std::enable_if<!has_Tstep<PreconditionerType, VectorType>::value,
                               PreconditionerType>::type * = nullptr>
     void
-    Tstep(const PreconditionerType &preconditioner,
+    Tstep(const MatrixType &        A,
+          const PreconditionerType &preconditioner,
           VectorType &              dst,
-          const VectorType &        src)
+          const VectorType &        src,
+          const double              relaxation,
+          VectorType &              residual,
+          VectorType &              tmp)
     {
-      Assert(false, ExcNotImplemented());
-      (void)preconditioner;
-      (void)dst;
-      (void)src;
+      residual.reinit(dst, true);
+      tmp.reinit(dst, true);
+
+      A.Tvmult(residual, dst);
+      residual.sadd(-1.0, 1.0, src);
+
+      preconditioner.Tvmult(tmp, residual);
+      dst.add(relaxation, tmp);
+    }
+
+    template <typename MatrixType,
+              typename PreconditionerType,
+              typename VectorType>
+    void
+    step_operations(const MatrixType &        A,
+                    const PreconditionerType &preconditioner,
+                    VectorType &              dst,
+                    const VectorType &        src,
+                    const double              relaxation,
+                    VectorType &              tmp1,
+                    VectorType &              tmp2,
+                    const unsigned int        i,
+                    const bool                transposed)
+    {
+      if (i == 0)
+        {
+          if (transposed)
+            preconditioner.Tvmult(dst, src);
+          else
+            preconditioner.vmult(dst, src);
+
+          if (relaxation != 1.0)
+            dst *= relaxation;
+        }
+      else
+        {
+          if (transposed)
+            Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
+          else
+            step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
+        }
+    }
+
+    template <typename MatrixType,
+              typename VectorType,
+              typename std::enable_if<!IsBlockVector<VectorType>::value,
+                                      VectorType>::type * = nullptr>
+    void
+    step_operations(const MatrixType &                A,
+                    const DiagonalMatrix<VectorType> &preconditioner,
+                    VectorType &                      dst,
+                    const VectorType &                src,
+                    const double                      relaxation,
+                    VectorType &                      residual,
+                    VectorType &,
+                    const unsigned int i,
+                    const bool         transposed)
+    {
+      if (i == 0)
+        {
+          const auto dst_ptr  = dst.begin();
+          const auto src_ptr  = src.begin();
+          const auto diag_ptr = preconditioner.get_vector().begin();
+
+          for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
+            dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
+        }
+      else
+        {
+          residual.reinit(src, true);
+
+          const auto dst_ptr      = dst.begin();
+          const auto src_ptr      = src.begin();
+          const auto residual_ptr = residual.begin();
+          const auto diag_ptr     = preconditioner.get_vector().begin();
+
+          if (transposed)
+            A.Tvmult(residual, dst);
+          else
+            A.vmult(residual, dst);
+
+          for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
+            dst_ptr[i] +=
+              relaxation * (src_ptr[i] - residual_ptr[i]) * diag_ptr[i];
+        }
     }
 
   } // namespace PreconditionRelaxation
@@ -1752,9 +2001,13 @@ PreconditionRelaxation<MatrixType, PreconditionerType>::initialize(
   const MatrixType &    rA,
   const AdditionalData &parameters)
 {
-  A              = &rA;
-  relaxation     = parameters.relaxation;
+  A          = &rA;
+  relaxation = parameters.relaxation;
+
+  Assert(parameters.preconditioner, ExcNotInitialized());
+
   preconditioner = parameters.preconditioner;
+  n_iterations   = parameters.n_iterations;
 }
 
 
@@ -1762,7 +2015,8 @@ template <typename MatrixType, typename PreconditionerType>
 inline void
 PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
 {
-  A = nullptr;
+  A              = nullptr;
+  preconditioner = nullptr;
 }
 
 template <typename MatrixType, typename PreconditionerType>
@@ -1790,10 +2044,14 @@ PreconditionRelaxation<MatrixType, PreconditionerType>::vmult(
   VectorType &      dst,
   const VectorType &src) const
 {
-  preconditioner->vmult(dst, src);
+  Assert(this->A != nullptr, ExcNotInitialized());
+  Assert(this->preconditioner != nullptr, ExcNotInitialized());
 
-  if (this->relaxation != 1.0)
-    dst *= this->relaxation;
+  VectorType tmp1, tmp2;
+
+  for (unsigned int i = 0; i < n_iterations; ++i)
+    internal::PreconditionRelaxation::step_operations(
+      *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
 }
 
 template <typename MatrixType, typename PreconditionerType>
@@ -1803,10 +2061,14 @@ PreconditionRelaxation<MatrixType, PreconditionerType>::Tvmult(
   VectorType &      dst,
   const VectorType &src) const
 {
-  preconditioner->Tvmult(dst, src);
+  Assert(this->A != nullptr, ExcNotInitialized());
+  Assert(this->preconditioner != nullptr, ExcNotInitialized());
+
+  VectorType tmp1, tmp2;
 
-  if (this->relaxation != 1.0)
-    dst *= this->relaxation;
+  for (unsigned int i = 0; i < n_iterations; ++i)
+    internal::PreconditionRelaxation::step_operations(
+      *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
 }
 
 template <typename MatrixType, typename PreconditionerType>
@@ -1816,7 +2078,14 @@ PreconditionRelaxation<MatrixType, PreconditionerType>::step(
   VectorType &      dst,
   const VectorType &src) const
 {
-  internal::PreconditionRelaxation::step(*preconditioner, dst, src);
+  Assert(this->A != nullptr, ExcNotInitialized());
+  Assert(this->preconditioner != nullptr, ExcNotInitialized());
+
+  VectorType tmp1, tmp2;
+
+  for (unsigned int i = 1; i <= n_iterations; ++i)
+    internal::PreconditionRelaxation::step_operations(
+      *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
 }
 
 template <typename MatrixType, typename PreconditionerType>
@@ -1826,7 +2095,14 @@ PreconditionRelaxation<MatrixType, PreconditionerType>::Tstep(
   VectorType &      dst,
   const VectorType &src) const
 {
-  internal::PreconditionRelaxation::Tstep(*preconditioner, dst, src);
+  Assert(this->A != nullptr, ExcNotInitialized());
+  Assert(this->preconditioner != nullptr, ExcNotInitialized());
+
+  VectorType tmp1, tmp2;
+
+  for (unsigned int i = 1; i <= n_iterations; ++i)
+    internal::PreconditionRelaxation::step_operations(
+      *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
 }
 
 //---------------------------------------------------------------------------
@@ -1839,7 +2115,8 @@ PreconditionJacobi<MatrixType>::initialize(const MatrixType &    A,
   Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
 
   AdditionalData parameters;
-  parameters.relaxation = 1.0;
+  parameters.relaxation   = 1.0;
+  parameters.n_iterations = parameters_in.n_iterations;
   parameters.preconditioner =
     std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
 
@@ -1856,7 +2133,8 @@ PreconditionSOR<MatrixType>::initialize(const MatrixType &    A,
   Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
 
   AdditionalData parameters;
-  parameters.relaxation = 1.0;
+  parameters.relaxation   = 1.0;
+  parameters.n_iterations = parameters_in.n_iterations;
   parameters.preconditioner =
     std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
 
@@ -1873,7 +2151,8 @@ PreconditionSSOR<MatrixType>::initialize(const MatrixType &    A,
   Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
 
   AdditionalData parameters;
-  parameters.relaxation = 1.0;
+  parameters.relaxation   = 1.0;
+  parameters.n_iterations = parameters_in.n_iterations;
   parameters.preconditioner =
     std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
 
@@ -1895,7 +2174,8 @@ PreconditionPSOR<MatrixType>::initialize(
   Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
 
   typename BaseClass::AdditionalData parameters;
-  parameters.relaxation = 1.0;
+  parameters.relaxation   = 1.0;
+  parameters.n_iterations = parameters_in.n_iterations;
   parameters.preconditioner =
     std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
 
@@ -1952,8 +2232,9 @@ PreconditionUseMatrix<MatrixType, VectorType>::vmult(
 
 template <typename MatrixType, typename PreconditionerType>
 inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
-  AdditionalData(const double relaxation)
+  AdditionalData(const double relaxation, const unsigned int n_iterations)
   : relaxation(relaxation)
+  , n_iterations(n_iterations)
 {}
 
 
index 64f18fae96a94da7a64a4ec2bdf1ac541949c154..dd618bcd9a40200f1fafe09ffc4979ad8283bc60 100644 (file)
@@ -385,7 +385,6 @@ namespace MatrixFreeOperators
     const std::shared_ptr<DiagonalMatrix<VectorType>> &
     get_matrix_diagonal() const;
 
-
     /**
      * Apply the Jacobi preconditioner, which multiplies every element of the
      * <tt>src</tt> vector by the inverse of the respective diagonal element and
diff --git a/tests/lac/precondition_relaxation_01.cc b/tests/lac/precondition_relaxation_01.cc
new file mode 100644 (file)
index 0000000..cc0a890
--- /dev/null
@@ -0,0 +1,221 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test PreconditionRelaxation for different Jacobi preconditioners.
+
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+
+#include "../tests.h"
+
+template <typename VectorType>
+class MyDiagonalMatrix
+{
+public:
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    diagonal_matrix.vmult(dst, src);
+  }
+
+  void
+  Tvmult(VectorType &dst, const VectorType &src) const
+  {
+    diagonal_matrix.vmult(dst, src);
+  }
+
+  VectorType &
+  get_vector()
+  {
+    return diagonal_matrix.get_vector();
+  }
+
+private:
+  DiagonalMatrix<VectorType> diagonal_matrix;
+};
+
+
+
+template <typename PreconditionerType, typename VectorType>
+std::tuple<double, double, double, double>
+test(const PreconditionerType &preconditioner, const VectorType &src)
+{
+  VectorType dst;
+  dst.reinit(src);
+
+  preconditioner.vmult(dst, src);
+  const double norm_0 = dst.l2_norm();
+
+  preconditioner.step(dst, src);
+  const double norm_1 = dst.l2_norm();
+
+  preconditioner.Tvmult(dst, src);
+  const double norm_2 = dst.l2_norm();
+
+  preconditioner.Tstep(dst, src);
+  const double norm_3 = dst.l2_norm();
+
+  return std::tuple<double, double, double, double>{norm_0,
+                                                    norm_1,
+                                                    norm_2,
+                                                    norm_3};
+}
+
+
+
+int
+main()
+{
+  initlog();
+
+  using Number     = double;
+  using VectorType = Vector<Number>;
+  using MatrixType = SparseMatrix<Number>;
+
+  const unsigned int dim    = 2;
+  const unsigned int degree = 1;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(3);
+
+  FE_Q<dim>      fe(degree);
+  QGauss<dim>    quad(degree + 1);
+  MappingQ1<dim> mapping;
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MatrixType system_matrix;
+
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp);
+
+  SparsityPattern sparsity_pattern;
+  sparsity_pattern.copy_from(dsp);
+  system_matrix.reinit(sparsity_pattern);
+
+  MatrixCreator::create_laplace_matrix(mapping,
+                                       dof_handler,
+                                       quad,
+                                       system_matrix);
+
+  VectorType diagonal(dof_handler.n_dofs());
+  VectorType src(dof_handler.n_dofs());
+
+  for (const auto &entry : system_matrix)
+    if (entry.row() == entry.column())
+      diagonal[entry.row()] = 1.0 / entry.value();
+
+  src = 1.0;
+
+  std::vector<double>       relaxations{1.0, 0.9, 1.1};
+  std::vector<unsigned int> n_iterationss{1, 2, 3};
+
+  for (const auto relaxation : relaxations)
+    for (const auto n_iterations : n_iterationss)
+      {
+        std::vector<std::tuple<double, double, double, double>> results;
+
+        {
+          // Test PreconditionJacobi
+          PreconditionJacobi<MatrixType> preconditioner;
+
+          PreconditionJacobi<MatrixType>::AdditionalData ad;
+          ad.relaxation   = relaxation;
+          ad.n_iterations = n_iterations;
+
+          preconditioner.initialize(system_matrix, ad);
+
+          results.emplace_back(test(preconditioner, src));
+        }
+
+        {
+          // Test PreconditionRelaxation + DiagonalMatrix: optimized path
+          using PreconditionerType = DiagonalMatrix<VectorType>;
+
+          PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
+
+          PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData
+            ad;
+          ad.relaxation     = relaxation;
+          ad.n_iterations   = n_iterations;
+          ad.preconditioner = std::make_shared<PreconditionerType>();
+          ad.preconditioner->get_vector() = diagonal;
+
+          preconditioner.initialize(system_matrix, ad);
+
+          results.emplace_back(test(preconditioner, src));
+        }
+
+        {
+          // Test PreconditionRelaxation + wrapper around DiagonalMatrix:
+          // optimized path cannot be used
+          using PreconditionerType = MyDiagonalMatrix<VectorType>;
+
+          PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
+
+          PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData
+            ad;
+          ad.relaxation     = relaxation;
+          ad.n_iterations   = n_iterations;
+          ad.preconditioner = std::make_shared<PreconditionerType>();
+          ad.preconditioner->get_vector() = diagonal;
+
+          preconditioner.initialize(system_matrix, ad);
+
+          results.emplace_back(test(preconditioner, src));
+        }
+
+        if (std::equal(results.begin(),
+                       results.end(),
+                       results.begin(),
+                       [](const auto &a, const auto &b) {
+                         if (std::abs(std::get<0>(a) - std::get<0>(b)) > 1e-6)
+                           return false;
+
+                         if (std::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6)
+                           return false;
+
+                         if (std::abs(std::get<2>(a) - std::get<2>(b)) > 1e-6)
+                           return false;
+
+                         if (std::abs(std::get<3>(a) - std::get<3>(b)) > 1e-6)
+                           return false;
+
+                         return true;
+                       }))
+          deallog << "OK!" << std::endl;
+        else
+          deallog << "ERROR!" << std::endl;
+      }
+}
diff --git a/tests/lac/precondition_relaxation_01.output b/tests/lac/precondition_relaxation_01.output
new file mode 100644 (file)
index 0000000..f738e07
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!
+DEAL::OK!

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.