# 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
)
#
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}
+}
--- /dev/null
+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
#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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 ¶meter_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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 ¶meter_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
#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>
# 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
* local ordinate (processor local indices).
*/
using LO = int;
+
+
/**
* global ordinate (global indices).
*/
typename MemorySpace::kokkos_space>;
# endif
-
/**
* Communication between processors.
*/
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
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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+DEAL::Matrix nonzeros: 13 13 13
+DEAL::Vector norm: 5.84509 2.26323
+DEAL::Vector norm: 5.88058 2.26323
+DEAL::OK
#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>
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;
// 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;
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
#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>
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;
// 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;
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
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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;
+ };
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff