trilinos_tpetra_precondition.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2024 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#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>
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
 
#  include <deal.II/base/enable_observer_pointer.h>
 
#  include <deal.II/lac/la_parallel_vector.h>
#  include <deal.II/lac/trilinos_tpetra_sparse_matrix.h>
 
#  include <Teuchos_BLAS_types.hpp>
#  include <Teuchos_ConfigDefs.hpp>
#  include <Teuchos_ParameterList.hpp>
#  include <Teuchos_RCPDecl.hpp>
#  include <Tpetra_Operator.hpp>
 
#endif // DEAL_II_TRILINOS_WITH_TPETRA
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_TRILINOS_WITH_TPETRA
namespace LinearAlgebra
{
  namespace TpetraWrappers
  {
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionBase : public EnableObserverPointer
    {
    public:
      struct AdditionalData
      {};
 
      PreconditionBase() = default;
 
      void
      clear();
 
      virtual void
      vmult(Vector<Number, MemorySpace>       &dst,
            const Vector<Number, MemorySpace> &src) const;
 
      virtual void
      Tvmult(Vector<Number, MemorySpace>       &dst,
             const Vector<Number, MemorySpace> &src) const;
 
      virtual void
      vmult(::Vector<Number> &dst, ::Vector<Number> &src) const;
 
      virtual void
      Tvmult(::Vector<Number> &dst, ::Vector<Number> &src) const;
 
      const TpetraTypes::LinearOperator<Number, MemorySpace> &
      trilinos_operator() const;
 
      Teuchos::RCP<TpetraTypes::LinearOperator<Number, MemorySpace>>
      trilinos_rcp() const;
 
      IndexSet
      locally_owned_domain_indices() const;
 
      IndexSet
      locally_owned_range_indices() const;
 
      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.");
 
      DeclExceptionMsg(
        ExcTransposeNotSupported,
        "The chosen preconditioner does not support transposing the matrix.");
    protected:
      Teuchos::RCP<TpetraTypes::LinearOperator<Number, MemorySpace>>
        preconditioner;
 
      Teuchos::ParameterList parameter_list;
    };
 
 
 
#  ifdef DEAL_II_TRILINOS_WITH_IFPACK2
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionIdentity : public PreconditionBase<Number, MemorySpace>
    {
    public:
      PreconditionIdentity() = default;
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A);
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionIfpackBase : public PreconditionBase<Number, MemorySpace>
    {
    public:
      PreconditionIfpackBase(const std::string &preconditioner_type);
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A);
 
      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:
      std::string preconditioner_type;
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionIfpack
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      PreconditionIfpack(const std::string &preconditioner_type);
 
      void
      set_parameter_list(Teuchos::ParameterList &parameter_list);
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionJacobi
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const double omega        = 1.,
                       const bool   fix_diagonal = false,
                       const double min_diagonal = 0.,
                       const int    n_sweeps     = 1);
 
        double omega;
 
        bool fix_diagonal;
 
        double min_diagonal;
 
        int n_sweeps;
      };
      PreconditionJacobi();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionL1Jacobi
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const double omega        = 1.,
                       const bool   fix_diagonal = false,
                       const double min_diagonal = 0.,
                       const int    n_sweeps     = 1);
 
        double omega;
 
        bool fix_diagonal;
 
        double min_diagonal;
 
        int n_sweeps;
      };
 
      PreconditionL1Jacobi();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionL1GaussSeidel
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        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);
 
        double omega;
 
        double eta;
 
        bool fix_diagonal;
 
        double min_diagonal;
 
        int n_sweeps;
      };
 
      PreconditionL1GaussSeidel();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionSOR : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const double omega        = 1.,
                       const int    overlap      = 0,
                       const bool   fix_diagonal = false,
                       const double min_diagonal = 0.,
                       const int    n_sweeps     = 1);
 
        double omega;
 
        int overlap;
 
        bool fix_diagonal;
 
        double min_diagonal;
 
        int n_sweeps;
      };
 
      PreconditionSOR();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionSSOR : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const double omega        = 1.,
                       const int    overlap      = 0,
                       const bool   fix_diagonal = false,
                       const double min_diagonal = 0.,
                       const int    n_sweeps     = 1);
 
        double omega;
 
        int overlap;
 
        bool fix_diagonal;
 
        double min_diagonal;
 
        int n_sweeps;
      };
 
 
      PreconditionSSOR();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionChebyshev
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        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);
        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;
 
        double eigenvalue_ratio;
 
        double min_diagonal;
 
        bool nonzero_starting;
      };
 
      PreconditionChebyshev();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionILU : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const int    ilu_fill = 0,
                       const double ilu_atol = 0.,
                       const double ilu_rtol = 1.,
                       const int    overlap  = 0);
 
        int ilu_fill;
 
        double ilu_atol;
 
        double ilu_rtol;
 
        int overlap;
      };
 
      PreconditionILU();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionILUT : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        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);
 
        double ilut_drop;
 
        double ilut_fill;
 
        double ilut_atol;
 
        double ilut_rtol;
 
        int overlap;
      };
 
 
 
      PreconditionILUT();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionBlockJacobi
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const int    n_local_parts = 1,
                       const double omega         = 1.,
                       const int    block_overlap = 0,
                       const int    n_sweeps      = 1);
 
        int n_local_parts;
 
        double omega;
 
        int block_overlap;
 
        int n_sweeps;
      };
 
 
      PreconditionBlockJacobi();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionBlockSOR
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const int    n_local_parts = 1,
                       const double omega         = 1,
                       const int    overlap       = 0,
                       const int    n_sweeps      = 1);
 
        int n_local_parts;
 
        double omega;
 
        int overlap;
 
        int n_sweeps;
      };
 
      PreconditionBlockSOR();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
 
 
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionBlockSSOR
      : public PreconditionIfpackBase<Number, MemorySpace>
    {
    public:
      struct AdditionalData
      {
        AdditionalData(const int    n_local_parts = 1,
                       const double omega         = 1,
                       const int    overlap       = 1,
                       const int    n_sweeps      = 1);
 
        int n_local_parts;
 
        double omega;
 
        int overlap;
 
        int n_sweeps;
      };
 
      PreconditionBlockSSOR();
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
#    ifdef DEAL_II_TRILINOS_WITH_TPETRA_MUELU
 
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionAMGMueLu : public PreconditionBase<Number, MemorySpace>
    {
    public:
      DeclException1(
        ExcTrilinosMueLuSmootherUnsupported,
        std::string,
        << "You tried to select the smoother type <" << arg1 << ">\n"
        << "but this is not supported by our interface\n"
        << "due to one of the following reasons:\n"
        << "* This smoother does not exist\n"
        << "* This smoother is not (yet) supported by our interface\n"
        << "* This smoother is not (yet) supported by MueLu.\n");
      DeclException1(
        ExcTrilinosMueLuCoarseSolverUnsupported,
        std::string,
        << "You tried to select the coarse solver type <" << arg1 << ">\n"
        << "but this is not supported by our interface\n"
        << "due to one of the following reasons:\n"
        << "* This coarse solver does not exist\n"
        << "* This coarse solver is not (yet) supported by our interface\n"
        << "* This coarse solver is not (yet) supported by MueLu.\n");
 
      struct AdditionalData
      {
        AdditionalData(const bool   elliptic              = true,
                       const bool   symmetric             = true,
                       const bool   w_cycle               = false,
                       const double aggregation_threshold = 1e-4,
                       const std::vector<std::vector<bool>> &constant_modes =
                         std::vector<std::vector<bool>>(0),
                       const int          smoother_sweeps  = 2,
                       const int          smoother_overlap = 0,
                       const bool         output_details   = false,
                       const std::string &smoother_type    = "Chebyshev",
                       const std::string &coarse_type      = "KLU2");
 
        bool elliptic;
 
        bool symmetric;
 
        bool w_cycle;
 
        double aggregation_threshold;
 
        std::vector<std::vector<bool>> constant_modes;
 
        int smoother_sweeps;
 
        int smoother_overlap;
 
        bool output_details;
 
        std::string smoother_type;
 
        std::string coarse_type;
      };
 
      PreconditionAMGMueLu() = default;
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const Teuchos::ParameterList            &parameters);
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
    };
#    endif
#  endif // DEAL_II_TRILINOS_WITH_IFPACK2
 
#  ifdef DEAL_II_TRILINOS_WITH_SHYLU_DDFROSCH
    template <typename Number, typename MemorySpace = ::MemorySpace::Host>
    class PreconditionFROSch : public PreconditionBase<Number, MemorySpace>
    {
    public:
      enum PreconditionerType
      {
        OneLevel,
        TwoLevel
      };
 
      enum SolverName
      {
        KLU,
        UMFPACK,
        MUMPS,
        SuperLU_dist
      };
 
      enum CombineMethod
      {
        Restricted,
        /*
         * Take the average on the overlap.
         */
        Averaging,
        /*
         * Add up the values on the overlap.
         */
        Full
      };
 
      enum CoarseType
      {
        IPOUHarmonicCoarseOperator,
        GDSWCoarseOperator,
        RGDSWCoarseOperator
      };
 
      struct AdditionalData
      {
        AdditionalData(
          const unsigned int       overlap                   = 1,
          const unsigned int       dimension                 = 2,
          const enum CombineMethod combine_values_in_overlap = Restricted,
          const enum SolverName    subdomain_solver          = KLU,
          const enum CoarseType    coarse_operator_type =
            IPOUHarmonicCoarseOperator,
          const enum SolverName extension_solver = KLU,
          const enum SolverName coarse_solver    = KLU);
 
        unsigned int overlap;
 
        unsigned int dimension;
 
        enum CombineMethod combine_values_in_overlap;
 
        enum SolverName subdomain_solver;
 
        enum CoarseType coarse_operator_type;
 
        enum SolverName extension_solver;
 
        enum SolverName coarse_solver;
      };
 
      PreconditionFROSch(
        const enum PreconditionerType precondition_type = OneLevel);
 
      void
      set_parameter_list(const Teuchos::ParameterList &parameter_list);
 
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const AdditionalData &additional_data = AdditionalData());
 
      template <int dim, int spacedim = dim>
      void
      initialize(const SparseMatrix<Number, MemorySpace> &A,
                 const DoFHandler<dim, spacedim>         &dof_handler,
                 const AdditionalData &additional_data = AdditionalData());
 
    private:
      const enum PreconditionerType precondition_type;
 
      bool user_provided_parameter_list;
    };
#  endif // DEAL_II_TRILINOS_WITH_SHYLU_DDFROSCH
 
  } // namespace TpetraWrappers
} // namespace LinearAlgebra
 
#endif // DEAL_II_TRILINOS_WITH_TPETRA
 
DEAL_II_NAMESPACE_CLOSE
 
#endif