]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add Tpetra wrappers for Trilinos Ifpack2 preconditioners 17183/head
authorJan Philipp Thiele <thiele@wias-berlin.de>
Thu, 18 Jul 2024 16:42:37 +0000 (18:42 +0200)
committerJan Philipp Thiele <thiele@wias-berlin.de>
Fri, 19 Jul 2024 11:42:49 +0000 (13:42 +0200)
17 files changed:
cmake/configure/configure_20_trilinos.cmake
doc/doxygen/references.bib
doc/news/changes/major/20240718Thiele [new file with mode: 0644]
include/deal.II/base/config.h.in
include/deal.II/lac/trilinos_tpetra_precondition.h [new file with mode: 0644]
include/deal.II/lac/trilinos_tpetra_precondition.templates.h [new file with mode: 0644]
include/deal.II/lac/trilinos_tpetra_types.h
source/lac/CMakeLists.txt
source/lac/trilinos_tpetra_precondition.cc [new file with mode: 0644]
tests/trilinos_tpetra/add_matrices_06.cc [new file with mode: 0644]
tests/trilinos_tpetra/add_matrices_06.output [new file with mode: 0644]
tests/trilinos_tpetra/direct_solver_01.cc
tests/trilinos_tpetra/direct_solver_01.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output
tests/trilinos_tpetra/direct_solver_02.cc
tests/trilinos_tpetra/direct_solver_02.with_p4est=true.with_trilinos_with_amesos2=true.mpirun=1.output
tests/trilinos_tpetra/precondition.cc [new file with mode: 0644]
tests/trilinos_tpetra/precondition.output [new file with mode: 0644]

index c80291be2453072b700c937c18276759548fa9a1..0c26dfa1469867985c09ce9c1947a5c7e8a5018c 100644 (file)
@@ -22,7 +22,7 @@ set(FEATURE_TRILINOS_DEPENDS MPI)
 # A list of optional Trilinos modules we use:
 #
 set(_deal_ii_trilinos_optional_modules
-  Amesos2 Belos EpetraExt Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan 
+  Amesos2 Belos EpetraExt Ifpack2 Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan 
 )
 
 #
index 8892274e69597bcca47b08345005998d190c4358..ebb2b052724b3a404841774f0569bdf6b2350f5a 100644 (file)
    url     = {http://sepwww.stanford.edu/oldsep/joe/Reprints/8IWSA.pdf},
 }
 
+@article{BFKY2011,
+  title={Multigrid smoothers for ultraparallel computing},
+  author={Baker, Allison H and Falgout, Robert D and Kolev, Tzanio V and Yang, Ulrike Meier},
+  journal={SIAM Journal on Scientific Computing},
+  volume={33},
+  number={5},
+  pages={2864--2887},
+  year={2011},
+  publisher={SIAM}
+}
diff --git a/doc/news/changes/major/20240718Thiele b/doc/news/changes/major/20240718Thiele
new file mode 100644 (file)
index 0000000..670190c
--- /dev/null
@@ -0,0 +1,30 @@
+New: LinearAlgebra::TpetraWrappers preconditioner interfaces
+for Trilinos Ifpack2, as listed below:
+<ul>
+    <li>Preconditioners that were also available in the old interface.
+        <ul>
+            <li>PreconditionIdentity</li>
+            <li>PreconditionJacobi</li>
+            <li>PreconditionSOR</li>
+            <li>PreconditionSSOR</li>
+            <li>PreconditionChebyshev</li>
+            <li>PreconditionILU</li>
+            <li>PreconditionILUT</li>
+            <li>PreconditionBlockJacobi</li>
+            <li>PreconditionBlockSOR</li>
+            <li>PreconditionBlockSSOR</li>
+        </ul>
+    </li>
+    <li>Newly introduced preconditioners
+        <ul>
+            <li>PreconditionL1Jacobi<br>
+                - A new variant to improve coupling in MPI-parallel applications</li>
+            <li>PreconditionL1GaussSeidel<br>
+                - A new variant to improve coupling in MPI-parallel applications</li>
+            <li>PreconditionIfpack<br>
+                - A generic Ifpack2 preconditioner for any
+                classes or parameters not interfaced (yet)</li>
+        </ul>
+    </li>
+</ul>
+(Jan Philipp Thiele, 2024/07/18)
\ No newline at end of file
index 5546e052c85448d0d90a4eb05071427f229e0ccc..e09209db4d1e6bd82b8dd2dd4fc028f8a7cf1f1d 100644 (file)
 #cmakedefine DEAL_II_TRILINOS_WITH_AMESOS2
 #cmakedefine DEAL_II_TRILINOS_WITH_BELOS
 #cmakedefine DEAL_II_TRILINOS_WITH_EPETRAEXT
+#cmakedefine DEAL_II_TRILINOS_WITH_IFPACK2
 #cmakedefine DEAL_II_TRILINOS_WITH_MUELU
 #cmakedefine DEAL_II_TRILINOS_WITH_NOX
 #cmakedefine DEAL_II_TRILINOS_WITH_ROL
