trilinos_solver.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2023 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_solver_h
#define dealii_trilinos_solver_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_TRILINOS
 
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/la_parallel_vector.h>
#  include <deal.II/lac/solver_control.h>
#  include <deal.II/lac/vector.h>
 
// for AztecOO solvers
#  include <Amesos.h>
#  include <AztecOO.h>
#  include <Epetra_LinearProblem.h>
#  include <Epetra_Operator.h>
 
// for Belos solvers
#  ifdef DEAL_II_TRILINOS_WITH_BELOS
#    include <BelosBlockCGSolMgr.hpp>
#    include <BelosBlockGmresSolMgr.hpp>
#    include <BelosEpetraAdapter.hpp>
#    include <BelosIteration.hpp>
#    include <BelosMultiVec.hpp>
#    include <BelosOperator.hpp>
#    include <BelosSolverManager.hpp>
#  endif
 
#  include <memory>
 
 
DEAL_II_NAMESPACE_OPEN
 
namespace TrilinosWrappers
{
  // forward declarations
#  ifndef DOXYGEN
  class SparseMatrix;
  class PreconditionBase;
#  endif
 
 
  class SolverBase
  {
  public:
    enum SolverName
    {
      cg,
      cgs,
      gmres,
      bicgstab,
      tfqmr
    } solver_name;
 
    struct AdditionalData
    {
      explicit AdditionalData(const bool         output_solver_details = false,
                              const unsigned int gmres_restart_parameter = 30);
 
      const bool output_solver_details;
 
      const unsigned int gmres_restart_parameter;
    };
 
    SolverBase(SolverControl &       cn,
               const AdditionalData &data = AdditionalData());
 
    SolverBase(const enum SolverName solver_name,
               SolverControl &       cn,
               const AdditionalData &data = AdditionalData());
 
    virtual ~SolverBase() = default;
 
    void
    solve(const SparseMatrix &    A,
          MPI::Vector &           x,
          const MPI::Vector &     b,
          const PreconditionBase &preconditioner);
 
    void
    solve(const Epetra_Operator & A,
          MPI::Vector &           x,
          const MPI::Vector &     b,
          const PreconditionBase &preconditioner);
 
    void
    solve(const Epetra_Operator &A,
          MPI::Vector &          x,
          const MPI::Vector &    b,
          const Epetra_Operator &preconditioner);
 
    void
    solve(const Epetra_Operator &   A,
          Epetra_MultiVector &      x,
          const Epetra_MultiVector &b,
          const PreconditionBase &  preconditioner);
 
    void
    solve(const Epetra_Operator &   A,
          Epetra_MultiVector &      x,
          const Epetra_MultiVector &b,
          const Epetra_Operator &   preconditioner);
 
 
 
    void
    solve(const SparseMatrix &          A,
          ::Vector<double> &      x,
          const ::Vector<double> &b,
          const PreconditionBase &      preconditioner);
 
    void
    solve(Epetra_Operator &             A,
          ::Vector<double> &      x,
          const ::Vector<double> &b,
          const PreconditionBase &      preconditioner);
 
    void
    solve(const SparseMatrix &                                      A,
          ::LinearAlgebra::distributed::Vector<double> &      x,
          const ::LinearAlgebra::distributed::Vector<double> &b,
          const PreconditionBase &preconditioner);
 
    void
    solve(Epetra_Operator &                                         A,
          ::LinearAlgebra::distributed::Vector<double> &      x,
          const ::LinearAlgebra::distributed::Vector<double> &b,
          const PreconditionBase &preconditioner);
 
 
    SolverControl &
    control() const;
 
    DeclException1(ExcTrilinosError,
                   int,
                   << "An error with error number " << arg1
                   << " occurred while calling a Trilinos function");
 
  protected:
    SolverControl &solver_control;
 
  private:
    template <typename Preconditioner>
    void
    do_solve(const Preconditioner &preconditioner);
 
    template <typename Preconditioner>
    void
    set_preconditioner(AztecOO &solver, const Preconditioner &preconditioner);
 
