From 98913f6a4fb21f4c45fcd65d360b624381d5cd1e Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 25 Oct 2022 13:14:55 +0200 Subject: [PATCH] Review changes --- include/deal.II/lac/trilinos_solver.h | 155 ++++++++++++++++---------- tests/trilinos/solver_belos_01.cc | 2 + 2 files changed, 96 insertions(+), 61 deletions(-) diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h index 6276035254..578dd39660 100644 --- a/include/deal.II/lac/trilinos_solver.h +++ b/include/deal.II/lac/trilinos_solver.h @@ -33,6 +33,7 @@ # include // for Belos solvers +# include # include # include # include @@ -739,6 +740,22 @@ namespace TrilinosWrappers class SolverBelos { public: + /** + * Enumeration object that Trilinos which solver to use. + * Currently enabled options are: + */ + enum SolverName + { + /** + * Use the conjugate gradient (CG) algorithm. + */ + cg, + /** + * Use the generalized minimum residual (GMRES) algorithm. + */ + gmres, + } solver_name; + /** * Struct that helps to configure SolverBelos. More advanced * parameters are passed to the constructor SolverBelos @@ -749,10 +766,16 @@ namespace TrilinosWrappers /** * Constructor. */ - AdditionalData(const bool right_preconditioning = false) + AdditionalData(const SolverName &solver_name = SolverName::cg, + const bool right_preconditioning = false) : right_preconditioning(right_preconditioning) {} + /** + * Solver type; + */ + SolverName solver_name; + /** * Flag for right preconditioning. */ @@ -888,8 +911,14 @@ namespace TrilinosWrappers 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]; + AssertThrow(static_cast(index[i]) < + this->vectors.size(), + ExcInternalError()); + + new_multi_vec->vectors[i] = std::make_shared(); + + AssertIndexRange(index[i], this->vectors.size()); + *new_multi_vec->vectors[i] = *this->vectors[index[i]]; } return new_multi_vec; @@ -908,10 +937,17 @@ namespace TrilinosWrappers 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 - *) { /*Nothing to do, since we are creating only a view.*/ }); + { + AssertThrow(static_cast(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; } @@ -951,6 +987,10 @@ namespace TrilinosWrappers 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(); } @@ -972,7 +1012,7 @@ namespace TrilinosWrappers const Teuchos::SerialDenseMatrix &B, const value_type beta) { - const auto &A = cast(A_); + const auto &A = try_to_get_underlying_vector(A_); const unsigned int n_rows = B.numRows(); const unsigned int n_cols = B.numCols(); @@ -999,8 +1039,8 @@ namespace TrilinosWrappers const value_type beta, const Belos::MultiVec &B_) { - const auto &A = cast(A_); - const auto &B = cast(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()); @@ -1043,7 +1083,7 @@ namespace TrilinosWrappers const Belos::MultiVec & A_, Teuchos::SerialDenseMatrix &B) const { - const auto &A = cast(A_); + const auto &A = try_to_get_underlying_vector(A_); const unsigned int n_rows = B.numRows(); const unsigned int n_cols = B.numCols(); @@ -1054,15 +1094,8 @@ namespace TrilinosWrappers 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])); - } - } + for (unsigned int j = 0; j < n_cols; ++j) + B(i, j) = alpha * ((*A.vectors[i]) * (*this->vectors[j])); } /** @@ -1073,7 +1106,7 @@ namespace TrilinosWrappers MvDot(const Belos::MultiVec &A_, std::vector & b) const { - const auto &A = cast(A_); + const auto &A = try_to_get_underlying_vector(A_); AssertThrow(this->vectors.size() == A.vectors.size(), ExcInternalError()); @@ -1166,7 +1199,7 @@ namespace TrilinosWrappers * Cast Belos::MultiVec to MultiVecWrapper. */ static MultiVecWrapper & - cast(Belos::MultiVec &vec_in) + try_to_get_underlying_vector(Belos::MultiVec &vec_in) { auto vec = dynamic_cast *>(&vec_in); @@ -1179,7 +1212,7 @@ namespace TrilinosWrappers * Cast Belos::MultiVec to MultiVecWrapper (const version). */ const static MultiVecWrapper & - cast(const Belos::MultiVec &vec_in) + try_to_get_underlying_vector(const Belos::MultiVec &vec_in) { auto vec = dynamic_cast *>(&vec_in); @@ -1191,31 +1224,19 @@ namespace TrilinosWrappers # ifdef HAVE_BELOS_TSQR virtual void - factorExplicit(Belos::MultiVec & Q, - Teuchos::SerialDenseMatrix &R, - const bool forceNonnegativeDiagonal = false) + factorExplicit(Belos::MultiVec &, + Teuchos::SerialDenseMatrix &, + const bool = 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()."); + Asser(false, ExcNotImplemented()); } virtual int revealRank( - Teuchos::SerialDenseMatrix & R, - const typename Teuchos::ScalarTraits::magnitudeType &tol) + Teuchos::SerialDenseMatrix &, + const typename Teuchos::ScalarTraits::magnitudeType &) { - 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()."); + Asser(false, ExcNotImplemented()); } # endif @@ -1262,7 +1283,8 @@ namespace TrilinosWrappers * Constructor. */ OperatorWrapper(const OperatorType &op) - : op(op){}; + : op(op) + {} /** * Destructor. @@ -1280,8 +1302,10 @@ namespace TrilinosWrappers // TODO: check for Tvmult AssertThrow(trans == Belos::NOTRANS, ExcNotImplemented()); - op.vmult(MultiVecWrapper::cast(y).genericVector(), - MultiVecWrapper::cast(x).genericVector()); + op.vmult(MultiVecWrapper::try_to_get_underlying_vector(y) + .genericVector(), + MultiVecWrapper::try_to_get_underlying_vector(x) + .genericVector()); } /** @@ -1320,24 +1344,24 @@ namespace TrilinosWrappers template template void - SolverBelos::solve(const OperatorType & a, - VectorType & x, - const VectorType & b, - const PreconditionerType &p) + SolverBelos::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; using OP = Belos::Operator; - Teuchos::RCP A = - Teuchos::rcp(new internal::OperatorWrapper(a)); + Teuchos::RCP A = Teuchos::rcp( + new internal::OperatorWrapper(A_dealii)); Teuchos::RCP P = Teuchos::rcp( - new internal::OperatorWrapper(p)); + new internal::OperatorWrapper(P_dealii)); Teuchos::RCP X = - Teuchos::rcp(new internal::MultiVecWrapper(x)); + Teuchos::rcp(new internal::MultiVecWrapper(x_dealii)); Teuchos::RCP B = - Teuchos::rcp(new internal::MultiVecWrapper(b)); + Teuchos::rcp(new internal::MultiVecWrapper(b_dealii)); Teuchos::RCP> problem = Teuchos::rcp(new Belos::LinearProblem(A, X, B)); @@ -1352,9 +1376,9 @@ namespace TrilinosWrappers // compute initial residal VectorType r; - r.reinit(x, true); - a.vmult(r, x); - r.sadd(-1., 1., b); + 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) @@ -1362,7 +1386,7 @@ namespace TrilinosWrappers double relative_tolerance_to_be_achieved = solver_control.tolerance() / norm_0; - const int max_steps = solver_control.max_steps(); + const unsigned int max_steps = solver_control.max_steps(); if (const auto *reduction_control = dynamic_cast(&solver_control)) @@ -1375,12 +1399,21 @@ namespace TrilinosWrappers belos_parameters_copy->set("Convergence Tolerance", relative_tolerance_to_be_achieved); - belos_parameters_copy->set("Maximum Iterations", max_steps); + belos_parameters_copy->set("Maximum Iterations", + static_cast(max_steps)); - Teuchos::RCP> solver = - Teuchos::rcp( + Teuchos::RCP> solver; + + if (additional_data.solver_name == SolverName::cg) + solver = Teuchos::rcp( + new Belos::BlockCGSolMgr(problem, + belos_parameters_copy)); + else if (additional_data.solver_name == SolverName::gmres) + solver = Teuchos::rcp( new Belos::BlockGmresSolMgr(problem, belos_parameters_copy)); + else + AssertThrow(false, ExcNotImplemented()); const auto flag = solver->solve(); diff --git a/tests/trilinos/solver_belos_01.cc b/tests/trilinos/solver_belos_01.cc index 439f5660f5..586687fe12 100644 --- a/tests/trilinos/solver_belos_01.cc +++ b/tests/trilinos/solver_belos_01.cc @@ -152,6 +152,8 @@ main(int argc, char *argv[]) typename TrilinosWrappers::SolverBelos::AdditionalData additional_data; + additional_data.solver_name = + TrilinosWrappers::SolverBelos::SolverName::gmres; additional_data.right_preconditioning = false; TrilinosWrappers::SolverBelos solver(solver_control, -- 2.39.5