diff --git a/include/deal.II/lac/trilinos_tpetra_precondition.h b/include/deal.II/lac/trilinos_tpetra_precondition.h
new file mode 100644 (file)
index 0000000..9650d73
--- /dev/null
@@ -0,0 +1,1329 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+#ifndef dealii_trilinos_tpetra_precondition_h
+#define dealii_trilinos_tpetra_precondition_h
+
+
+#include <deal.II/base/config.h>
+
+#include "deal.II/base/exceptions.h"
+#include "deal.II/base/memory_space.h"
+
+#include "deal.II/lac/vector.h"
+
+#include <Teuchos_BLAS_types.hpp>
+
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+
+#  include <deal.II/base/subscriptor.h>
+
+#  include <deal.II/lac/la_parallel_vector.h>
+#  include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#  include <Teuchos_ConfigDefs.hpp>
+#  include <Teuchos_ParameterList.hpp>
+#  include <Teuchos_RCPDecl.hpp>
+#  include <Tpetra_Operator.hpp>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+
+    /**
+     * The base class for all preconditioners based on Tpetra sparse matrices.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionBase : public Subscriptor
+    {
+    public:
+      /**
+       * @brief Standardized data struct to pipe additional flags to the
+       * preconditioner.
+       *
+       */
+      struct AdditionalData
+      {};
+
+      /**
+       * @brief Constructor. Does not do anything. The <tt>initialize</tt> function of
+       * the derived classes will have to create the preconditioner from a given
+       * sparse matrix.
+       *
+       */
+      PreconditionBase() = default;
+
+      /**
+       * @brief Desctructor.
+       * Destroys the preconditioner, leaving an object like just after having
+       * called the constructor.
+       */
+      void
+      clear();
+
+      /**
+       * @brief Apply the preconditioner.
+       *
+       * @param dst Input vector to apply the preconditioner to
+       * @param src Result vector
+       */
+      virtual void
+      vmult(Vector<Number, MemorySpace>       &dst,
+            const Vector<Number, MemorySpace> &src) const;
+
+      /**
+       * @brief Apply the transpose preconditioner
+       *
+       * @param dst Input vector to apply the preconditioner to
+       * @param src Result vector
+       */
+      virtual void
+      Tvmult(Vector<Number, MemorySpace>       &dst,
+             const Vector<Number, MemorySpace> &src) const;
+
+      /**
+       * @brief Apply the preconditioner
+       *
+       * @param dst Input vector to apply the preconditioner to
+       * @param src Result vector
+       */
+      virtual void
+      vmult(dealii::Vector<Number> &dst, dealii::Vector<Number> &src) const;
+
+      /**
+       * @brief Apply the transpose preconditioner
+       */
+      virtual void
+      Tvmult(dealii::Vector<Number> &dst, dealii::Vector<Number> &src) const;
+
+      /**
+       * @brief Access to underlying Trilinos data
+       *
+       * Calling this function from an uninitialized object will cause an
+       * exception.
+       */
+      const TpetraTypes::LinearOperator<Number, MemorySpace> &
+      trilinos_operator() const;
+
+      /**
+       * @brief Access to underlying Trilinos data
+       *
+       * Calling this function from an uninitialized object will cause an
+       * exception.
+       */
+      Teuchos::RCP<TpetraTypes::LinearOperator<Number, MemorySpace>>
+      trilinos_rcp() const;
+
+      /**
+       * @name Partitioners
+       */
+      /** @{ */
+
+      /**
+       * Return the partitioning of the domain space of this matrix, i.e., the
+       * partitioning of the vectors this matrix has to be multiplied with.
+       *
+       * @return IndexSet of the domain of this preconditioner
+       */
+      IndexSet
+      locally_owned_domain_indices() const;
+
+      /**
+       * Return the partitioning of the range space of this matrix, i.e., the
+       * partitioning of the vectors that are result from matrix-vector
+       * products.
+       *
+       * @return IndexSet of the range of this preconditioner
+       */
+      IndexSet
+      locally_owned_range_indices() const;
+
+      /** @} */
+
+      /**
+       * @addtogroup Exceptions
+       */
+      /** @{ */
+      /**
+       * The maps of the underlying matrix and one of the vectors
+       * do not match.
+       */
+      DeclException1(
+        ExcNonMatchingMaps,
+        std::string,
+        << "The sparse matrix the preconditioner is based on "
+        << "uses a map that is not compatible to the one in vector " << arg1
+        << ". Check preconditioner and matrix setup.");
+
+      /**
+       * The chosen preconditioner does not support a transposed apply.
+       */
+      DeclExceptionMsg(
+        ExcTransposeNotSupported,
+        "The chosen preconditioner does not support transposing the matrix.");
+      /** @} */
+
+    protected:
+      /**
+       * This is a pointer to the preconditioner object that is used when
+       * applying the preconditioner.
+       */
+      Teuchos::RCP<TpetraTypes::LinearOperator<Number, MemorySpace>>
+        preconditioner;
+
+      /**
+       * @brief The list of preconditioner parameters.
+       *
+       * This structure is Trilinos counterpart to the AdditionalData structures
+       * in deal.II. Therefore any initialize will at some point pass this to
+       * the preconditioner. Most derived classes will handle building this list
+       * based on an AdditionalData object that exposes and defaults the most
+       * sensible parameters. But some classes will also offer full
+       * customization for experienced Trilinos users.
+       *
+       */
+      Teuchos::ParameterList parameter_list;
+    };
+
+
+
+#  ifdef DEAL_II_TRILINOS_WITH_IFPACK2
+    /**
+     * @brief Wrapper class for the IdentitySolver preconditioner of Ifpack2.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     *
+     * @tparam Number scalar type of the preconditioner
+     * @tparam MemorySpace
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionIdentity : public PreconditionBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief Construct identity preconditioner.
+       *
+       */
+      PreconditionIdentity() = default;
+
+      /**
+       * @brief Initializes the preconditioner for the matrix <tt>A</tt>.
+       *
+       * Note that this only needs the matrix for information on the parallel
+       * distribution of the matrix and will return the original vector on
+       * <tt>vmult</tt> or <tt>Tvmult</tt>.
+       *
+       * @param A Matrix to base the preconditioner on.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A);
+    };
+
+
+
+    /**
+     * @brief The base class for all Ifpack2 preconditioners which are handled through its Factory.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionIfpackBase : public PreconditionBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief Constructor.
+       *
+       * This constructor will only save the type name given.
+       * The actual preconditioner has to be constructed through the
+       * initialize method of one of the underlying classes.
+       *
+       * @param preconditioner_type
+       */
+      PreconditionIfpackBase(const std::string &preconditioner_type);
+
+      /**
+       * Initializes the preconditioner for the matrix <tt>A</tt> based on
+       * the <tt>parameter_set</tt>.
+       *
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A);
+
+      /**
+       * The chosen preconditioner is not supported or configured with Ifpack2.
+       */
+      DeclException1(
+        ExcTrilinosIpack2PreconditionerUnsupported,
+        std::string,
+        << "You tried to select the preconditioner type <" << arg1 << ">\n"
+        << "but this preconditioner is not supported by Trilinos/Ifpack22\n"
+        << "due to one of the following reasons:\n"
+        << "* This preconditioner does not exist\n"
+        << "* This preconditioner has a specialized constructor not supported by the Ifpack2 Factory.\n"
+        << "* This preconditioner is not (yet) supported by Trilinos/Ifpack2\n"
+        << "* Trilinos/Ifpack2 was not configured for its use.");
+
+    protected:
+      /**
+       * The set preconditioner type to be handed to
+       * <tt>Ifpack2::Factory::create()</tt>
+       */
+      std::string preconditioner_type;
+    };
+
+
+
+    /**
+     * @brief Wrapper to create custom Ifpack2 preconditioners.
+     *
+     * This general purpose class supports any preconditioner of Ifpack2
+     * that can be constructed through its <tt>Ifpack2::Factory::create()</tt>
+     * function.
+     *
+     * @note This is a class for users familiar with Trilinos.
+     *
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionIfpack
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief Construct a new custom Ifpack2 preconditioner.
+       *
+       * The currently tested options for the type are
+       * <ul>
+       * <li> "RELAXATION" </li>
+       * <li> "CHEBYSHEV" </li>
+       * <li> "RILUK" </li>
+       * <li> "ILUT" </li>
+       * <li> "SCHWARZ" </li>
+       * </ul>
+       *
+       * But please refer to the User's Guide of Ifpack2 for more information.
+       *
+       * @param preconditioner_type the type based on Ifpack2 notation
+       */
+      PreconditionIfpack(const std::string &preconditioner_type);
+
+      /**
+       * @brief Set the parameter list for the preconditioner.
+       *
+       * This list will be passed to the Ifpack2 preconditioner during
+       * initialization. For parameter options based on the chosen
+       * preconditioner type see the User's Guide of Ifpack2.
+       *
+       * @param parameter_list
+       */
+      void
+      set_parameter_list(Teuchos::ParameterList &parameter_list);
+    };
+
+
+
+    /**
+     * @brief The classical Jacobi preconditioner within Ifpack2.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     *
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionJacobi
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         * The defaults yield a textbook Jacobi preconditioner.
+         *
+         * @param omega The damping factor for the relaxation.
+         * @param fix_diagonal Fix small diagonal entries for the inversion?
+         * @param min_diagonal Threshold for fix_diagonal.
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult.
+         */
+        AdditionalData(const double omega        = 1.,
+                       const bool   fix_diagonal = false,
+                       const double min_diagonal = 0.,
+                       const int    n_sweeps     = 1);
+
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Whether or not to enlarge entries below a threshold
+         *
+         * If the matrix diagonal contains entries close to zero, their
+         * inversion will be close to division by zero. This can be corrected by
+         * setting this value to <tt>true</tt>. This will result in additional
+         * computational work during initialization.
+         */
+        bool fix_diagonal;
+        /**
+         * @brief Threshold below which entries will be fixed.
+         *
+         * If the threshold is zero (default) only entries which are exactly
+         * zero will be replaced will small nonzero values.
+         *
+         */
+        double min_diagonal;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionJacobi();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * @brief l1 variant of the Jacobi preconditioner.
+     *
+     * This variant adds the l1-Norm of all off-processor entries of a local row
+     * i to its diagonal entry d_ii to improve coupling in MPI-parallel
+     * applications, as described and introduced in @cite BFKY2011.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionL1Jacobi
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         * The defaults yield the l1 Jacobi preconditioner specified in the
+         * original publication.
+         *
+         * @param omega The damping factor for the relaxation.
+         * @param fix_diagonal Fix small diagonal entries for the inversion?
+         * @param min_diagonal Threshold for fix_diagonal.
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult.
+         */
+        AdditionalData(const double omega        = 1.,
+                       const bool   fix_diagonal = false,
+                       const double min_diagonal = 0.,
+                       const int    n_sweeps     = 1);
+
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Whether or not to enlarge entries below a threshold
+         *
+         * If the matrix diagonal contains entries close to zero, their
+         * inversion will be close to division by zero. This can be corrected by
+         * setting this value to <tt>true</tt>. This will result in additional
+         * computational work during initialization.
+         */
+        bool fix_diagonal;
+        /**
+         * @brief Threshold below which entries will be fixed.
+         *
+         * If the threshold is zero (default) only entries which are exactly
+         * zero will be replaced will small nonzero values.
+         *
+         */
+        double min_diagonal;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionL1Jacobi();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * @brief l1 variant of the Gauss-Seidel preconditioner.
+     *
+     * This variant adds the l1-Norm of all off-processor entries of a single
+     * row i of the upper triangular part of A to its diagonal entry d_ii to
+     * improve coupling in MPI-parallel applications, , as described and
+     * introduced in @cite BFKY2011.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionL1GaussSeidel
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         * The defaults yield the l1 Gauss-Seidel preconditioner specified in
+         * the original publication.
+         *
+         * @param omega The damping factor for the relaxation.
+         * @param eta Threshold parameter for diagonal correction.
+         * @param fix_diagonal Fix small diagonal entries for the inversion?
+         * @param min_diagonal Threshold for fix_diagonal.
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult.
+         */
+        AdditionalData(const double omega        = 1,
+                       const double eta          = 1.5,
+                       const bool   fix_diagonal = false,
+                       const double min_diagonal = 0,
+                       const int    n_sweeps     = 1);
+
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+
+        /**
+         * @brief Threshold parameter for diagonal correction.
+         *
+         * The l1 correction for a row i will only be added if it is more than
+         * eta times larger than the original diagonal entry.
+         *
+         */
+        double eta;
+        /**
+         * @brief Whether or not to enlarge entries below a threshold
+         *
+         * If the matrix diagonal contains entries close to zero, their
+         * inversion will be close to division by zero. This can be corrected by
+         * setting this value to <tt>true</tt>. This will result in additional
+         * computational work during initialization.
+         */
+        bool fix_diagonal;
+        /**
+         * @brief Threshold below which entries will be fixed.
+         *
+         * If the threshold is zero (default) only entries which are exactly
+         * zero will be replaced will small nonzero values.
+         *
+         */
+        double min_diagonal;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionL1GaussSeidel();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * @brief The class for the SOR preconditioner within Ifpack2.
+     *
+     * If the code is executed MPI-parallel the individual processors will be
+     * connected with an additive Schwarz method.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionSOR : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         * The defaults yield a textbook Jacobi preconditioner.
+         *
+         * @param omega The damping factor for the relaxation.
+         * @param overlap Overlap between processor local matrices.
+         * @param fix_diagonal Fix small diagonal entries for the inversion?
+         * @param min_diagonal Threshold for fix_diagonal.
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult.
+         */
+        AdditionalData(const double omega        = 1.,
+                       const int    overlap      = 0,
+                       const bool   fix_diagonal = false,
+                       const double min_diagonal = 0.,
+                       const int    n_sweeps     = 1);
+
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Overlap between processor local matrices.
+         *
+         * The amount of overlap between the processor local matrix blocks
+         * used in the underlying additive Schwarz method.
+         *
+         * The default will yield a block Jacobi preconditioner with each
+         * processor forming its own local block.
+         */
+        int overlap;
+        /**
+         * @brief Whether or not to enlarge entries below a threshold
+         *
+         * If the matrix diagonal contains entries close to zero, their
+         * inversion will be close to division by zero. This can be corrected by
+         * setting this value to <tt>true</tt>. This will result in additional
+         * computational work during initialization.
+         */
+        bool fix_diagonal;
+        /**
+         * @brief Threshold below which entries will be fixed.
+         *
+         * If the threshold is zero (default) only entries which are exactly
+         * zero will be replaced will small nonzero values.
+         *
+         */
+        double min_diagonal;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionSOR();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * The class for the classical SSOR preconditioner within Ifpack2.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionSSOR : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         * The defaults yield a textbook Jacobi preconditioner.
+         *
+         * @param omega The damping factor for the relaxation.
+         * @param overlap Overlap between processor local matrices.
+         * @param fix_diagonal Fix small diagonal entries for the inversion?
+         * @param min_diagonal Threshold for fix_diagonal.
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult.
+         */
+        AdditionalData(const double omega        = 1.,
+                       const int    overlap      = 0,
+                       const bool   fix_diagonal = false,
+                       const double min_diagonal = 0.,
+                       const int    n_sweeps     = 1);
+
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Overlap between processor local matrices.
+         *
+         * The amount of overlap between the processor local matrix blocks
+         * used in the underlying additive Schwarz method.
+         *
+         * The default will yield a block Jacobi preconditioner with each
+         * processor forming its own local block.
+         */
+        int overlap;
+        /**
+         * @brief Whether or not to enlarge entries below a threshold
+         *
+         * If the matrix diagonal contains entries close to zero, their
+         * inversion will be close to division by zero. This can be corrected by
+         * setting this value to <tt>true</tt>. This will result in additional
+         * computational work during initialization.
+         */
+        bool fix_diagonal;
+        /**
+         * @brief Threshold below which entries will be fixed.
+         *
+         * If the threshold is zero (default) only entries which are exactly
+         * zero will be replaced will small nonzero values.
+         *
+         */
+        double min_diagonal;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionSSOR();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * The class for the Chebyshev preconditioner within Ifpack2.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionChebyshev
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * @param degree Degree of the Chebyshev polynomial.
+         * @param max_eigenvalue Upper bound for the maximum eigenvalue of the matrix.
+         * @param min_eigenvalue Lower bound for the minimum eigenvalue of the matrix.
+         * @param eigenvalue_ratio Ratio between maximum and minimum eigenvalue.
+         * @param min_diagonal Threshold for increasing diagonal entries.
+         * @param nonzero_starting Do not zero starting entries of solution.
+         */
+        AdditionalData(const int    degree           = 1,
+                       const double max_eigenvalue   = 10.,
+                       const double min_eigenvalue   = 1.,
+                       const double eigenvalue_ratio = 30.,
+                       const double min_diagonal     = 1e-12,
+                       const bool   nonzero_starting = false);
+        /**
+         * @brief Degree of the Chebyshev polynomial.
+         *
+         * The degree directly corresponds to the number of matrix-vector
+         * products that have to be performed during a single application of the
+         * vmult() and Tvmult() operations.
+         *
+         */
+        int degree;
+        /*
+         * @brief Upper bound for the maximum eigenvalue of the matrix.
+         *
+         * This needs to be set properly for the appropriate performance of this
+         * preconditioner.
+         */
+        double max_eigenvalue;
+        /*
+         * @brief Lower bound for the minimum eigenvalue of the matrix.
+         *
+         */
+        double min_eigenvalue;
+        /**
+         * @brief Estimated ratio between maximum and mimimum eigenvalue.
+         *
+         */
+        double eigenvalue_ratio;
+        /**
+         * @brief Threshold below which entries will be fixed.
+         *
+         * If the threshold is zero (default) only entries which are exactly
+         * zero will be replaced will small nonzero values.
+         *
+         */
+        double min_diagonal;
+        /**
+         * @brief Do not zero starting entries of solution vector.
+         *
+         * The default (false) zeroes out the entries of dst during vmult() and
+         * Tvmult() which is the recommended setting.
+         *
+         * However, in some situations (e.g. high-frequency error smoothing) it
+         * can be useful to append previous data to the Chebyshev corrections.
+         * The user should really know what they are doing when setting this
+         * flag to true.
+         *
+         */
+        bool nonzero_starting;
+      };
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionChebyshev();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * @brief The ILU/ILU(K)/RILU(K) preconditioner.
+     *
+     * The class for the modified incomplete LU factorization preconditioner
+     * RILUK within Ifpack2. This preconditioner works both in serial and
+     * parallel, depending on the matrix it is based on.
+     *
+     * In general, an incomplete factorization only has values on the sparsity
+     * pattern of the matrix, but one can gradually increase <tt>ilu_fill</tt>
+     * to eventually obtain a full LU factorization.
+     *
+     * For parallel applications this preconditioner implements an additive
+     * Schwarz preconditioner with RILUK as local smoothing on each processor.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionILU : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Before factorization each diagonal entry will be modified by the
+         * following fomula \f[a_{ii}^{new} = \alpha\sign(a_{ii})+\beta a_{ii}
+         * \f] with \f[\alpha\f] given by ilu_atol and \f[\beta\f] given by
+         * ilu_rtol.
+         *
+         * @param ilu_fill Amount of additional fill-in.
+         * @param ilu_atol Constant to be to each diagonal entry before factorization.
+         * @param ilu_rtol Factor to scale all diagonal entries by before factorization.
+         * @param overlapOverlap between processor local matrices.
+         */
+        AdditionalData(const int    ilu_fill = 0,
+                       const double ilu_atol = 0.,
+                       const double ilu_rtol = 1.,
+                       const int    overlap  = 0);
+
+        /**
+         * @brief Amount of additional fill-in.
+         *
+         * Level-of-fill to increase the sparsity pattern of the preconditioner.
+         * A large enough value will result in a complete LU factorization.
+         *
+         */
+        int ilu_fill;
+        /**
+         * @brief Constant to be added to each diagonal entry before factorization.
+         *
+         */
+        double ilu_atol;
+        /**
+         * @brief Factor to scale all diagonal entries by before factorization.
+         *
+         */
+        double ilu_rtol;
+        /**
+         * @brief Overlap between processor local matrices.
+         *
+         * The amount of overlap between the processor local matrix blocks
+         * used in the underlying additive Schwarz method.
+         *
+         * The default will yield a block Jacobi preconditioner with each
+         * processor forming its own local block.
+         */
+        int overlap;
+      };
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionILU();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * The class for the thresholded incomplete LU (ILUT) preconditioner within
+     * Ifpack2.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionILUT : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Before factorization each diagonal entry will be modified by the
+         * following fomula \f[a_{ii}^{new} = \alpha\sign(a_{ii})+\beta a_{ii}
+         * \f] with \f[\alpha\f] given by ilut_atol and \f[\beta\f] given by
+         * ilut_rtol.
+         *
+         * @param ilut_drop Threshold for dropping entries.
+         * @param ilut_fill Amount of additional fill-in.
+         * @param ilut_atol Constant to be to each diagonal entry before factorization.
+         * @param ilut_rtol Factor to scale all diagonal entries by before factorization.
+         * @param overlapOverlap between processor local matrices.
+         */
+        AdditionalData(const double ilut_drop = 0.,
+                       const double ilut_fill = 0.,
+                       const double ilut_atol = 0.,
+                       const double ilut_rtol = 1.,
+                       const int    overlap   = 0);
+
+        /**
+         * @brief Threshold for dropping entries.
+         *
+         * Together with <tt>ilut_fill</tt> this controls the amount of fill-in
+         * and the actual values to be used or dropped.
+         */
+        double ilut_drop;
+        /**
+         * @brief Amount of additional fill-in.
+         *
+         * Level-of-fill to increase the sparsity pattern of the preconditioner.
+         * A large enough value will result in a complete LU factorization.
+         * Note, however, that this will drastically increase the memory
+         * requirement, especially for 3d simulations.
+         *
+         */
+        double ilut_fill;
+        /**
+         * @brief Constant to be added to each diagonal entry before factorization.
+         *
+         */
+        double ilut_atol;
+        /**
+         * @brief Factor to scale all diagonal entries by before factorization.
+         *
+         */
+        double ilut_rtol;
+        /**
+         * @brief Overlap between processor local matrices.
+         *
+         * The amount of overlap between the processor local matrix blocks
+         * used in the underlying additive Schwarz method.
+         *
+         * The default will yield a block Jacobi preconditioner with each
+         * processor forming its own local block.
+         */
+        int overlap;
+      };
+
+
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionILUT();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * The class for the BlockJacobi preconditioner within Ifpack2.
+     *
+     * @note This preconditioner always uses linear partitioning.
+     * There are other partitioners available that need additional user provided
+     * data, so we suggest using PreconditionIfpack with a custom
+     * Teuchos::ParameterList. Details on the parameters can be found in the
+     * User's Guide of Ifpack2.
+     *
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionBlockJacobi
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         * The defaults yield a textbook Jacobi preconditioner.
+         *
+         * @param n_local_parts The number of blocks per processor.
+         * @param omega The damping factor for the relaxation.
+         * @param block_overlap Amount of overlap between blocks
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult
+         */
+        AdditionalData(const int    n_local_parts = 1,
+                       const double omega         = 1.,
+                       const int    block_overlap = 0,
+                       const int    n_sweeps      = 1);
+
+        /**
+         * @brief Number of blocks per processor.
+         *
+         * The default of 1 results in a single block per processor
+         * which corresponds to Preconditioner from PreconditionJacobi
+         */
+        int n_local_parts;
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Amount of overlap between blocks
+         *
+         */
+        int block_overlap;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+
+      /**
+       * @brief Constructor.
+       *
+       */
+      PreconditionBlockJacobi();
+
+      /**
+       * @brief Compute the preconditioner based on the given matrix and parameters.
+       *
+       * @param A The matrix to base the preconditioner on.
+       * @param additional_data The set of parameters to tune the preconditioner.
+       */
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * @brief The class for the Block SOR preconditioner within Ifpack2.
+     *
+     * If the code is executed MPI-parallel the individual processors will be
+     * connected with an additive Schwarz method.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionBlockSOR
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         *
+         * @param n_local_parts The number of blocks per processor.
+         * @param omega The damping factor for the relaxation.
+         * @param overlap Overlap between processors
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult
+         */
+        AdditionalData(const int    n_local_parts = 1,
+                       const double omega         = 1,
+                       const int    overlap       = 0,
+                       const int    n_sweeps      = 1);
+
+        /**
+         * @brief Number of blocks per processor.
+         *
+         * The default of 1 results in a single block per processor
+         * which corresponds to Preconditioner from PreconditionJacobi
+         */
+        int n_local_parts;
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Overlap between processor local matrices.
+         *
+         * The amount of overlap between the processor local matrix blocks
+         * used in the underlying additive Schwarz method.
+         *
+         * The default will yield a block Jacobi preconditioner with each
+         * processor forming its own local block.
+         */
+        int overlap;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+      PreconditionBlockSOR();
+
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+
+
+
+    /**
+     * @brief The class for the Block SSOR preconditioner within Ifpack2.
+     *
+     * If the code is executed MPI-parallel the individual processors will be
+     * connected with an additive Schwarz method.
+     *
+     * @ingroup TpetraWrappers
+     * @ingroup Preconditioners
+     */
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
+    class PreconditionBlockSSOR
+      : public PreconditionIfpackBase<Number, MemorySpace>
+    {
+    public:
+      /**
+       * @brief The set of additional parameters to tune the preconditioner.
+       *
+       */
+      struct AdditionalData
+      {
+        /**
+         * @brief Constructor.
+         *
+         * Set the parameters to be used in the preconditioner.
+         *
+         * @param n_local_parts The number of blocks per processor.
+         * @param omega The damping factor for the relaxation.
+         * @param overlap Overlap between processors
+         * @param n_sweeps Number of relaxation sweeps per call to vmult or Tvmult
+         */
+        AdditionalData(const int    n_local_parts = 1,
+                       const double omega         = 1,
+                       const int    overlap       = 1,
+                       const int    n_sweeps      = 1);
+        /**
+         * @brief Number of blocks per processor.
+         *
+         * The default of 1 results in a single block per processor
+         * which corresponds to Preconditioner from PreconditionJacobi
+         */
+        int n_local_parts;
+        /**
+         * @brief Relaxation damping factor.
+         *
+         */
+        double omega;
+        /**
+         * @brief Overlap between processor local matrices.
+         *
+         * The amount of overlap between the processor local matrix blocks
+         * used in the underlying additive Schwarz method.
+         *
+         * The default will yield a block Jacobi preconditioner with each
+         * processor forming its own local block.
+         */
+        int overlap;
+        /**
+         * @brief Set how often the preconditioner should be applied during vmult() or Tvmult().
+         *
+         */
+        int n_sweeps;
+      };
+
+      PreconditionBlockSSOR();
+
+      void
+      initialize(const SparseMatrix<Number, MemorySpace> &A,
+                 const AdditionalData &additional_data = AdditionalData());
+    };
+#  endif // DEAL_II_TRILINOS_WITH_IFPACK2
+
+  } // namespace TpetraWrappers
+} // namespace LinearAlgebra
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_TRILINOS_WITH_TPETRA
+
+#endif
diff --git a/include/deal.II/lac/trilinos_tpetra_precondition.templates.h b/include/deal.II/lac/trilinos_tpetra_precondition.templates.h
new file mode 100644 (file)
index 0000000..554f623
--- /dev/null
@@ -0,0 +1,805 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+#ifndef dealii_trilinos_tpetra_precondition_templates_h
+#define dealii_trilinos_tpetra_precondition_templates_h
+
+
+#include <deal.II/base/config.h>
+
+#include "deal.II/base/exceptions.h"
+#include <deal.II/base/index_set.h>
+
+#include "deal.II/lac/trilinos_tpetra_types.h"
+
+#include <string>
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+#  ifdef DEAL_II_TRILINOS_WITH_IFPACK2
+#    include <deal.II/lac/trilinos_tpetra_precondition.h>
+#    include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#    include <Ifpack2_Factory.hpp>
+#    include <Ifpack2_IdentitySolver.hpp>
+#    include <Ifpack2_IdentitySolver_decl.hpp>
+#    include <Ifpack2_Preconditioner.hpp>
+#    include <Kokkos_Core_fwd.hpp>
+#    include <Kokkos_DualView.hpp>
+#    include <Teuchos_RCPDecl.hpp>
+#    include <Tpetra_CrsMatrix_decl.hpp>
+#    include <Tpetra_Vector_decl.hpp>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    /* ---------------------- PreconditionBase ------------------------ */
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionBase<Number, MemorySpace>::clear()
+    {
+      preconditioner.reset();
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    const TpetraTypes::LinearOperator<Number, MemorySpace> &
+    PreconditionBase<Number, MemorySpace>::trilinos_operator() const
+    {
+      return *preconditioner;
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    Teuchos::RCP<TpetraTypes::LinearOperator<Number, MemorySpace>>
+    PreconditionBase<Number, MemorySpace>::trilinos_rcp() const
+    {
+      return preconditioner;
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    IndexSet
+    PreconditionBase<Number, MemorySpace>::locally_owned_domain_indices() const
+    {
+      return IndexSet(preconditioner->getDomainMap());
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    IndexSet
+    PreconditionBase<Number, MemorySpace>::locally_owned_range_indices() const
+    {
+      return IndexSet(preconditioner->getRangeMap());
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    PreconditionBase<Number, MemorySpace>::vmult(
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &src) const
+    {
+      Assert(dst.trilinos_vector().getMap()->isSameAs(
+               *(preconditioner->getRangeMap())),
+             ExcNonMatchingMaps("dst"));
+      Assert(src.trilinos_rcp()->getMap()->isSameAs(
+               *(preconditioner->getDomainMap())),
+             ExcNonMatchingMaps("src"));
+
+      preconditioner->apply(src.trilinos_vector(), dst.trilinos_vector());
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    PreconditionBase<Number, MemorySpace>::Tvmult(
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &src) const
+    {
+      Assert(preconditioner->hasTransposeApply(), ExcTransposeNotSupported());
+
+      Assert(dst.trilinos_rcp()->getMap()->isSameAs(
+               *(preconditioner->getDomainMap())),
+             ExcNonMatchingMaps("dst"));
+
+      Assert(src.trilinos_rcp()->getMap()->isSameAs(
+               *(preconditioner->getRangeMap())),
+             ExcNonMatchingMaps("src"));
+
+      preconditioner->apply(src.trilinos_vector(),
+                            dst.trilinos_vector(),
+                            Teuchos::TRANS);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    PreconditionBase<Number, MemorySpace>::vmult(
+      dealii::Vector<Number> &dst,
+      dealii::Vector<Number> &src) const
+    {
+      AssertDimension(dst.size(),
+                      preconditioner->getRangeMap()->getGlobalNumElements());
+      AssertDimension(src.size(),
+                      preconditioner->getDomainMap()->getGlobalNumElements());
+
+      // Construct Tpetra Vectors as views to the actual dealii::Vector data
+      // Since all of Tpetra is based on Kokkos::DualView this is a bit of a
+      // process
+
+      // 1. Make host view with underlying dealii::Vector data.
+      //    (This is always a n x 1 matrix since multivector expects that.)
+      TpetraTypes::HostViewTypeUnmanaged<Number> host_src(src.data(),
+                                                          src.size(),
+                                                          1);
+      TpetraTypes::HostViewTypeUnmanaged<Number> host_dst(dst.data(),
+                                                          dst.size(),
+                                                          1);
+
+      // 2. Create device mirror (if Device == Host, no alloc )
+      auto dev_src = Kokkos::create_mirror_view_and_copy(
+        typename MemorySpace::kokkos_space::execution_space(), host_src);
+
+      auto dev_dst = Kokkos::create_mirror_view_and_copy(
+        typename MemorySpace::kokkos_space::execution_space(), host_dst);
+
+      // 3. Create dual views from the previous ones
+      TpetraTypes::DualViewTypeUnmanaged<Number, MemorySpace> dual_src(
+        dev_src, host_src);
+      TpetraTypes::DualViewTypeUnmanaged<Number, MemorySpace> dual_dst(
+        dev_dst, host_dst);
+
+      // 4. Create Tpetra Vector from
+      TpetraTypes::VectorType<Number, MemorySpace> tpetra_dst(
+        preconditioner->getRangeMap(), dual_dst);
+      TpetraTypes::VectorType<Number, MemorySpace> tpetra_src(
+        preconditioner->getDomainMap(), dual_src);
+
+      // 5. Finally, apply our preconditioner
+      preconditioner->apply(tpetra_src, tpetra_dst);
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    inline void
+    PreconditionBase<Number, MemorySpace>::Tvmult(
+      dealii::Vector<Number> &dst,
+      dealii::Vector<Number> &src) const
+    {
+      AssertDimension(dst.size(),
+                      preconditioner->getDomainMap()->getGlobalNumElements());
+      AssertDimension(src.size(),
+                      preconditioner->getRangeMap()->getGlobalNumElements());
+
+      Assert(preconditioner->hasTransposeApply(), ExcTransposeNotSupported());
+      // Construct Tpetra Vectors as views to the actual dealii::Vector data
+      // Since all of Tpetra is based on Kokkos::DualView this is a bit of a
+      // process
+
+      // 1. Make host view with underlying dealii::Vector data.
+      //    (This is always a n x 1 matrix since multivector expects that.)
+      TpetraTypes::HostViewTypeUnmanaged<Number> host_src(src.data(),
+                                                          src.size(),
+                                                          1);
+      TpetraTypes::HostViewTypeUnmanaged<Number> host_dst(dst.data(),
+                                                          dst.size(),
+                                                          1);
+
+      // 2. Create device mirror (if Device == Host, no alloc )
+      auto dev_src = Kokkos::create_mirror_view_and_copy(
+        typename MemorySpace::kokkos_space::execution_space(), host_src);
+
+      auto dev_dst = Kokkos::create_mirror_view_and_copy(
+        typename MemorySpace::kokkos_space::execution_space(), host_dst);
+
+      // 3. Create dual views from the previous ones
+      TpetraTypes::DualViewTypeUnmanaged<Number, MemorySpace> dual_src(
+        dev_src, host_src);
+      TpetraTypes::DualViewTypeUnmanaged<Number, MemorySpace> dual_dst(
+        dev_dst, host_dst);
+
+      // 4. Create Tpetra Vector from
+      TpetraTypes::VectorType<Number, MemorySpace> tpetra_dst(
+        preconditioner->getRangeMap(), dual_dst);
+      TpetraTypes::VectorType<Number, MemorySpace> tpetra_src(
+        preconditioner->getDomainMap(), dual_src);
+
+      // 5. Finally, apply our preconditioner
+      preconditioner->apply(tpetra_src, tpetra_dst, Teuchos::TRANS);
+    }
+
+
+
+    /* ---------------------- PreconditionIdentity ------------------------ */
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionIdentity<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A)
+    {
+      Teuchos::RCP<Ifpack2::IdentitySolver<
+        TpetraTypes::RowMatrixType<Number, MemorySpace>>>
+        ifpack2Preconditioner;
+      ifpack2Preconditioner =
+        rcp(new Ifpack2::IdentitySolver<
+            TpetraTypes::RowMatrixType<Number, MemorySpace>>(A.trilinos_rcp()));
+
+      ifpack2Preconditioner->initialize();
+      ifpack2Preconditioner->compute();
+      this->preconditioner = ifpack2Preconditioner;
+    }
+
+
+
+    /* ---------------------- PreconditionIfpackBase ------------------------ */
+    template <typename Number, typename MemorySpace>
+    PreconditionIfpackBase<Number, MemorySpace>::PreconditionIfpackBase(
+      const std::string &preconditioner_type)
+      : preconditioner_type(preconditioner_type)
+    {
+      Ifpack2::Factory factory;
+      bool             supported =
+        factory.isSupported<TpetraTypes::MatrixType<Number, MemorySpace>>(
+          preconditioner_type);
+      AssertThrow(supported,
+                  ExcTrilinosIpack2PreconditionerUnsupported(
+                    preconditioner_type));
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionIfpackBase<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A)
+    {
+      Ifpack2::Factory factory;
+      Teuchos::RCP<TpetraTypes::Ifpack2PreconType<Number, MemorySpace>>
+        ifpack2Preconditioner;
+      ifpack2Preconditioner =
+        factory.create(preconditioner_type, A.trilinos_rcp());
+      ifpack2Preconditioner->setParameters(this->parameter_list);
+      ifpack2Preconditioner->initialize();
+      ifpack2Preconditioner->compute();
+      this->preconditioner = ifpack2Preconditioner;
+    }
+
+
+
+    /* ---------------------- PreconditionIfpack ------------------------ */
+    template <typename Number, typename MemorySpace>
+    PreconditionIfpack<Number, MemorySpace>::PreconditionIfpack(
+      const std::string &preconditioner_type)
+      : PreconditionIfpackBase<Number, MemorySpace>(preconditioner_type)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionIfpack<Number, MemorySpace>::set_parameter_list(
+      Teuchos::ParameterList &parameter_list)
+    {
+      this->parameter_list.setParameters(parameter_list);
+    }
+
+
+
+    /* ---------------------- PreconditionJacobi ------------------------ */
+    template <typename Number, typename MemorySpace>
+    PreconditionJacobi<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const double omega,
+      const bool   fix_diagonal,
+      const double min_diagonal,
+      const int    n_sweeps)
+      : omega(omega)
+      , fix_diagonal(fix_diagonal)
+      , min_diagonal(min_diagonal)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionJacobi<Number, MemorySpace>::PreconditionJacobi()
+      : PreconditionIfpackBase<Number, MemorySpace>("RELAXATION")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionJacobi<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("relaxation: type", "Jacobi");
+      this->parameter_list.set("relaxation: damping factor", ad.omega);
+      this->parameter_list.set("relaxation: sweeps", ad.n_sweeps);
+      this->parameter_list.set("relaxation: fix tiny diagonal entries",
+                               ad.fix_diagonal);
+      this->parameter_list.set("relaxation: min diagonal value",
+                               ad.min_diagonal);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionL1Jacobi ------------------------ */
+    template <typename Number, typename MemorySpace>
+    PreconditionL1Jacobi<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const double omega,
+      const bool   fix_diagonal,
+      const double min_diagonal,
+      const int    n_sweeps)
+      : omega(omega)
+      , fix_diagonal(fix_diagonal)
+      , min_diagonal(min_diagonal)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionL1Jacobi<Number, MemorySpace>::PreconditionL1Jacobi()
+      : PreconditionIfpackBase<Number, MemorySpace>("RELAXATION")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionL1Jacobi<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("relaxation: type", "Jacobi");
+      this->parameter_list.set("relaxation: use l1", true);
+      this->parameter_list.set("relaxation: damping factor", ad.omega);
+      this->parameter_list.set("relaxation: sweeps", ad.n_sweeps);
+      this->parameter_list.set("relaxation: fix tiny diagonal entries",
+                               ad.fix_diagonal);
+      this->parameter_list.set("relaxation: min diagonal value",
+                               ad.min_diagonal);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionL1GaussSeidel -------------------- */
+    template <typename Number, typename MemorySpace>
+    PreconditionL1GaussSeidel<Number, MemorySpace>::AdditionalData::
+      AdditionalData(const double omega,
+                     const double eta,
+                     const bool   fix_diagonal,
+                     const double min_diagonal,
+                     const int    n_sweeps)
+      : omega(omega)
+      , eta(eta)
+      , fix_diagonal(fix_diagonal)
+      , min_diagonal(min_diagonal)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionL1GaussSeidel<Number, MemorySpace>::PreconditionL1GaussSeidel()
+      : PreconditionIfpackBase<Number, MemorySpace>("RELAXATION")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionL1GaussSeidel<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("relaxation: type", "Gauss-Seidel");
+      this->parameter_list.set("relaxation: use l1", true);
+      this->parameter_list.set("relaxation: damping factor", ad.omega);
+      this->parameter_list.set("relaxation: l1 eta", ad.eta);
+      this->parameter_list.set("relaxation: sweeps", ad.n_sweeps);
+      this->parameter_list.set("relaxation: fix tiny diagonal entries",
+                               ad.fix_diagonal);
+      this->parameter_list.set("relaxation: min diagonal value",
+                               ad.min_diagonal);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionSOR ------------------------ */
+    template <typename Number, typename MemorySpace>
+    PreconditionSOR<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const double omega,
+      const int    overlap,
+      const bool   fix_diagonal,
+      const double min_diagonal,
+      const int    n_sweeps)
+      : omega(omega)
+      , overlap(overlap)
+      , fix_diagonal(fix_diagonal)
+      , min_diagonal(min_diagonal)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionSOR<Number, MemorySpace>::PreconditionSOR()
+      : PreconditionIfpackBase<Number, MemorySpace>("SCHWARZ")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionSOR<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("schwarz: subdomain solver name", "RELAXATION");
+      {
+        Teuchos::ParameterList &relaxParams =
+          this->parameter_list.sublist("schwarz: subdomain solver parameters",
+                                       false);
+        relaxParams.set("relaxation: type", "Symmetric Gauss-Seidel");
+        relaxParams.set("relaxation: damping factor", ad.omega);
+        relaxParams.set("relaxation: sweeps", ad.n_sweeps);
+        relaxParams.set("relaxation: fix tiny diagonal entries",
+                        ad.fix_diagonal);
+        relaxParams.set("relaxation: min diagonal value", ad.min_diagonal);
+      }
+      this->parameter_list.set("schwarz: combine mode", "ADD");
+      this->parameter_list.set("schwarz: overlap level", ad.overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionSSOR ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionSSOR<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const double omega,
+      const int    overlap,
+      const bool   fix_diagonal,
+      const double min_diagonal,
+      const int    n_sweeps)
+      : omega(omega)
+      , overlap(overlap)
+      , fix_diagonal(fix_diagonal)
+      , min_diagonal(min_diagonal)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionSSOR<Number, MemorySpace>::PreconditionSSOR()
+      : PreconditionIfpackBase<Number, MemorySpace>("SCHWARZ")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionSSOR<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("schwarz: subdomain solver name", "RELAXATION");
+      {
+        Teuchos::ParameterList &relaxParams =
+          this->parameter_list.sublist("schwarz: subdomain solver parameters",
+                                       false);
+        relaxParams.set("relaxation: type", "Symmetric Gauss-Seidel");
+        relaxParams.set("relaxation: damping factor", ad.omega);
+        relaxParams.set("relaxation: sweeps", ad.n_sweeps);
+        relaxParams.set("relaxation: fix tiny diagonal entries",
+                        ad.fix_diagonal);
+        relaxParams.set("relaxation: min diagonal value", ad.min_diagonal);
+      }
+      this->parameter_list.set("schwarz: combine mode", "ADD");
+      this->parameter_list.set("schwarz: overlap level", ad.overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionBlockJacobi ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionBlockJacobi<Number, MemorySpace>::AdditionalData::
+      AdditionalData(const int    n_local_parts,
+                     const double omega,
+                     const int    block_overlap,
+                     const int    n_sweeps)
+      : n_local_parts(n_local_parts)
+      , omega(omega)
+      , block_overlap(block_overlap)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionBlockJacobi<Number, MemorySpace>::PreconditionBlockJacobi()
+      : PreconditionIfpackBase<Number, MemorySpace>("BLOCK_RELAXATION")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionBlockJacobi<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("relaxation: type", "Jacobi");
+      this->parameter_list.set("relaxation: damping factor", ad.omega);
+      this->parameter_list.set("relaxation: sweeps", ad.n_sweeps);
+      this->parameter_list.set("partitioner: local parts", ad.n_local_parts);
+      this->parameter_list.set("partitioner: overlap", ad.block_overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionBlockSOR ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionBlockSOR<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const int    n_local_parts,
+      const double omega,
+      const int    overlap,
+      const int    n_sweeps)
+      : n_local_parts(n_local_parts)
+      , omega(omega)
+      , overlap(overlap)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionBlockSOR<Number, MemorySpace>::PreconditionBlockSOR()
+      : PreconditionIfpackBase<Number, MemorySpace>("SCHWARZ")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionBlockSOR<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("schwarz: subdomain solver name",
+                               "BLOCK_RELAXATION");
+      {
+        Teuchos::ParameterList &relaxParams =
+          this->parameter_list.sublist("schwarz: subdomain solver parameters",
+                                       false);
+        relaxParams.set("relaxation: type", "Gauss-Seidel");
+        relaxParams.set("relaxation: damping factor", ad.omega);
+        relaxParams.set("relaxation: sweeps", ad.n_sweeps);
+        relaxParams.set("partitioner: local parts", ad.n_local_parts);
+      }
+      this->parameter_list.set("schwarz: combine mode", "ADD");
+      this->parameter_list.set("schwarz: overlap level", ad.overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionBlockSSOR ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionBlockSSOR<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const int    n_local_parts,
+      const double omega,
+      const int    overlap,
+      const int    n_sweeps)
+      : n_local_parts(n_local_parts)
+      , omega(omega)
+      , overlap(overlap)
+      , n_sweeps(n_sweeps)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionBlockSSOR<Number, MemorySpace>::PreconditionBlockSSOR()
+      : PreconditionIfpackBase<Number, MemorySpace>("SCHWARZ")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionBlockSSOR<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("schwarz: subdomain solver name",
+                               "BLOCK_RELAXATION");
+      {
+        Teuchos::ParameterList &relaxParams =
+          this->parameter_list.sublist("schwarz: subdomain solver parameters",
+                                       false);
+        relaxParams.set("relaxation: type", "Symmetric Gauss-Seidel");
+        relaxParams.set("relaxation: damping factor", ad.omega);
+        relaxParams.set("relaxation: sweeps", ad.n_sweeps);
+        relaxParams.set("partitioner: local parts", ad.n_local_parts);
+      }
+      this->parameter_list.set("schwarz: combine mode", "ADD");
+      this->parameter_list.set("schwarz: overlap level", ad.overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionChebyshev ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionChebyshev<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const int    degree,
+      const double max_eigenvalue,
+      const double eigenvalue_ratio,
+      const double min_eigenvalue,
+      const double min_diagonal,
+      const bool   nonzero_starting)
+      : degree(degree)
+      , max_eigenvalue(max_eigenvalue)
+      , min_eigenvalue(min_eigenvalue)
+      , eigenvalue_ratio(eigenvalue_ratio)
+      , min_diagonal(min_diagonal)
+      , nonzero_starting(nonzero_starting)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionChebyshev<Number, MemorySpace>::PreconditionChebyshev()
+      : PreconditionIfpackBase<Number, MemorySpace>("CHEBYSHEV")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionChebyshev<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("chebyshev: degree", ad.degree);
+      this->parameter_list.set("chebyshev: max eigenvalue", ad.max_eigenvalue);
+      this->parameter_list.set("chebyshev: min eigenvalue", ad.min_eigenvalue);
+      this->parameter_list.set("chebyshev: ratio eigenvalue",
+                               ad.eigenvalue_ratio);
+      this->parameter_list.set("chebyshev: min diagonal value",
+                               ad.min_diagonal);
+      this->parameter_list.set("chebyshev: zero starting solution",
+                               !ad.nonzero_starting);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionILU ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionILU<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const int    ilu_fill,
+      const double ilu_atol,
+      const double ilu_rtol,
+      const int    overlap)
+      : ilu_fill(ilu_fill)
+      , ilu_atol(ilu_atol)
+      , ilu_rtol(ilu_rtol)
+      , overlap(overlap)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionILU<Number, MemorySpace>::PreconditionILU()
+      : PreconditionIfpackBase<Number, MemorySpace>("SCHWARZ")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionILU<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("schwarz: subdomain solver name", "RILUK");
+      {
+        Teuchos::ParameterList &rilukParams =
+          this->parameter_list.sublist("schwarz: subdomain solver parameters",
+                                       false);
+        rilukParams.set("fact: iluk level-of-fill", ad.ilu_fill);
+        rilukParams.set("fact: absolute threshold", ad.ilu_atol);
+        rilukParams.set("fact: relative threshold", ad.ilu_rtol);
+      }
+      this->parameter_list.set("schwarz: combine mode", "ADD");
+      this->parameter_list.set("schwarz: overlap level", ad.overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+
+
+
+    /* ---------------------- PreconditionILUT ----------------------*/
+    template <typename Number, typename MemorySpace>
+    PreconditionILUT<Number, MemorySpace>::AdditionalData::AdditionalData(
+      const double ilut_drop,
+      const double ilut_fill,
+      const double ilut_atol,
+      const double ilut_rtol,
+      const int    overlap)
+      : ilut_drop(ilut_drop)
+      , ilut_fill(ilut_fill)
+      , ilut_atol(ilut_atol)
+      , ilut_rtol(ilut_rtol)
+      , overlap(overlap)
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    PreconditionILUT<Number, MemorySpace>::PreconditionILUT()
+      : PreconditionIfpackBase<Number, MemorySpace>("SCHWARZ")
+    {}
+
+
+
+    template <typename Number, typename MemorySpace>
+    void
+    PreconditionILUT<Number, MemorySpace>::initialize(
+      const SparseMatrix<Number, MemorySpace> &A,
+      const AdditionalData                    &ad)
+    {
+      this->parameter_list.set("schwarz: subdomain solver name", "ILUT");
+      {
+        Teuchos::ParameterList &ilutParams =
+          this->parameter_list.sublist("schwarz: subdomain solver parameters",
+                                       false);
+        ilutParams.set("fact: drop tolerance", ad.ilut_drop);
+        ilutParams.set("fact: ilut level-of-fill", ad.ilut_fill);
+        ilutParams.set("fact: absolute threshold", ad.ilut_atol);
+        ilutParams.set("fact: relative threshold", ad.ilut_rtol);
+      }
+      this->parameter_list.set("schwarz: combine mode", "ADD");
+      this->parameter_list.set("schwarz: overlap level", ad.overlap);
+      PreconditionIfpackBase<Number, MemorySpace>::initialize(A);
+    }
+  } // namespace TpetraWrappers
+} // namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif // DEAL_II_TRILINOS_WITH_IFPACK2
+#endif   // DEAL_II_TRILINOS_WITH_TPETRA
+
+#endif
index 5dda027e86d57b8a859530a23ebe0a0e0b00e248..bf38dc778442b61d36ab4f4319cadcbb3cb7c217 100644 (file)
 
 #ifdef DEAL_II_TRILINOS_WITH_TPETRA
 
+#  include <Kokkos_DualView.hpp>
 #  include <Tpetra_Details_DefaultTypes.hpp>
 
+
 // Forward declarations
 #  ifndef DOXYGEN
+#    include <Ifpack2_Preconditioner.hpp>
 #    include <Tpetra_CrsGraph_fwd.hpp>
 #    include <Tpetra_CrsMatrix_fwd.hpp>
 #    include <Tpetra_Export_fwd.hpp>
@@ -32,8 +35,8 @@
 #    include <Tpetra_Map_fwd.hpp>
 #    include <Tpetra_MultiVector_fwd.hpp>
 #    include <Tpetra_Operator_fwd.hpp>
+#    include <Tpetra_RowMatrix_fwd.hpp>
 #    include <Tpetra_Vector_fwd.hpp>
-
 #  endif
 
 DEAL_II_NAMESPACE_OPEN
@@ -49,6 +52,8 @@ namespace LinearAlgebra
        * local ordinate (processor local indices).
        */
       using LO = int;
+
+
       /**
        * global ordinate (global indices).
        */
@@ -72,7 +77,6 @@ namespace LinearAlgebra
         typename MemorySpace::kokkos_space>;
 #  endif
 
-
       /**
        * Communication between processors.
        */
@@ -128,6 +132,55 @@ namespace LinearAlgebra
       template <typename Number, typename MemorySpace>
       using MatrixType =
         Tpetra::CrsMatrix<Number, LO, GO, NodeType<MemorySpace>>;
+
+
+      /**
+       * Tpetra type for a parallel distributed row matrix.
+       */
+      template <typename Number, typename MemorySpace>
+      using RowMatrixType =
+        Tpetra::RowMatrix<Number, LO, GO, NodeType<MemorySpace>>;
+
+
+      /**
+       * @brief Typedef for the Kokkos::View type of unmanaged host memory.
+       *
+       * Unmanaged means that the view is not allocating data itself, but only
+       * uses an external pointer. Consequently, the memory is not freed upon
+       * destruction of the view.
+       * This is needed for shallow copies of deal.II LA structures
+       * to Trilinos LA structures.
+       *
+       */
+      template <typename Number>
+      using HostViewTypeUnmanaged = Kokkos::View<Number **,
+                                                 Kokkos::LayoutLeft,
+                                                 Kokkos::HostSpace,
+                                                 Kokkos::MemoryUnmanaged>;
+
+      /**
+       * @brief Typedef for the Kokkos::View type of unmanaged memory.
+       *
+       * Unmanaged means that the view is not allocating data itself, but only
+       * uses an external pointer. Consequently, the memory is not freed upon
+       * destruction of the view.
+       * This is needed for shallow copies of deal.II LA structures
+       * to Trilinos LA structures.
+       */
+      template <typename Number, typename MemorySpace>
+      using DualViewTypeUnmanaged =
+        Kokkos::DualView<Number **,
+                         Kokkos::LayoutLeft,
+                         typename MemorySpace::kokkos_space::execution_space,
+                         Kokkos::MemoryUnmanaged>;
+
+      /**
+       * Type for a Trilinos preconditioner from the Ifpack2 package.
+       */
+      template <typename Number, typename MemorySpace>
+      using Ifpack2PreconType =
+        Ifpack2::Preconditioner<Number, LO, GO, NodeType<MemorySpace>>;
+
     } // namespace TpetraTypes
   }   // namespace TpetraWrappers
 } // namespace LinearAlgebra