    std::unique_ptr<Epetra_LinearProblem> linear_problem;
 
    std::unique_ptr<AztecOO_StatusTest> status_test;
 
    AztecOO solver;
 
    const AdditionalData additional_data;
  };
 
 
  // provide a declaration for two explicit specializations
  template <>
  void
  SolverBase::set_preconditioner(AztecOO &               solver,
                                 const PreconditionBase &preconditioner);
 
  template <>
  void
  SolverBase::set_preconditioner(AztecOO &              solver,
                                 const Epetra_Operator &preconditioner);
 
 
  class SolverCG : public SolverBase
  {
  public:
    SolverCG(SolverControl &cn, const AdditionalData &data = AdditionalData());
  };
 
 
 
  class SolverCGS : public SolverBase
  {
  public:
    SolverCGS(SolverControl &cn, const AdditionalData &data = AdditionalData());
  };
 
 
 
  class SolverGMRES : public SolverBase
  {
  public:
    SolverGMRES(SolverControl &       cn,
                const AdditionalData &data = AdditionalData());
  };
 
 
 
  class SolverBicgstab : public SolverBase
  {
  public:
    SolverBicgstab(SolverControl &       cn,
                   const AdditionalData &data = AdditionalData());
  };
 
 
 
  class SolverTFQMR : public SolverBase
  {
  public:
    SolverTFQMR(SolverControl &       cn,
                const AdditionalData &data = AdditionalData());
  };
 
 
 
  class SolverDirect
  {
  public:
    struct AdditionalData
    {
      explicit AdditionalData(const bool         output_solver_details = false,
                              const std::string &solver_type = "Amesos_Klu");
 
      bool output_solver_details;
 
      std::string solver_type;
    };
 
    SolverDirect(SolverControl &       cn,
                 const AdditionalData &data = AdditionalData());
 
    virtual ~SolverDirect() = default;
 
    void
    initialize(const SparseMatrix &A);
 
    void
    solve(MPI::Vector &x, const MPI::Vector &b);
 
    void
    solve(::LinearAlgebra::distributed::Vector<double> &      x,
          const ::LinearAlgebra::distributed::Vector<double> &b);
 
    void
    solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b);
 
    void
    solve(const SparseMatrix &          A,
          ::Vector<double> &      x,
          const ::Vector<double> &b);
 
    void
    solve(const SparseMatrix &                                      A,
          ::LinearAlgebra::distributed::Vector<double> &      x,
          const ::LinearAlgebra::distributed::Vector<double> &b);
 
    SolverControl &
    control() const;
 
    DeclException1(ExcTrilinosError,
                   int,
                   << "An error with error number " << arg1
                   << " occurred while calling a Trilinos function");
 
  private:
    void
    do_solve();
 
    SolverControl &solver_control;
 
    std::unique_ptr<Epetra_LinearProblem> linear_problem;
 
    std::unique_ptr<Amesos_BaseSolver> solver;
 
    const AdditionalData additional_data;
  };
 
 
 
#  ifdef DEAL_II_TRILINOS_WITH_BELOS
  template <typename VectorType>
  class SolverBelos
  {
  public:
    enum SolverName
    {
      cg,
      gmres,
    } solver_name;
 
    struct AdditionalData
    {
      AdditionalData(const SolverName solver_name           = SolverName::cg,
                     const bool       right_preconditioning = false)
        : solver_name(solver_name)
        , right_preconditioning(right_preconditioning)
      {}
 
      SolverName solver_name;
 
      bool right_preconditioning;
    };
 
    SolverBelos(SolverControl &                             solver_control,
                const AdditionalData &                      additional_data,
                const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters);
 
    template <typename OperatorType, typename PreconditionerType>
    void
    solve(const OperatorType &      a,
          VectorType &              x,
          const VectorType &        b,
          const PreconditionerType &p);
 
  private:
    SolverControl &                             solver_control;
    const AdditionalData                        additional_data;
    const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters;
  };
