From: Peter Munch 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 # include +// for AztecOO solvers # include # include # include # include +// for Belos solvers +# include +# include +# include +# include +# include +# include + # include @@ -718,8 +727,469 @@ namespace TrilinosWrappers const AdditionalData additional_data; }; + + + template + class SolverBelos + { + public: + SolverBelos(const Teuchos::RCP &belos_parameters); + + template + void + solve(const OperatorType & a, + VectorType & x, + const VectorType & b, + const PreconditionerType &p); + + private: + const Teuchos::RCP &belos_parameters; + }; + } // namespace TrilinosWrappers + + +# ifndef DOXYGEN + +namespace TrilinosWrappers +{ + namespace internal + { + template + class MultiVecWrapper + : public Belos::MultiVec + { + 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(vector), + [](auto *) {}); + } + + virtual ~MultiVecWrapper() = default; + + virtual Belos::MultiVec * + Clone(const int numvecs) const + { + auto new_multi_vec = new MultiVecWrapper; + + new_multi_vec->vectors.resize(numvecs); + + for (auto &vec : new_multi_vec->vectors) + { + vec = std::make_shared(); + + AssertThrow(this->vectors.size() > 0, ExcInternalError()); + vec->reinit(*this->vectors[0]); + } + + return new_multi_vec; + } + + virtual Belos::MultiVec * + CloneCopy() const + { + AssertThrow(false, ExcNotImplemented()); + } + + virtual Belos::MultiVec * + CloneCopy(const std::vector &index) const + { + auto new_multi_vec = new MultiVecWrapper; + + new_multi_vec->vectors.resize(index.size()); + + for (unsigned int i = 0; i < index.size(); ++i) + { + new_multi_vec->vectors[i] = std::make_shared(); + *new_multi_vec->vectors[i] = *this->vectors[i]; + } + + return new_multi_vec; + } + + virtual Belos::MultiVec * + CloneViewNonConst(const std::vector &index) + { + auto new_multi_vec = new MultiVecWrapper; + + 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 * + CloneView(const std::vector &index) const + { + auto new_multi_vec = new MultiVecWrapper; + + new_multi_vec->vectors.resize(index.size()); + + for (unsigned int i = 0; i < index.size(); ++i) + { + AssertThrow(static_cast(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 & A_, + const Teuchos::SerialDenseMatrix &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(A.GetNumberVecs()), + ExcInternalError()); + AssertThrow(n_cols == static_cast(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 &A_, + const value_type beta, + const Belos::MultiVec &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 &alpha) + { + AssertThrow(false, ExcNotImplemented()); + (void)alpha; + } + + virtual void + MvTransMv(const value_type alpha, + const Belos::MultiVec & A_, + Teuchos::SerialDenseMatrix &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(A.GetNumberVecs()), + ExcInternalError()); + AssertThrow(n_cols == static_cast(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 &A_, + std::vector & 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::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 &A, + const std::vector & 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 & + cast(Belos::MultiVec &vec_in) + { + auto vec = dynamic_cast *>(&vec_in); + + AssertThrow(vec, ExcInternalError()); + + return *vec; + } + + const static MultiVecWrapper & + cast(const Belos::MultiVec &vec_in) + { + auto vec = dynamic_cast *>(&vec_in); + + AssertThrow(vec, ExcInternalError()); + + return *vec; + } + + +# ifdef HAVE_BELOS_TSQR + virtual void + factorExplicit(Belos::MultiVec & Q, + Teuchos::SerialDenseMatrix &R, + const bool forceNonnegativeDiagonal = false) + { + TEUCHOS_TEST_FOR_EXCEPTION( + true, + std::logic_error, + "The Belos::MultiVec<" + << Teuchos::TypeNameTraits::name() + << "> subclass which you " + "are using does not implement the TSQR-related method factorExplicit()."); + } + + virtual int + revealRank( + Teuchos::SerialDenseMatrix & R, + const typename Teuchos::ScalarTraits::magnitudeType &tol) + { + TEUCHOS_TEST_FOR_EXCEPTION( + true, + std::logic_error, + "The Belos::MultiVec<" + << Teuchos::TypeNameTraits::name() + << "> subclass which you " + "are using does not implement the TSQR-related method revealRank()."); + } + +# endif + + private: + std::vector> vectors; + + MultiVecWrapper() = default; + }; + + template + class OperatorWrapper + : public Belos::Operator + { + 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 &x, + Belos::MultiVec & y, + Belos::ETrans trans = Belos::NOTRANS) const + { + AssertThrow(trans == Belos::NOTRANS, ExcNotImplemented()); + + op.vmult(MultiVecWrapper::cast(y).genericVector(), + MultiVecWrapper::cast(x).genericVector()); + } + + virtual bool + HasApplyTranspose() const + { + return false; + } + + private: + const OperatorType &op; + }; + + } // namespace internal + + + + template + SolverBelos::SolverBelos( + const Teuchos::RCP &belos_parameters) + : belos_parameters(belos_parameters) + {} + + + + template + template + void + SolverBelos::solve(const OperatorType & a, + VectorType & x, + const VectorType & b, + const PreconditionerType &p) + { + using value_type = typename VectorType::value_type; + + using MV = Belos::MultiVec; + using OP = Belos::Operator; + + Teuchos::RCP A = + Teuchos::rcp(new internal::OperatorWrapper(a)); + Teuchos::RCP P = Teuchos::rcp( + new internal::OperatorWrapper(p)); + Teuchos::RCP X = + Teuchos::rcp(new internal::MultiVecWrapper(x)); + Teuchos::RCP B = + Teuchos::rcp(new internal::MultiVecWrapper(b)); + + Teuchos::RCP> problem = + Teuchos::rcp(new Belos::LinearProblem(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> solver = + Teuchos::rcp( + new Belos::BlockGmresSolMgr(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 +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include + +#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; + + const unsigned int dim = 2; + const unsigned int fe_degree = 1; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(3); + + FE_Q fe(fe_degree); + QGauss quad(fe_degree + 1); + MappingQ1 mapping; + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints 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( + 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( + 1.0), + r, + affine_constraints); + + Teuchos::RCP 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 A = + Teuchos::rcp(const_cast( + &system_matrix.trilinos_matrix()), + false); + Teuchos::RCP B, X; + + LinearAlgebra::EpetraWrappers::Vector x_(dof_handler.locally_owned_dofs(), + dof_handler.get_communicator()); + LinearAlgebra::ReadWriteVector 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 r_temp( + dof_handler.locally_owned_dofs()); + r_temp.import(r, VectorOperation::insert); + r_.import(r_temp, VectorOperation::insert); + + X = Teuchos::rcp(&x_.trilinos_vector(), false); + B = Teuchos::rcp(&r_.trilinos_vector(), false); + + Belos::LinearProblem problem(A, X, B); + bool set = problem.setProblem(); + + AssertThrow(set, ExcInternalError()); + + Teuchos::RCP> newSolver = + Teuchos::rcp(new Belos::BlockGmresSolMgr( + 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 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