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 -------------------------*/