#  endif
 
} // namespace TrilinosWrappers
 
 
 
#  ifndef DOXYGEN
 
#    ifdef DEAL_II_TRILINOS_WITH_BELOS
namespace TrilinosWrappers
{
  namespace internal
  {
    template <class VectorType>
    class MultiVecWrapper
      : public Belos::MultiVec<typename VectorType::value_type>
    {
    public:
      using value_type = typename VectorType::value_type;
 
      static bool
      this_type_is_missing_a_specialization()
      {
        return false;
      }
 
      MultiVecWrapper(VectorType &vector)
      {
        this->vectors.resize(1);
        this->vectors[0].reset(
          &vector,
          [](auto *) { /*Nothing to do, since vector is owned outside.*/ });
      }
 
      MultiVecWrapper(const VectorType &vector)
      {
        this->vectors.resize(1);
        this->vectors[0].reset(
          &const_cast<VectorType &>(vector),
          [](auto *) { /*Nothing to do, since vector is owned outside.*/ });
      }
 
      virtual ~MultiVecWrapper() = default;
 
      virtual Belos::MultiVec<value_type> *
      Clone(const int numvecs) const
      {
        auto new_multi_vec = new MultiVecWrapper<VectorType>;
 
        new_multi_vec->vectors.resize(numvecs);
 
        for (auto &vec : new_multi_vec->vectors)
          {
            vec = std::make_shared<VectorType>();
 
            AssertThrow(this->vectors.size() > 0, ExcInternalError());
            vec->reinit(*this->vectors[0]);
          }
 
        return new_multi_vec;
      }
 
      virtual Belos::MultiVec<value_type> *
      CloneCopy() const
      {
        AssertThrow(false, ExcNotImplemented());
      }
 
      virtual Belos::MultiVec<value_type> *
      CloneCopy(const std::vector<int> &index) const
      {
        auto new_multi_vec = new MultiVecWrapper<VectorType>;
 
        new_multi_vec->vectors.resize(index.size());
 
        for (unsigned int i = 0; i < index.size(); ++i)
          {
            AssertThrow(static_cast<unsigned int>(index[i]) <
                          this->vectors.size(),
                        ExcInternalError());
 
            new_multi_vec->vectors[i] = std::make_shared<VectorType>();
 
            AssertIndexRange(index[i], this->vectors.size());
            *new_multi_vec->vectors[i] = *this->vectors[index[i]];
          }
 
        return new_multi_vec;
      }
 
      virtual Belos::MultiVec<value_type> *
      CloneViewNonConst(const std::vector<int> &index)
      {
        auto new_multi_vec = new MultiVecWrapper<VectorType>;
 
        new_multi_vec->vectors.resize(index.size());
 
        for (unsigned int i = 0; i < index.size(); ++i)
          {
            AssertThrow(static_cast<unsigned int>(index[i]) <
                          this->vectors.size(),
                        ExcInternalError());
 
            new_multi_vec->vectors[i].reset(
              this->vectors[index[i]].get(),
              [](
                auto
                  *) { /*Nothing to do, since we are creating only a view.*/ });
          }
 
        return new_multi_vec;
      }
 
      virtual const Belos::MultiVec<value_type> *
      CloneView(const std::vector<int> &index) const
      {
        auto new_multi_vec = new MultiVecWrapper<VectorType>;
 
        new_multi_vec->vectors.resize(index.size());
 
        for (unsigned int i = 0; i < index.size(); ++i)
          {
            AssertThrow(static_cast<unsigned int>(index[i]) <
                          this->vectors.size(),
                        ExcInternalError());
 
            new_multi_vec->vectors[i].reset(
              this->vectors[index[i]].get(),
              [](
                auto
                  *) { /*Nothing to do, since we are creating only a view.*/ });
          }
 
        return new_multi_vec;
      }
 
      virtual ptrdiff_t
      GetGlobalLength() const
      {
        AssertThrow(this->vectors.size() > 0, ExcInternalError());
 
        for (unsigned int i = 1; i < this->vectors.size(); ++i)
          AssertDimension(this->vectors[0]->size(), this->vectors[i]->size());
 
        return this->vectors[0]->size();
      }
 
      virtual int
      GetNumberVecs() const
      {
        return vectors.size();
      }
 
      virtual void
      MvTimesMatAddMv(const value_type                                   alpha,
                      const Belos::MultiVec<value_type> &                A_,
                      const Teuchos::SerialDenseMatrix<int, value_type> &B,
                      const value_type                                   beta)
      {
        const auto &A = try_to_get_underlying_vector(A_);
 
        const unsigned int n_rows = B.numRows();
        const unsigned int n_cols = B.numCols();
 
        AssertThrow(n_rows == static_cast<unsigned int>(A.GetNumberVecs()),
                    ExcInternalError());
        AssertThrow(n_cols == static_cast<unsigned int>(this->GetNumberVecs()),
                    ExcInternalError());
 
        for (unsigned int i = 0; i < n_cols; ++i)
          (*this->vectors[i]) *= beta;
 
        for (unsigned int i = 0; i < n_cols; ++i)
          for (unsigned int j = 0; j < n_rows; ++j)
            this->vectors[i]->add(alpha * B(j, i), *A.vectors[j]);
      }
 
      virtual void
      MvAddMv(const value_type                   alpha,
              const Belos::MultiVec<value_type> &A_,
              const value_type                   beta,
              const Belos::MultiVec<value_type> &B_)
      {
        const auto &A = try_to_get_underlying_vector(A_);
        const auto &B = try_to_get_underlying_vector(B_);
 
        AssertThrow(this->vectors.size() == A.vectors.size(),
                    ExcInternalError());
        AssertThrow(this->vectors.size() == B.vectors.size(),
                    ExcInternalError());
 
        for (unsigned int i = 0; i < this->vectors.size(); ++i)
          {
            this->vectors[i]->equ(alpha, *A.vectors[i]);
            this->vectors[i]->add(beta, *B.vectors[i]);
          }
      }
 
      virtual void
      MvScale(const value_type alpha)
      {
        for (unsigned int i = 0; i < this->vectors.size(); ++i)
          (*this->vectors[i]) *= alpha;
      }
 
      virtual void
      MvScale(const std::vector<value_type> &alpha)
      {
        AssertThrow(false, ExcNotImplemented());
        (void)alpha;
      }
 
      virtual void
      MvTransMv(const value_type                             alpha,
                const Belos::MultiVec<value_type> &          A_,
                Teuchos::SerialDenseMatrix<int, value_type> &B) const
      {
        const auto &A = try_to_get_underlying_vector(A_);
 
        const unsigned int n_rows = B.numRows();
        const unsigned int n_cols = B.numCols();
 
        AssertThrow(n_rows == static_cast<unsigned int>(A.GetNumberVecs()),
                    ExcInternalError());
        AssertThrow(n_cols == static_cast<unsigned int>(this->GetNumberVecs()),
                    ExcInternalError());
 
        for (unsigned int i = 0; i < n_rows; ++i)
          for (unsigned int j = 0; j < n_cols; ++j)
            B(i, j) = alpha * ((*A.vectors[i]) * (*this->vectors[j]));
      }
 
      virtual void
      MvDot(const Belos::MultiVec<value_type> &A_,
            std::vector<value_type> &          b) const
      {
        const auto &A = try_to_get_underlying_vector(A_);
 
        AssertThrow(this->vectors.size() == A.vectors.size(),
                    ExcInternalError());
        AssertThrow(this->vectors.size() == b.size(), ExcInternalError());
 
        for (unsigned int i = 0; i < this->vectors.size(); ++i)
          b[i] = (*this->vectors[i]) * (*A.vectors[i]);
      }
 
      virtual void
      MvNorm(
        std::vector<typename Teuchos::ScalarTraits<value_type>::magnitudeType>
          &             normvec,
        Belos::NormType type = Belos::TwoNorm) const
      {
        AssertThrow(type == Belos::TwoNorm, ExcNotImplemented());
        AssertThrow(this->vectors.size() == normvec.size(), ExcInternalError());
 
        for (unsigned int i = 0; i < this->vectors.size(); ++i)
          normvec[i] = this->vectors[i]->l2_norm();
      }
 
      virtual void
      SetBlock(const Belos::MultiVec<value_type> &A,
               const std::vector<int> &           index)
      {
        AssertThrow(false, ExcNotImplemented());
        (void)A;
        (void)index;
      }
 
      virtual void
      MvRandom()
      {
        AssertThrow(false, ExcNotImplemented());
      }
 
      virtual void
      MvInit(const value_type alpha)
      {
        AssertThrow(false, ExcNotImplemented());
        (void)alpha;
      }
 
      virtual void
      MvPrint(std::ostream &os) const
      {
        AssertThrow(false, ExcNotImplemented());
        (void)os;
      }
 
      VectorType &
      genericVector()
      {
        AssertThrow(GetNumberVecs() == 1, ExcNotImplemented());
 
        return *vectors[0];
      }
 
      const VectorType &
      genericVector() const
      {
        AssertThrow(GetNumberVecs() == 1, ExcNotImplemented());
 
        return *vectors[0];
      }
 
      static MultiVecWrapper<VectorType> &
      try_to_get_underlying_vector(Belos::MultiVec<value_type> &vec_in)
      {
        auto vec = dynamic_cast<MultiVecWrapper<VectorType> *>(&vec_in);
 
        AssertThrow(vec, ExcInternalError());
 
        return *vec;
      }
 
      const static MultiVecWrapper<VectorType> &
      try_to_get_underlying_vector(const Belos::MultiVec<value_type> &vec_in)
      {
        auto vec = dynamic_cast<const MultiVecWrapper<VectorType> *>(&vec_in);
 
        AssertThrow(vec, ExcInternalError());
 
        return *vec;
      }
 
 
#      ifdef HAVE_BELOS_TSQR
      virtual void
      factorExplicit(Belos::MultiVec<value_type> &,
                     Teuchos::SerialDenseMatrix<int, value_type> &,
                     const bool = false)
      {
        Assert(false, ExcNotImplemented());
      }
 
      virtual int
      revealRank(
        Teuchos::SerialDenseMatrix<int, value_type> &,
        const typename Teuchos::ScalarTraits<value_type>::magnitudeType &)
      {
        Assert(false, ExcNotImplemented());
      }
 
#      endif
 
    private:
      std::vector<std::shared_ptr<VectorType>> vectors;
 
      MultiVecWrapper() = default;
    };
 
