From: Peter Munch <peterrmuench@gmail.com>
Date: Mon, 17 Oct 2022 07:50:59 +0000 (+0200)
Subject: Add TrilinosWrappers::SolverBelos
X-Git-Tag: v9.5.0-rc1~874^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b4f2b2833c532e9084ae378c9bb8ec13ab8a570f;p=dealii.git

Add TrilinosWrappers::SolverBelos
---

diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h
index c58a25a330..651c610cf6 100644
--- a/include/deal.II/lac/trilinos_solver.h
+++ b/include/deal.II/lac/trilinos_solver.h
@@ -26,11 +26,20 @@
 #  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
+#  include <BelosBlockGmresSolMgr.hpp>
+#  include <BelosEpetraAdapter.hpp>
+#  include <BelosIteration.hpp>
+#  include <BelosMultiVec.hpp>
+#  include <BelosOperator.hpp>
+#  include <BelosSolverManager.hpp>
+
 #  include <memory>
 
 
@@ -718,8 +727,469 @@ namespace TrilinosWrappers
     const AdditionalData additional_data;
   };
 
+
+
+  template <typename VectorType>
+  class SolverBelos
+  {
+  public:
+    SolverBelos(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:
+    const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters;
+  };
+
 } // namespace TrilinosWrappers
 
