trilinos_precondition.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2019 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_trilinos_precondition_h
#  define dealii_trilinos_precondition_h
 
 
#  include <deal.II/base/config.h>
 
#  ifdef DEAL_II_WITH_TRILINOS
 
#    include <deal.II/base/subscriptor.h>
 
#    include <deal.II/lac/la_parallel_vector.h>
#    include <deal.II/lac/trilinos_vector.h>
 
#    include <memory>
 
#    ifdef DEAL_II_WITH_MPI
#      include <Epetra_MpiComm.h>
#    else
#      include <Epetra_SerialComm.h>
#    endif
#    include <Epetra_Map.h>
#    include <Epetra_MultiVector.h>
#    include <Epetra_RowMatrix.h>
#    include <Epetra_Vector.h>
#    include <Teuchos_ParameterList.hpp>
 
// forward declarations
class Ifpack_Preconditioner;
class Ifpack_Chebyshev;
namespace ML_Epetra
{
  class MultiLevelPreconditioner;
}
 
 
DEAL_II_NAMESPACE_OPEN
 
// forward declarations
template <typename number>
class SparseMatrix;
template <typename number>
class Vector;
class SparsityPattern;
 
namespace TrilinosWrappers
{
  // forward declarations
  class SparseMatrix;
  class BlockSparseMatrix;
  class SolverBase;
 
  class PreconditionBase : public Subscriptor
  {
  public:
    using size_type = ::types::global_dof_index;
 
    struct AdditionalData
    {};
 
    PreconditionBase();
 
    PreconditionBase(const PreconditionBase &);
 
    ~PreconditionBase() override = default;
 
    void
    clear();
 
    MPI_Comm
    get_mpi_communicator() const;
 
    void
    transpose();
 
    virtual void
    vmult(MPI::Vector &dst, const MPI::Vector &src) const;
 
    virtual void
    Tvmult(MPI::Vector &dst, const MPI::Vector &src) const;
 
    virtual void
    vmult(::Vector<double> &dst, const ::Vector<double> &src) const;
 
    virtual void
    Tvmult(::Vector<double> &      dst,
           const ::Vector<double> &src) const;
 
    virtual void
    vmult(::LinearAlgebra::distributed::Vector<double> &      dst,
          const ::LinearAlgebra::distributed::Vector<double> &src) const;
 
    virtual void
    Tvmult(::LinearAlgebra::distributed::Vector<double> &      dst,
           const ::LinearAlgebra::distributed::Vector<double> &src) const;
 
    Epetra_Operator &
    trilinos_operator() 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.");
 
    friend class SolverBase;
 
  protected:
    Teuchos::RCP<Epetra_Operator> preconditioner;
 
#    ifdef DEAL_II_WITH_MPI
    Epetra_MpiComm communicator;
#    else
    Epetra_SerialComm communicator;
#    endif
 