    template <class OperatorType, class VectorType>
    class OperatorWrapper
      : public Belos::Operator<typename VectorType::value_type>
    {
    public:
      using value_type = typename VectorType::value_type;
 
      static bool
      this_type_is_missing_a_specialization()
      {
        return false;
      }
 
      OperatorWrapper(const OperatorType &op)
        : op(op)
      {}
 
      virtual ~OperatorWrapper() = default;
 
      virtual void
      Apply(const Belos::MultiVec<value_type> &x,
            Belos::MultiVec<value_type> &      y,
            Belos::ETrans                      trans = Belos::NOTRANS) const
      {
        // TODO: check for Tvmult
        AssertThrow(trans == Belos::NOTRANS, ExcNotImplemented());
 
        op.vmult(MultiVecWrapper<VectorType>::try_to_get_underlying_vector(y)
                   .genericVector(),
                 MultiVecWrapper<VectorType>::try_to_get_underlying_vector(x)
                   .genericVector());
      }
 
      virtual bool
      HasApplyTranspose() const
      {
        // TODO: check for Tvmult
        return false;
      }
 
    private:
      const OperatorType &op;
    };
 
  } // namespace internal
 
 
 
  template <typename VectorType>
  SolverBelos<VectorType>::SolverBelos(
    SolverControl &                             solver_control,
    const AdditionalData &                      additional_data,
    const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters)
    : solver_control(solver_control)
    , additional_data(additional_data)
    , belos_parameters(belos_parameters)
  {}
 
 
 