+
+
+#  ifndef DOXYGEN
+
+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 *) {});
+      }
+
+      MultiVecWrapper(const VectorType &vector)
+      {
+        this->vectors.resize(1);
+        this->vectors[0].reset(&const_cast<VectorType &>(vector),
+                               [](auto *) {});
+      }
+
+      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)
+          {
+            new_multi_vec->vectors[i]  = std::make_shared<VectorType>();
+            *new_multi_vec->vectors[i] = *this->vectors[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)
+          new_multi_vec->vectors[i].reset(this->vectors[index[i]].get(),
+                                          [](auto *) {});
+
+        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 *) {});
+          }
+
+        return new_multi_vec;
+      }
+
+      virtual ptrdiff_t
+      GetGlobalLength() const
+      {
+        AssertThrow(this->vectors.size() > 0, ExcInternalError());
+        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 = cast(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 = cast(A_);
+        const auto &B = cast(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 = cast(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)
+              {
+                AssertThrow(A.vectors[i], ExcNotImplemented());
+                AssertThrow(this->vectors[j], ExcNotImplemented());
+
+                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 = cast(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> &
+      cast(Belos::MultiVec<value_type> &vec_in)
+      {
+        auto vec = dynamic_cast<MultiVecWrapper<VectorType> *>(&vec_in);
+
+        AssertThrow(vec, ExcInternalError());
+
+        return *vec;
+      }
+
+      const static MultiVecWrapper<VectorType> &
+      cast(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> &                Q,
+                     Teuchos::SerialDenseMatrix<int, value_type> &R,
+                     const bool forceNonnegativeDiagonal = false)
+      {
+        TEUCHOS_TEST_FOR_EXCEPTION(
+          true,
+          std::logic_error,
+          "The Belos::MultiVec<"
+            << Teuchos::TypeNameTraits<value_type>::name()
+            << "> subclass which you "
+               "are using does not implement the TSQR-related method factorExplicit().");
+      }
+
+      virtual int
+      revealRank(
+        Teuchos::SerialDenseMatrix<int, value_type> &                    R,
+        const typename Teuchos::ScalarTraits<value_type>::magnitudeType &tol)
+      {
+        TEUCHOS_TEST_FOR_EXCEPTION(
+          true,
+          std::logic_error,
+          "The Belos::MultiVec<"
+            << Teuchos::TypeNameTraits<value_type>::name()
+            << "> subclass which you "
+               "are using does not implement the TSQR-related method revealRank().");
+      }
+
+#    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(){};
+
+      virtual void
+      Apply(const Belos::MultiVec<value_type> &x,
+            Belos::MultiVec<value_type> &      y,
+            Belos::ETrans                      trans = Belos::NOTRANS) const
+      {
+        AssertThrow(trans == Belos::NOTRANS, ExcNotImplemented());
+
+        op.vmult(MultiVecWrapper<VectorType>::cast(y).genericVector(),
+                 MultiVecWrapper<VectorType>::cast(x).genericVector());
+      }
+
+      virtual bool
+      HasApplyTranspose() const
+      {
+        return false;
+      }
+
+    private:
+      const OperatorType &op;
+    };
+
+  } // namespace internal
+
+
+
+  template <typename VectorType>
+  SolverBelos<VectorType>::SolverBelos(
+    const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters)
+    : belos_parameters(belos_parameters)
+  {}
+
+
+
+  template <typename VectorType>
+  template <typename OperatorType, typename PreconditionerType>
+  void
+  SolverBelos<VectorType>::solve(const OperatorType &      a,
+                                 VectorType &              x,
+                                 const VectorType &        b,
+                                 const PreconditionerType &p)
+  {
+    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));
+    Teuchos::RCP<OP> P = Teuchos::rcp(
+      new internal::OperatorWrapper<PreconditionerType, VectorType>(p));
+    Teuchos::RCP<MV> X =
+      Teuchos::rcp(new internal::MultiVecWrapper<VectorType>(x));
+    Teuchos::RCP<MV> B =
+      Teuchos::rcp(new internal::MultiVecWrapper<VectorType>(b));
+
+    Teuchos::RCP<Belos::LinearProblem<value_type, MV, OP>> problem =
+      Teuchos::rcp(new Belos::LinearProblem<value_type, MV, OP>(A, X, B));
+
+    if (true /*TODO*/)
+      problem->setLeftPrec(P);
+    else
+      problem->setRightPrec(P);
+
+    const bool problem_set = problem->setProblem();
+    AssertThrow(problem_set, ExcInternalError());
+
+    Teuchos::RCP<Belos::SolverManager<value_type, MV, OP>> solver =
+      Teuchos::rcp(
+        new Belos::BlockGmresSolMgr<value_type, MV, OP>(problem,
+                                                        belos_parameters));
+
+    const auto flag = solver->solve();
+
+    AssertThrow(flag == Belos::ReturnType::Converged, ExcInternalError());
+
+    unsigned int n_iterations = solver->getNumIters();
+
+    (void)n_iterations;
+  }
+
+} // namespace TrilinosWrappers
+
+#  endif
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif // DEAL_II_WITH_TRILINOS
diff --git a/tests/trilinos/solver_belos_01.cc b/tests/trilinos/solver_belos_01.cc
new file mode 100644
index 0000000000..dae643dac1
--- /dev/null
+++ b/tests/trilinos/solver_belos_01.cc
@@ -0,0 +1,156 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Check TrilinosWrappers::SolverBelos for GMRES.
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  using Number     = double;
+  using VectorType = Vector<Number>;
+
+  const unsigned int dim       = 2;
+  const unsigned int fe_degree = 1;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(3);
+
+  FE_Q<dim>      fe(fe_degree);
+  QGauss<dim>    quad(fe_degree + 1);
+  MappingQ1<dim> mapping;
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<Number> affine_constraints;
+  DoFTools::make_zero_boundary_constraints(dof_handler, affine_constraints);
+  affine_constraints.close();
+
+  TrilinosWrappers::SparsityPattern dsp(dof_handler.locally_owned_dofs(),
+                                        dof_handler.get_communicator());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp, affine_constraints);
+  dsp.compress();
+
+  TrilinosWrappers::SparseMatrix system_matrix;
+  system_matrix.reinit(dsp);
+
+  MatrixCreator::
+    create_laplace_matrix<dim, dim, TrilinosWrappers::SparseMatrix>(
+      mapping, dof_handler, quad, system_matrix, nullptr, affine_constraints);
+
+  TrilinosWrappers::PreconditionILU ilu;
+  ilu.initialize(system_matrix);
+
+  VectorType x(dof_handler.n_dofs());
+  VectorType r(dof_handler.n_dofs());
+
+  VectorTools::create_right_hand_side(mapping,
+                                      dof_handler,
+                                      quad,
+                                      Functions::ConstantFunction<dim, Number>(
+                                        1.0),
+                                      r,
+                                      affine_constraints);
+
+  Teuchos::RCP<Teuchos::ParameterList> belos_parameters =
+    Teuchos::rcp(new Teuchos::ParameterList);
+
+  belos_parameters->set("Num Blocks", 20);
+  belos_parameters->set("Block Size", 10);
+  belos_parameters->set("Verbosity", 0);
+
+  if (true)
+    {
+      x = 0.0;
+
+      using MV = Epetra_MultiVector;
+      using OP = Epetra_Operator;
+
+      Teuchos::RCP<Epetra_CrsMatrix> A =
+        Teuchos::rcp(const_cast<Epetra_CrsMatrix *>(
+                       &system_matrix.trilinos_matrix()),
+                     false);
+      Teuchos::RCP<Epetra_MultiVector> B, X;
+
+      LinearAlgebra::EpetraWrappers::Vector x_(dof_handler.locally_owned_dofs(),
+                                               dof_handler.get_communicator());
+      LinearAlgebra::ReadWriteVector<Number> x_temp(
+        dof_handler.locally_owned_dofs());
+      x_temp.import(x, VectorOperation::insert);
+      x_.import(x_temp, VectorOperation::insert);
+
+      LinearAlgebra::EpetraWrappers::Vector r_(dof_handler.locally_owned_dofs(),
+                                               dof_handler.get_communicator());
+      LinearAlgebra::ReadWriteVector<Number> r_temp(
+        dof_handler.locally_owned_dofs());
+      r_temp.import(r, VectorOperation::insert);
+      r_.import(r_temp, VectorOperation::insert);
+
+      X = Teuchos::rcp<Epetra_MultiVector>(&x_.trilinos_vector(), false);
+      B = Teuchos::rcp<Epetra_MultiVector>(&r_.trilinos_vector(), false);
+
+      Belos::LinearProblem<double, MV, OP> problem(A, X, B);
+      bool                                 set = problem.setProblem();
+
+      AssertThrow(set, ExcInternalError());
+
+      Teuchos::RCP<Belos::SolverManager<double, MV, OP>> newSolver =
+        Teuchos::rcp(new Belos::BlockGmresSolMgr<double, MV, OP>(
+          Teuchos::rcp(&problem, false), belos_parameters));
+      Belos::ReturnType flag = newSolver->solve();
+
+      AssertThrow(flag == Belos::ReturnType::Converged, ExcInternalError());
+
+      deallog << x_.l2_norm() << std::endl;
+    }
+
+  if (true)
+    {
+      x = 0.0;
+
+      TrilinosWrappers::SolverBelos<VectorType> solver(belos_parameters);
+      solver.solve(system_matrix, x, r, ilu);
+
+      deallog << x.l2_norm() << std::endl;
+    }
+}
diff --git a/tests/trilinos/solver_belos_01.with_trilinos=true.output b/tests/trilinos/solver_belos_01.with_trilinos=true.output
new file mode 100644
index 0000000000..9ff02fecb6
--- /dev/null
+++ b/tests/trilinos/solver_belos_01.with_trilinos=true.output
@@ -0,0 +1,3 @@
+
+DEAL::0.334271
+DEAL::0.334271