    std::shared_ptr<Epetra_Map> vector_distributor;
  };
 
 
  class PreconditionJacobi : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const double       omega        = 1,
                     const double       min_diagonal = 0,
                     const unsigned int n_sweeps     = 1);
 
      double omega;
 
      double min_diagonal;
 
      unsigned int n_sweeps;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionSSOR : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const double       omega        = 1,
                     const double       min_diagonal = 0,
                     const unsigned int overlap      = 0,
                     const unsigned int n_sweeps     = 1);
 
      double omega;
 
      double min_diagonal;
 
      unsigned int overlap;
 
      unsigned int n_sweeps;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionSOR : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const double       omega        = 1,
                     const double       min_diagonal = 0,
                     const unsigned int overlap      = 0,
                     const unsigned int n_sweeps     = 1);
 
      double omega;
 
      double min_diagonal;
 
      unsigned int overlap;
 
      unsigned int n_sweeps;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionBlockJacobi : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int block_size          = 1,
                     const std::string &block_creation_type = "linear",
                     const double       omega               = 1,
                     const double       min_diagonal        = 0,
                     const unsigned int n_sweeps            = 1);
 
      unsigned int block_size;
 
      std::string block_creation_type;
 
      double omega;
 
      double min_diagonal;
 
      unsigned int n_sweeps;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionBlockSSOR : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int block_size          = 1,
                     const std::string &block_creation_type = "linear",
                     const double       omega               = 1,
                     const double       min_diagonal        = 0,
                     const unsigned int overlap             = 0,
                     const unsigned int n_sweeps            = 1);
 
      unsigned int block_size;
 
      std::string block_creation_type;
 
      double omega;
 
      double min_diagonal;
 
      unsigned int overlap;
 
      unsigned int n_sweeps;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionBlockSOR : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int block_size          = 1,
                     const std::string &block_creation_type = "linear",
                     const double       omega               = 1,
                     const double       min_diagonal        = 0,
                     const unsigned int overlap             = 0,
                     const unsigned int n_sweeps            = 1);
 
      unsigned int block_size;
 
      std::string block_creation_type;
 
      double omega;
 
      double min_diagonal;
 
      unsigned int overlap;
 
      unsigned int n_sweeps;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionIC : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int ic_fill = 0,
                     const double       ic_atol = 0.,
                     const double       ic_rtol = 1.,
                     const unsigned int overlap = 0);
 
      unsigned int ic_fill;
 
      double ic_atol;
 
      double ic_rtol;
 
      unsigned int overlap;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionILU : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int ilu_fill = 0,
                     const double       ilu_atol = 0.,
                     const double       ilu_rtol = 1.,
                     const unsigned int overlap  = 0);
 
      unsigned int ilu_fill;
 
      double ilu_atol;
 
      double ilu_rtol;
 
      unsigned int overlap;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionILUT : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const double       ilut_drop = 0.,
                     const unsigned int ilut_fill = 0,
                     const double       ilut_atol = 0.,
                     const double       ilut_rtol = 1.,
                     const unsigned int overlap   = 0);
 
      double ilut_drop;
 
      unsigned int ilut_fill;
 
      double ilut_atol;
 
      double ilut_rtol;
 
      unsigned int overlap;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionBlockwiseDirect : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int overlap = 0);
 
 
      unsigned int overlap;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionChebyshev : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int degree           = 1,
                     const double       max_eigenvalue   = 10.,
                     const double       eigenvalue_ratio = 30.,
                     const double       min_eigenvalue   = 1.,
                     const double       min_diagonal     = 1e-12,
                     const bool         nonzero_starting = false);
 
      unsigned int degree;
 
      double max_eigenvalue;
 
      double eigenvalue_ratio;
 
      double min_eigenvalue;
 
      double min_diagonal;
 
      bool nonzero_starting;
    };
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
  };
 
 
 
  class PreconditionAMG : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const bool         elliptic              = true,
                     const bool         higher_order_elements = false,
                     const unsigned int n_cycles              = 1,
                     const bool         w_cyle                = false,
                     const double       aggregation_threshold = 1e-4,
                     const std::vector<std::vector<bool>> &constant_modes =
                       std::vector<std::vector<bool>>(0),
                     const unsigned int smoother_sweeps  = 2,
                     const unsigned int smoother_overlap = 0,
                     const bool         output_details   = false,
                     const char *       smoother_type    = "Chebyshev",
                     const char *       coarse_type      = "Amesos-KLU");
 
      void
      set_parameters(
        Teuchos::ParameterList &             parameter_list,
        std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
        const Epetra_RowMatrix &             matrix) const;
 
      void
      set_parameters(
        Teuchos::ParameterList &             parameter_list,
        std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
        const SparseMatrix &                 matrix) const;
 
      void
      set_operator_null_space(
        Teuchos::ParameterList &             parameter_list,
        std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
        const Epetra_RowMatrix &             matrix) const;
 
      void
      set_operator_null_space(
        Teuchos::ParameterList &             parameter_list,
        std::unique_ptr<Epetra_MultiVector> &distributed_constant_modes,
        const SparseMatrix &                 matrix) const;
 
      bool elliptic;
 
      bool higher_order_elements;
 
      unsigned int n_cycles;
 
      bool w_cycle;
 
      double aggregation_threshold;
 
      std::vector<std::vector<bool>> constant_modes;
 
      unsigned int smoother_sweeps;
 
      unsigned int smoother_overlap;
 
      bool output_details;
 
      const char *smoother_type;
 
      const char *coarse_type;
    };
 
    ~PreconditionAMG() override;
 
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const Epetra_RowMatrix &matrix,
               const AdditionalData &  additional_data = AdditionalData());
 
    void
    initialize(const SparseMatrix &          matrix,
               const Teuchos::ParameterList &ml_parameters);
 
    void
    initialize(const Epetra_RowMatrix &      matrix,
               const Teuchos::ParameterList &ml_parameters);
 
    template <typename number>
    void
    initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
               const AdditionalData &additional_data = AdditionalData(),
               const double          drop_tolerance  = 1e-13,
               const ::SparsityPattern *use_this_sparsity = nullptr);
 
    void
    reinit();
 
    void
    clear();
 
    size_type
    memory_consumption() const;
 
  private:
    std::shared_ptr<SparseMatrix> trilinos_matrix;
  };
 
 
 