  template <typename VectorType>
  template <typename OperatorType, typename PreconditionerType>
  void
  SolverBelos<VectorType>::solve(const OperatorType &      A_dealii,
                                 VectorType &              x_dealii,
                                 const VectorType &        b_dealii,
                                 const PreconditionerType &P_dealii)
  {
    using value_type = typename VectorType::value_type;
 
    using MV = Belos::MultiVec<value_type>;
    using OP = Belos::Operator<value_type>;
 
    Teuchos::RCP<OP> A = Teuchos::rcp(
      new internal::OperatorWrapper<OperatorType, VectorType>(A_dealii));
    Teuchos::RCP<OP> P = Teuchos::rcp(
      new internal::OperatorWrapper<PreconditionerType, VectorType>(P_dealii));
    Teuchos::RCP<MV> X =
      Teuchos::rcp(new internal::MultiVecWrapper<VectorType>(x_dealii));
    Teuchos::RCP<MV> B =
      Teuchos::rcp(new internal::MultiVecWrapper<VectorType>(b_dealii));
 
    Teuchos::RCP<Belos::LinearProblem<value_type, MV, OP>> problem =
      Teuchos::rcp(new Belos::LinearProblem<value_type, MV, OP>(A, X, B));
 
    if (additional_data.right_preconditioning == false)
      problem->setLeftPrec(P);
    else
      problem->setRightPrec(P);
 
    const bool problem_set = problem->setProblem();
    AssertThrow(problem_set, ExcInternalError());
 
    // compute initial residal
    VectorType r;
    r.reinit(x_dealii, true);
    A_dealii.vmult(r, x_dealii);
    r.sadd(-1., 1., b_dealii);
    const auto norm_0 = r.l2_norm();
 
    if (solver_control.check(0, norm_0) != SolverControl::iterate)
      return;
 
    double relative_tolerance_to_be_achieved =
      solver_control.tolerance() / norm_0;
    const unsigned int max_steps = solver_control.max_steps();
 
    if (const auto *reduction_control =
          dynamic_cast<ReductionControl *>(&solver_control))
      relative_tolerance_to_be_achieved =
        std::max(relative_tolerance_to_be_achieved,
                 reduction_control->reduction());
 
    Teuchos::RCP<Teuchos::ParameterList> belos_parameters_copy(
      Teuchos::rcp(new Teuchos::ParameterList(*belos_parameters)));
 
    belos_parameters_copy->set("Convergence Tolerance",
                               relative_tolerance_to_be_achieved);
    belos_parameters_copy->set("Maximum Iterations",
                               static_cast<int>(max_steps));
 
    Teuchos::RCP<Belos::SolverManager<value_type, MV, OP>> solver;
 
    if (additional_data.solver_name == SolverName::cg)
      solver = Teuchos::rcp(
        new Belos::BlockCGSolMgr<value_type, MV, OP>(problem,
                                                     belos_parameters_copy));
    else if (additional_data.solver_name == SolverName::gmres)
      solver = Teuchos::rcp(
        new Belos::BlockGmresSolMgr<value_type, MV, OP>(problem,
                                                        belos_parameters_copy));
    else
      AssertThrow(false, ExcNotImplemented());
 
    const auto flag = solver->solve();
 
    solver_control.check(solver->getNumIters(), solver->achievedTol() * norm_0);
 
    AssertThrow(flag == Belos::ReturnType::Converged ||
                  ((dynamic_cast<IterationNumberControl *>(&solver_control) !=
                    nullptr) &&
                   (solver_control.last_step() == max_steps)),
                SolverControl::NoConvergence(solver_control.last_step(),
                                             solver_control.last_value()));
  }
 
} // namespace TrilinosWrappers
#    endif
 
#  endif
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_TRILINOS
 
#endif