# 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>
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+}
projects / dealii.git / commitdiff