index 5dec8854307932a09d88240230447a3ae92182e5..f77a42be6a45b0a03db0a867fb788acd2685be71 100644 (file)
@@ -145,6 +145,7 @@ if(DEAL_II_WITH_TRILINOS)
     trilinos_tpetra_block_vector.cc
     trilinos_tpetra_block_sparse_matrix.cc
     trilinos_tpetra_communication_pattern.cc
+    trilinos_tpetra_precondition.cc
     trilinos_tpetra_solver_direct.cc
     trilinos_tpetra_sparse_matrix.cc
     trilinos_tpetra_sparsity_pattern.cc
diff --git a/source/lac/trilinos_tpetra_precondition.cc b/source/lac/trilinos_tpetra_precondition.cc
new file mode 100644 (file)
index 0000000..21f7ee2
--- /dev/null
@@ -0,0 +1,109 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+#include "deal.II/lac/trilinos_tpetra_precondition.templates.h"
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+#  ifdef DEAL_II_TRILINOS_WITH_IFPACK2
+
+DEAL_II_NAMESPACE_OPEN
+
+#    ifndef DOXYGEN
+// explicit instantiations
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+#      ifdef HAVE_TPETRA_INST_FLOAT
+    template class PreconditionBase<float>;
+    template class PreconditionIdentity<float>;
+    template class PreconditionIfpackBase<float>;
+    template class PreconditionIfpack<float>;
+    template class PreconditionJacobi<float>;
+    template class PreconditionL1Jacobi<float>;
+    template class PreconditionL1GaussSeidel<float>;
+    template class PreconditionSOR<float>;
+    template class PreconditionSSOR<float>;
+    template class PreconditionBlockJacobi<float>;
+    template class PreconditionBlockSOR<float>;
+    template class PreconditionBlockSSOR<float>;
+    template class PreconditionChebyshev<float>;
+    template class PreconditionILU<float>;
+    template class PreconditionILUT<float>;
+#      endif
+
+#      ifdef HAVE_TPETRA_INST_DOUBLE
+    template class PreconditionBase<double>;
+    template class PreconditionIdentity<double>;
+    template class PreconditionIfpackBase<double>;
+    template class PreconditionIfpack<double>;
+    template class PreconditionJacobi<double>;
+    template class PreconditionL1Jacobi<double>;
+    template class PreconditionL1GaussSeidel<double>;
+    template class PreconditionSOR<double>;
+    template class PreconditionSSOR<double>;
+    template class PreconditionBlockJacobi<double>;
+    template class PreconditionBlockSOR<double>;
+    template class PreconditionBlockSSOR<double>;
+    template class PreconditionChebyshev<double>;
+    template class PreconditionILU<double>;
+    template class PreconditionILUT<double>;
+#      endif
+#      ifdef DEAL_II_WITH_COMPLEX_VALUES
+#        ifdef HAVE_TPETRA_INST_COMPLEX_FLOAT
+    template class PreconditionBase<std::complex<float>>;
+    template class PreconditionIdentity<std::complex<float>>;
+    template class PreconditionIfpackBase<std::complex<float>>;
+    template class PreconditionIfpack<std::complex<float>>;
+    template class PreconditionJacobi<std::complex<float>>;
+    template class PreconditionL1Jacobi<std::complex<float>>;
+    template class PreconditionL1GaussSeidel<std::complex<float>>;
+    template class PreconditionSOR<std::complex<float>>;
+    template class PreconditionSSOR<std::complex<float>>;
+    template class PreconditionBlockJacobi<std::complex<float>>;
+    template class PreconditionBlockSOR<std::complex<float>>;
+    template class PreconditionBlockSSOR<std::complex<float>>;
+    template class PreconditionChebyshev<std::complex<float>>;
+    template class PreconditionILU<std::complex<float>>;
+    template class PreconditionILUT<std::complex<float>>;
+#        endif
+#        ifdef HAVE_TPETRA_INST_COMPLEX_DOUBLE
+    template class PreconditionBase<std::complex<double>>;
+    template class PreconditionIdentity<std::complex<double>>;
+    template class PreconditionIfpackBase<std::complex<double>>;
+    template class PreconditionIfpack<std::complex<double>>;
+    template class PreconditionJacobi<std::complex<double>>;
+    template class PreconditionL1Jacobi<std::complex<double>>;
+    template class PreconditionL1GaussSeidel<std::complex<double>>;
+    template class PreconditionSOR<std::complex<double>>;
+    template class PreconditionSSOR<std::complex<double>>;
+    template class PreconditionBlockJacobi<std::complex<double>>;
+    template class PreconditionBlockSOR<std::complex<double>>;
+    template class PreconditionBlockSSOR<std::complex<double>>;
+    template class PreconditionChebyshev<std::complex<double>>;
+    template class PreconditionILU<std::complex<double>>;
+    template class PreconditionILUT<std::complex<double>>;
+#        endif
+#      endif
+
+#    endif
+
+  } // namespace TpetraWrappers
+
+} // namespace LinearAlgebra
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif // DEAL_II_TRILINOS_WITH_IFPACK2
+#endif   // DEAL_II_TRILINOS_WITH_TPETRA
diff --git a/tests/trilinos_tpetra/add_matrices_06.cc b/tests/trilinos_tpetra/add_matrices_06.cc
new file mode 100644 (file)
index 0000000..180281b
--- /dev/null
@@ -0,0 +1,127 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2004 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// check that SparseMatrix::add(SparseMatrix) and
+// SparseMatrix::copy_from(SparseMatrix) do not destroy memory when called
+// from matrices with the same sparsity pattern and thus, structures depending
+// on these matrices such as PreconditionJacobi do not get destroyed
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/trilinos_tpetra_precondition.h>
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+
+void
+test(LinearAlgebra::TpetraWrappers::SparseMatrix<double> &m)
+{
+  LinearAlgebra::TpetraWrappers::SparseMatrix<double> m2(m.m(), m.n(), 3U),
+    m3(m.m(), m.n(), 3U);
+
+  // first set a few entries one-by-one
+  for (unsigned int i = 0; i < m.m(); ++i)
+    for (unsigned int j = 0; j < m.n(); ++j)
+      if ((i + 2 * j + 1) % 3 == 0 || i == j)
+        {
+          m.set(i, j, i * j * .5 + .5);
+          m2.set(i, j, 1.);
+          m3.set(i, j, -0.1);
+        }
+
+  m.compress(VectorOperation::insert);
+  m2.compress(VectorOperation::insert);
+  m3.compress(VectorOperation::insert);
+
+  deallog << "Matrix nonzeros: " << m.n_nonzero_elements() << ' '
+          << m2.n_nonzero_elements() << ' ' << m3.n_nonzero_elements()
+          << std::endl;
+
+  m.copy_from(m2);
+  m.add(0.12, m3);
+
+  Vector<double> src(m.m()), dst(m.m());
+  src = 1.;
+
+  LinearAlgebra::TpetraWrappers::PreconditionJacobi<double> prec;
+  prec.initialize(m);
+
+  m.vmult(dst, src);
+  deallog << "Vector norm: " << dst.l2_norm() << ' ';
+
+  prec.vmult(dst, src);
+  deallog << dst.l2_norm() << std::endl;
+
+  m.copy_from(m2);
+  m.add(0.06, m3);
+
+  m.vmult(dst, src);
+  deallog << "Vector norm: " << dst.l2_norm() << ' ';
+
+  prec.vmult(dst, src);
+  deallog << dst.l2_norm() << std::endl;
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  try
+    {
+      {
+        LinearAlgebra::TpetraWrappers::SparseMatrix<double> m(5U, 5U, 3U);
+
+        test(m);
+      }
+    }
+  catch (const std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos_tpetra/add_matrices_06.output b/tests/trilinos_tpetra/add_matrices_06.output
new file mode 100644 (file)
index 0000000..f1e4081
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Matrix nonzeros: 13 13 13
+DEAL::Vector norm: 5.84509 2.26323
+DEAL::Vector norm: 5.88058 2.26323
+DEAL::OK
index 04bf70831f7acf175f2ed4ce101f3a11303900b1..24f71ee43f43bbbe1cfad1d4d51633326f088acc 100644 (file)
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/tria.h>
 
-#include "deal.II/lac/exceptions.h"
 #include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/full_matrix.h>
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/sparsity_tools.h>
-// #include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_tpetra_precondition.h>
 #include <deal.II/lac/trilinos_tpetra_solver_direct.h>
 #include <deal.II/lac/vector.h>
 
@@ -300,20 +299,16 @@ template <int dim>
 void
 Step4<dim>::solve()
 {
-  // TODO: implement precondition SSOR Wrappers
-  // // Compute 'reference' solution with CG solver and SSOR preconditioner
-  // TrilinosWrappers::PreconditionSSOR preconditioner;
-  // preconditioner.initialize(system_matrix);
+  // Compute 'reference' solution with CG solver and SSOR preconditioner
+  LinearAlgebra::TpetraWrappers::PreconditionSSOR<double> preconditioner;
+  preconditioner.initialize(system_matrix);
   LinearAlgebra::TpetraWrappers::Vector<double> temp_solution(system_rhs);
   temp_solution = 0;
   SolverControl solver_control(1000, 1e-12);
   SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
     solver_control);
 
-  solver.solve(system_matrix,
-               temp_solution,
-               system_rhs,
-               PreconditionIdentity());
+  solver.solve(system_matrix, temp_solution, system_rhs, preconditioner);
 
   constraints.distribute(temp_solution);
   solution = temp_solution;
@@ -323,10 +318,7 @@ Step4<dim>::solve()
   // do CG solve for new rhs
   temp_solution = 0;
   solution      = 0;
-  solver.solve(system_matrix,
-               temp_solution,
-               system_rhs_two,
-               PreconditionIdentity());
+  solver.solve(system_matrix, temp_solution, system_rhs_two, preconditioner);
 
   constraints.distribute(temp_solution);
   solution = temp_solution;
index 1facc1d50ab350a0d8252ba99bbe536514207c40..dff97d268c27645111f0ba2b621ad61af1d00514 100644 (file)
@@ -1,8 +1,8 @@
 
 DEAL:cg::Starting value 21.8027
-DEAL:cg::Convergence step 109 value 0
+DEAL:cg::Convergence step 86 value 0
 DEAL:cg::Starting value 0.0940512
-DEAL:cg::Convergence step 98 value 0
+DEAL:cg::Convergence step 77 value 0
 DEAL:DirectKLU2::Starting value 0
 DEAL:DirectKLU2::Convergence step 0 value 0
 DEAL:DirectKLU2::Norm of error in direct solve 1: 0
index 5e01f256793205a6de69a5828dfac949078ae170..e384b2fe0bc603a44895abb908ea25842801f524 100644 (file)
@@ -35,7 +35,7 @@
 #include <deal.II/lac/precondition.h>
 #include <deal.II/lac/solver_cg.h>
 #include <deal.II/lac/sparsity_tools.h>
-// #include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_tpetra_precondition.h>
 #include <deal.II/lac/trilinos_tpetra_solver_direct.h>
 #include <deal.II/lac/vector.h>
 
@@ -298,20 +298,16 @@ template <int dim>
 void
 Step4<dim>::solve()
 {
-  // TODO: implement precondition SSOR Wrappers
-  // // Compute 'reference' solution with CG solver and SSOR preconditioner
-  // TrilinosWrappers::PreconditionSSOR preconditioner;
-  // preconditioner.initialize(system_matrix);
+  // Compute 'reference' solution with CG solver and SSOR preconditioner
+  LinearAlgebra::TpetraWrappers::PreconditionSSOR<double> preconditioner;
+  preconditioner.initialize(system_matrix);
   LinearAlgebra::TpetraWrappers::Vector<double> temp_solution(system_rhs);
   temp_solution = 0;
   SolverControl solver_control(1000, 1e-12);
   SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
     solver_control);
 
-  solver.solve(system_matrix,
-               temp_solution,
-               system_rhs,
-               PreconditionIdentity());
+  solver.solve(system_matrix, temp_solution, system_rhs, preconditioner);
 
   constraints.distribute(temp_solution);
   solution = temp_solution;
@@ -321,10 +317,7 @@ Step4<dim>::solve()
   // do CG solve for new rhs
   temp_solution = 0;
   solution      = 0;
-  solver.solve(system_matrix,
-               temp_solution,
-               system_rhs_two,
-               PreconditionIdentity());
+  solver.solve(system_matrix, temp_solution, system_rhs_two, preconditioner);
 
   constraints.distribute(temp_solution);
   solution = temp_solution;
index 2f8d30713d8c24287dd71e0c245679357f2640a2..12fce67a02642dfaa121eacc550306f816325f0b 100644 (file)
@@ -1,8 +1,8 @@
 
 DEAL:cg::Starting value 21.8027
-DEAL:cg::Convergence step 109 value 0
+DEAL:cg::Convergence step 86 value 0
 DEAL:cg::Starting value 0.0940512
-DEAL:cg::Convergence step 98 value 0
+DEAL:cg::Convergence step 77 value 0
 DEAL:DirectKLU2::Starting value 0
 DEAL:DirectKLU2::Convergence step 0 value 0
 DEAL:DirectKLU2::Norm of error in direct solve 1: 0
diff --git a/tests/trilinos_tpetra/precondition.cc b/tests/trilinos_tpetra/precondition.cc
new file mode 100644 (file)
index 0000000..efbdcc8
--- /dev/null
@@ -0,0 +1,584 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2013 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// solves a 2D Poisson equation for linear elements with the simple
+// (Ifpack-based) Trilinos preconditioners
+
+#include "deal.II/base/exceptions.h"
+#include "deal.II/base/logstream.h"
+#include <deal.II/base/function.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/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include "deal.II/lac/trilinos_tpetra_vector.h"
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/solver_bicgstab.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/trilinos_tpetra_precondition.h>
+#include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <Teuchos_ParameterList.hpp>
+
+#include "../tests.h"
+
+
+
+template <int dim>
+class Step4
+{
+public:
+  Step4();
+  void
+  run();
+
+private:
+  void
+  make_grid();
+  void
+  setup_system();
+  void
+  assemble_system();
+  void
+  solve(int cycle);
+
+  Triangulation<dim> triangulation;
+  FE_Q<dim>          fe;
+  DoFHandler<dim>    dof_handler;
+
+  AffineConstraints<double> constraints;
+
+  LinearAlgebra::TpetraWrappers::SparseMatrix<double> system_matrix;
+
+  LinearAlgebra::TpetraWrappers::Vector<double> solution;
+  LinearAlgebra::TpetraWrappers::Vector<double> system_rhs;
+};
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+public:
+  RightHandSide()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  BoundaryValues()
+    : Function<dim>()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value(const Point<dim> &p,
+                          const unsigned int /*component*/) const
+{
+  double return_value = 0;
+  for (unsigned int i = 0; i < dim; ++i)
+    return_value += 4 * std::pow(p[i], 4);
+
+  return return_value;
+}
+
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value(const Point<dim> &p,
+                           const unsigned int /*component*/) const
+{
+  return p.square();
+}
+
+
+
+template <int dim>
+Step4<dim>::Step4()
+  : fe(1)
+  , dof_handler(triangulation)
+{}
+
+
+template <int dim>
+void
+Step4<dim>::make_grid()
+{
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+  triangulation.refine_global(5);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  constraints.clear();
+  std::map<unsigned int, double> boundary_values;
+  VectorTools::interpolate_boundary_values(dof_handler,
+                                           0,
+                                           BoundaryValues<dim>(),
+                                           constraints);
+  constraints.close();
+
+  DynamicSparsityPattern c_sparsity(dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, c_sparsity, constraints, false);
+  system_matrix.reinit(c_sparsity);
+  IndexSet is(dof_handler.n_dofs());
+  is.add_range(0, dof_handler.n_dofs());
+  solution.reinit(is);
+  system_rhs.reinit(is);
+}
+
+
+template <int dim>
+void
+Step4<dim>::assemble_system()
+{
+  QGauss<dim> quadrature_formula(fe.degree + 1);
+
+  const RightHandSide<dim> right_hand_side;
+
+  FEValues<dim> fe_values(fe,
+                          quadrature_formula,
+                          update_values | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quadrature_formula.size();
+
+  FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+  Vector<double>     cell_rhs(dofs_per_cell);
+
+  std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+
+  for (; cell != endc; ++cell)
+    {
+      fe_values.reinit(cell);
+      cell_matrix = 0;
+      cell_rhs    = 0;
+
+      for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              cell_matrix(i, j) +=
+                (fe_values.shape_grad(i, q_point) *
+                 fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point));
+
+            cell_rhs(i) +=
+              (fe_values.shape_value(i, q_point) *
+               right_hand_side.value(fe_values.quadrature_point(q_point)) *
+               fe_values.JxW(q_point));
+          }
+
+      cell->get_dof_indices(local_dof_indices);
+      constraints.distribute_local_to_global(
+        cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
+    }
+  system_matrix.compress(VectorOperation::add);
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::solve(int cycle)
+{
+  deallog.push(Utilities::int_to_string(dof_handler.n_dofs(), 5));
+
+  {
+    deallog.push("Identity");
+    static constexpr std::array<unsigned int, 2> lower{{49, 100}};
+    LinearAlgebra::TpetraWrappers::PreconditionIdentity<double> preconditioner;
+    solution = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    Teuchos::ParameterList precon_params;
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    // Make sure the general Ifpack preconditioner works at all.
+    // For this we use the point relaxation (SSOR)
+    // with default parameters, resulting in symmetric Gauss-Seidel (SGS)
+    deallog.push("CustomSGS");
+    static constexpr std::array<unsigned int, 2>              lower{{49, 100}};
+    LinearAlgebra::TpetraWrappers::PreconditionIfpack<double> preconditioner(
+      "RELAXATION");
+    solution = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    Teuchos::ParameterList precon_params;
+    precon_params.set("relaxation: type", "Symmetric Gauss-Seidel");
+    preconditioner.set_parameter_list(precon_params);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("Jacobi");
+    static constexpr std::array<unsigned int, 2>              lower{{49, 100}};
+    LinearAlgebra::TpetraWrappers::PreconditionJacobi<double> preconditioner;
+    solution = 0;
+
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("l1Jacobi");
+    static constexpr std::array<unsigned int, 2> lower{{49, 100}};
+    LinearAlgebra::TpetraWrappers::PreconditionL1Jacobi<double> preconditioner;
+    solution = 0;
+
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("l1GaussSeidel");
+    static constexpr std::array<unsigned int, 2> lower{{31, 62}};
+    LinearAlgebra::TpetraWrappers::PreconditionL1GaussSeidel<double>
+      preconditioner;
+    solution = 0;
+
+    SolverControl solver_control(1000, 1e-10);
+    SolverBicgstab<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("SOR");
+    static constexpr std::array<unsigned int, 2>           lower{{31, 62}};
+    LinearAlgebra::TpetraWrappers::PreconditionSOR<double> preconditioner;
+    solution = 0;
+
+    SolverControl solver_control(1000, 1e-5);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("SSOR");
+    static constexpr std::array<unsigned int, 2>            lower{{40, 77}};
+    LinearAlgebra::TpetraWrappers::PreconditionSSOR<double> preconditioner;
+    solution = 0;
+
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("BlockJacobi");
+    static constexpr std::array<unsigned int, 2> lower{{73, 145}};
+    LinearAlgebra::TpetraWrappers::PreconditionBlockJacobi<double>
+      preconditioner;
+    LinearAlgebra::TpetraWrappers::PreconditionBlockJacobi<
+      double>::AdditionalData data;
+    data.n_local_parts = 16;
+    solution           = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix, data);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("BlockSOR");
+    static constexpr std::array<unsigned int, 2>                lower{{18, 37}};
+    LinearAlgebra::TpetraWrappers::PreconditionBlockSOR<double> preconditioner;
+    LinearAlgebra::TpetraWrappers::PreconditionBlockSOR<double>::AdditionalData
+      data;
+    data.n_local_parts = 16;
+    data.omega         = 0.8;
+    solution           = 0;
+    SolverControl solver_control(1000, 1e-5);
+    SolverBicgstab<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix, data);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("BlockSSOR");
+    static constexpr std::array<unsigned int, 2> lower{{30, 59}};
+    LinearAlgebra::TpetraWrappers::PreconditionBlockSSOR<double> preconditioner;
+    LinearAlgebra::TpetraWrappers::PreconditionBlockSSOR<double>::AdditionalData
+      data;
+    data.n_local_parts = 16;
+    data.omega         = 1.2;
+    solution           = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix, data);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    static constexpr std::array<unsigned int, 2> lower{{30, 56}};
+    deallog.push("ILU");
+    LinearAlgebra::TpetraWrappers::PreconditionILU<double> preconditioner;
+    solution = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("ILUT");
+    static constexpr std::array<unsigned int, 2>            lower{{11, 19}};
+    LinearAlgebra::TpetraWrappers::PreconditionILUT<double> preconditioner;
+    LinearAlgebra::TpetraWrappers::PreconditionILUT<double>::AdditionalData
+      data;
+    data.ilut_drop = 1e-6;
+    data.ilut_fill = 3;
+    solution       = 0;
+    SolverControl solver_control(1000, 1e-5);
+    SolverBicgstab<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix, data);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("Chebyshev");
+    static constexpr std::array<unsigned int, 2> lower{{23, 46}};
+    LinearAlgebra::TpetraWrappers::PreconditionChebyshev<double> preconditioner;
+    LinearAlgebra::TpetraWrappers::PreconditionChebyshev<double>::AdditionalData
+      data;
+    data.max_eigenvalue = 2.5;
+    data.degree         = 3;
+    solution            = 0;
+    SolverControl solver_control(1000, 1e-10);
+    SolverCG<LinearAlgebra::TpetraWrappers::Vector<double>> solver(
+      solver_control);
+    preconditioner.initialize(system_matrix, data);
+    check_solver_within_range(
+      solver.solve(system_matrix, solution, system_rhs, preconditioner),
+      solver_control.last_step(),
+      lower[cycle],
+      lower[cycle] + 2);
+    deallog.pop();
+  }
+
+  //   {
+  //     deallog.push("Direct");
+  //     TrilinosWrappers::PreconditionBlockwiseDirect preconditioner;
+  //     solution = 0;
+  //     SolverControl solver_control(1000, 1e-10);
+  //     SolverCG<>    solver(solver_control);
+  //     preconditioner.initialize(system_matrix);
+  //     check_solver_within_range(
+  //       solver.solve(system_matrix, solution, system_rhs, preconditioner),
+  //       solver_control.last_step(),
+  //       1U,
+  //       1U);
+  //     deallog.pop();
+  //   }
+
+  {
+    // Make sure the general Ifpack preconditioner throws an Exception for an
+    // unsupported solver.
+    // IGLOO (incomplete global least optimal option) is of course a made up
+    // preconditioner name.
+
+    deallog.push("CustomIGLOO");
+    try
+      {
+        LinearAlgebra::TpetraWrappers::PreconditionIfpack<double>
+          preconditioner("IGLOO");
+      }
+    catch (dealii::ExceptionBase &exc)
+      {
+        deallog << "Error: " << std::endl;
+        exc.print_info(deallog.get_file_stream());
+      }
+    deallog.pop();
+  }
+
+  deallog.pop();
+}
+
+
+
+template <int dim>
+void
+Step4<dim>::run()
+{
+  for (unsigned int cycle = 0; cycle < 2; ++cycle)
+    {
+      if (cycle == 0)
+        make_grid();
+      else
+        triangulation.refine_global(1);
+
+      setup_system();
+      assemble_system();
+      solve(cycle);
+    }
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  try
+    {
+      Step4<2> test;
+      test.run();
+    }
+  catch (const std::exception &exc)
+    {
+      deallog << std::endl
+              << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl
+              << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/trilinos_tpetra/precondition.output b/tests/trilinos_tpetra/precondition.output
new file mode 100644 (file)
index 0000000..30f6004
--- /dev/null
@@ -0,0 +1,43 @@
+
+DEAL:01089:Identity::Solver stopped within 49 - 51 iterations
+DEAL:01089:CustomSGS::Solver stopped after 41 iterations
+DEAL:01089:Jacobi::Solver stopped within 49 - 51 iterations
+DEAL:01089:l1Jacobi::Solver stopped within 49 - 51 iterations
+DEAL:01089:l1GaussSeidel::Solver stopped after 48 iterations
+DEAL:01089:SOR::Solver stopped after 25 iterations
+DEAL:01089:SSOR::Solver stopped within 40 - 42 iterations
+DEAL:01089:BlockJacobi::Solver stopped after 58 iterations
+DEAL:01089:BlockSOR::Solver stopped after 16 iterations
+DEAL:01089:BlockSSOR::Solver stopped after 28 iterations
+DEAL:01089:ILU::Solver stopped within 30 - 32 iterations
+DEAL:01089:ILUT::Solver stopped after 5 iterations
+DEAL:01089:Chebyshev::Solver stopped after 38 iterations
+DEAL:01089:CustomIGLOO::Error: 
+    You tried to select the preconditioner type <IGLOO>
+but this preconditioner is not supported by Trilinos/Ifpack22
+due to one of the following reasons:
+* This preconditioner does not exist
+* This preconditioner has a specialized constructor not supported by the Ifpack2 Factory.
+* This preconditioner is not (yet) supported by Trilinos/Ifpack2
+* Trilinos/Ifpack2 was not configured for its use.
+DEAL:04225:Identity::Solver stopped within 100 - 102 iterations
+DEAL:04225:CustomSGS::Solver stopped after 78 iterations
+DEAL:04225:Jacobi::Solver stopped within 100 - 102 iterations
+DEAL:04225:l1Jacobi::Solver stopped within 100 - 102 iterations
+DEAL:04225:l1GaussSeidel::Solver stopped after 98 iterations
+DEAL:04225:SOR::Solver stopped after 46 iterations
+DEAL:04225:SSOR::Solver stopped within 77 - 79 iterations
+DEAL:04225:BlockJacobi::Solver stopped after 82 iterations
+DEAL:04225:BlockSOR::Solver stopped after 20 iterations
+DEAL:04225:BlockSSOR::Solver stopped after 40 iterations
+DEAL:04225:ILU::Solver stopped within 56 - 58 iterations
+DEAL:04225:ILUT::Solver stopped after 10 iterations
+DEAL:04225:Chebyshev::Solver stopped after 75 iterations
+DEAL:04225:CustomIGLOO::Error: 
+    You tried to select the preconditioner type <IGLOO>
+but this preconditioner is not supported by Trilinos/Ifpack22
+due to one of the following reasons:
+* This preconditioner does not exist
+* This preconditioner has a specialized constructor not supported by the Ifpack2 Factory.
+* This preconditioner is not (yet) supported by Trilinos/Ifpack2
+* Trilinos/Ifpack2 was not configured for its use.
\ No newline at end of file

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.