trilinos_precondition.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2018 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/trilinos_vector.h>
#  include <deal.II/lac/la_parallel_vector.h>

#  include <memory>

#  ifdef DEAL_II_WITH_MPI
#    include <Epetra_MpiComm.h>
#  else
#    include <Epetra_SerialComm.h>
#  endif
#  include <Epetra_Map.h>

#  include <Teuchos_ParameterList.hpp>
#  include <Epetra_RowMatrix.h>
#  include <Epetra_Vector.h>

// 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:
    typedef ::types::global_dof_index size_type;

    struct AdditionalData
    {};

    PreconditionBase ();

    PreconditionBase (const PreconditionBase &);

    ~PreconditionBase () = 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:
    std::shared_ptr<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");

      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();


    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();

    ~PreconditionAMGMueLu();

    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;

    void Tvmult (MPI::Vector       &dst,
                 const MPI::Vector &src) const;

    void vmult (::Vector<double>       &dst,
                const ::Vector<double> &src) const;

    void Tvmult (::Vector<double>       &dst,
                 const ::Vector<double> &src) const;

    void vmult (LinearAlgebra::distributed::Vector<double>       &dst,
                const ::LinearAlgebra::distributed::Vector<double> &src) const;

    void Tvmult (LinearAlgebra::distributed::Vector<double>       &dst,
                 const ::LinearAlgebra::distributed::Vector<double> &src) const;
  };



// -------------------------- 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.vector_partitioner().SameAs(preconditioner->OperatorRangeMap()),
            ExcNonMatchingMaps("dst"));
    Assert (src.vector_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.vector_partitioner().SameAs(preconditioner->OperatorRangeMap()),
            ExcNonMatchingMaps("dst"));
    Assert (src.vector_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 (static_cast<TrilinosWrappers::types::int_type>(dst.size()),
                     preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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 (static_cast<TrilinosWrappers::types::int_type>(dst.size()),
                     preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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 (static_cast<TrilinosWrappers::types::int_type>(dst.local_size()),
                     preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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 (static_cast<TrilinosWrappers::types::int_type>(dst.local_size()),
                     preconditioner->OperatorDomainMap().NumMyElements());
    AssertDimension (static_cast<TrilinosWrappers::types::int_type>(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

}


DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_TRILINOS

/*----------------------------   trilinos_precondition.h     ---------------------------*/

#endif
/*----------------------------   trilinos_precondition.h     ---------------------------*/