#    if defined(DOXYGEN) || DEAL_II_TRILINOS_VERSION_GTE(11, 14, 0)
 
  class PreconditionAMGMueLu : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const bool         elliptic              = true,
                     const unsigned int n_cycles              = 1,
                     const bool         w_cyle                = false,
                     const double       aggregation_threshold = 1e-4,
                     const std::vector<std::vector<bool>> &constant_modes =
                       std::vector<std::vector<bool>>(0),
                     const unsigned int smoother_sweeps  = 2,
                     const unsigned int smoother_overlap = 0,
                     const bool         output_details   = false,
                     const char *       smoother_type    = "Chebyshev",
                     const char *       coarse_type      = "Amesos-KLU");
 
      bool elliptic;
 
      unsigned int n_cycles;
 
      bool w_cycle;
 
      double aggregation_threshold;
 
      std::vector<std::vector<bool>> constant_modes;
 
      unsigned int smoother_sweeps;
 
      unsigned int smoother_overlap;
 
      bool output_details;
 
      const char *smoother_type;
 
      const char *coarse_type;
    };
 
    PreconditionAMGMueLu();
 
    virtual ~PreconditionAMGMueLu() override = default;
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const Epetra_CrsMatrix &matrix,
               const AdditionalData &  additional_data = AdditionalData());
 
    void
    initialize(const SparseMatrix &    matrix,
               Teuchos::ParameterList &muelu_parameters);
 
    void
    initialize(const Epetra_CrsMatrix &matrix,
               Teuchos::ParameterList &muelu_parameters);
 
    template <typename number>
    void
    initialize(const ::SparseMatrix<number> &deal_ii_sparse_matrix,
               const AdditionalData &additional_data = AdditionalData(),
               const double          drop_tolerance  = 1e-13,
               const ::SparsityPattern *use_this_sparsity = nullptr);
 
    void
    clear();
 
    size_type
    memory_consumption() const;
 
  private:
    std::shared_ptr<SparseMatrix> trilinos_matrix;
  };
#    endif
 
 
 
  class PreconditionIdentity : public PreconditionBase
  {
  public:
    struct AdditionalData
    {};
 
    void
    initialize(const SparseMatrix &  matrix,
               const AdditionalData &additional_data = AdditionalData());
 
    void
    vmult(MPI::Vector &dst, const MPI::Vector &src) const override;
 
    void
    Tvmult(MPI::Vector &dst, const MPI::Vector &src) const override;
 
    void
    vmult(::Vector<double> &      dst,
          const ::Vector<double> &src) const override;
 
    void
    Tvmult(::Vector<double> &      dst,
           const ::Vector<double> &src) const override;
 
    void
    vmult(LinearAlgebra::distributed::Vector<double> &              dst,
          const ::LinearAlgebra::distributed::Vector<double> &src)
      const override;
 
    void
    Tvmult(LinearAlgebra::distributed::Vector<double> &              dst,
           const ::LinearAlgebra::distributed::Vector<double> &src)
      const override;
  };
 
 
 
  // ----------------------- inline and template functions --------------------
 
 
#    ifndef DOXYGEN
 
 
  inline void
  PreconditionBase::transpose()
  {
    // This only flips a flag that tells
    // Trilinos that any vmult operation
    // should be done with the
    // transpose. However, the matrix
    // structure is not reset.
    int ierr;
 
    if (!preconditioner->UseTranspose())
      {
        ierr = preconditioner->SetUseTranspose(true);
        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
      }
    else
      {
        ierr = preconditioner->SetUseTranspose(false);
        AssertThrow(ierr == 0, ExcTrilinosError(ierr));
      }
  }
 
 
  inline void
  PreconditionBase::vmult(MPI::Vector &dst, const MPI::Vector &src) const
  {
    Assert(dst.trilinos_partitioner().SameAs(
             preconditioner->OperatorRangeMap()),
           ExcNonMatchingMaps("dst"));
    Assert(src.trilinos_partitioner().SameAs(
             preconditioner->OperatorDomainMap()),
           ExcNonMatchingMaps("src"));
 
    const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
                                                  dst.trilinos_vector());
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
  }
 
  inline void
  PreconditionBase::Tvmult(MPI::Vector &dst, const MPI::Vector &src) const
  {
    Assert(dst.trilinos_partitioner().SameAs(
             preconditioner->OperatorRangeMap()),
           ExcNonMatchingMaps("dst"));
    Assert(src.trilinos_partitioner().SameAs(
             preconditioner->OperatorDomainMap()),
           ExcNonMatchingMaps("src"));
 
    preconditioner->SetUseTranspose(true);
    const int ierr = preconditioner->ApplyInverse(src.trilinos_vector(),
                                                  dst.trilinos_vector());
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
    preconditioner->SetUseTranspose(false);
  }
 
 
  // For the implementation of the <code>vmult</code> function with deal.II
  // data structures we note that invoking a call of the Trilinos
  // preconditioner requires us to use Epetra vectors as well. We do this by
  // providing a view, i.e., feed Trilinos with a pointer to the data, so we
  // avoid copying the content of the vectors during the iteration (this
  // function is only useful when used in serial anyway). In the declaration
  // of the right hand side, we need to cast the source vector (that is
  // <code>const</code> in all deal.II calls) to non-constant value, as this
  // is the way Trilinos wants to have them.
  inline void
  PreconditionBase::vmult(::Vector<double> &      dst,
                          const ::Vector<double> &src) const
  {
    AssertDimension(dst.size(),
                    preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension(src.size(),
                    preconditioner->OperatorRangeMap().NumMyElements());
    Epetra_Vector tril_dst(View,
                           preconditioner->OperatorDomainMap(),
                           dst.begin());
    Epetra_Vector tril_src(View,
                           preconditioner->OperatorRangeMap(),
                           const_cast<double *>(src.begin()));
 
    const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
  }
 
 
  inline void
  PreconditionBase::Tvmult(::Vector<double> &      dst,
                           const ::Vector<double> &src) const
  {
    AssertDimension(dst.size(),
                    preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension(src.size(),
                    preconditioner->OperatorRangeMap().NumMyElements());
    Epetra_Vector tril_dst(View,
                           preconditioner->OperatorDomainMap(),
                           dst.begin());
    Epetra_Vector tril_src(View,
                           preconditioner->OperatorRangeMap(),
                           const_cast<double *>(src.begin()));
 
    preconditioner->SetUseTranspose(true);
    const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
    preconditioner->SetUseTranspose(false);
  }
 
 
 
  inline void
  PreconditionBase::vmult(
    LinearAlgebra::distributed::Vector<double> &      dst,
    const LinearAlgebra::distributed::Vector<double> &src) const
  {
    AssertDimension(dst.local_size(),
                    preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension(src.local_size(),
                    preconditioner->OperatorRangeMap().NumMyElements());
    Epetra_Vector tril_dst(View,
                           preconditioner->OperatorDomainMap(),
                           dst.begin());
    Epetra_Vector tril_src(View,
                           preconditioner->OperatorRangeMap(),
                           const_cast<double *>(src.begin()));
 
    const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
  }
 
  inline void
  PreconditionBase::Tvmult(
    LinearAlgebra::distributed::Vector<double> &      dst,
    const LinearAlgebra::distributed::Vector<double> &src) const
  {
    AssertDimension(dst.local_size(),
                    preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension(src.local_size(),
                    preconditioner->OperatorRangeMap().NumMyElements());
    Epetra_Vector tril_dst(View,
                           preconditioner->OperatorDomainMap(),
                           dst.begin());
    Epetra_Vector tril_src(View,
                           preconditioner->OperatorRangeMap(),
                           const_cast<double *>(src.begin()));
 
    preconditioner->SetUseTranspose(true);
    const int ierr = preconditioner->ApplyInverse(tril_src, tril_dst);
    AssertThrow(ierr == 0, ExcTrilinosError(ierr));
    preconditioner->SetUseTranspose(false);
  }
 
#    endif
 
} // namespace TrilinosWrappers
 
 
DEAL_II_NAMESPACE_CLOSE
 
#  endif // DEAL_II_WITH_TRILINOS
 
#endif // trilinos_precondition_h
/*------------------------- trilinos_precondition.h -------------------------*/