15#ifndef dealii_trilinos_solver_h
16#define dealii_trilinos_solver_h
21#ifdef DEAL_II_WITH_TRILINOS
33# include <Epetra_LinearProblem.h>
34# include <Epetra_Operator.h>
37# ifdef DEAL_II_TRILINOS_WITH_BELOS
38# include <BelosBlockCGSolMgr.hpp>
39# include <BelosBlockGmresSolMgr.hpp>
40# include <BelosEpetraAdapter.hpp>
41# include <BelosIteration.hpp>
42# include <BelosMultiVec.hpp>
43# include <BelosOperator.hpp>
44# include <BelosSolverManager.hpp>
57 class PreconditionBase;
211 Epetra_MultiVector &x,
212 const Epetra_MultiVector &b,
227 Epetra_MultiVector &x,
228 const Epetra_MultiVector &b,
246 const ::Vector<double> &b,
263 const ::Vector<double> &b,
275 const ::LinearAlgebra::distributed::Vector<double> &b,
288 const ::LinearAlgebra::distributed::Vector<double> &b,
303 <<
"An error with error number " << arg1
304 <<
" occurred while calling a Trilinos function");
320 template <
typename Preconditioner>
322 do_solve(
const Preconditioner &preconditioner);
327 template <
typename Preconditioner>
577 const ::LinearAlgebra::distributed::Vector<double> &b);
594 const ::LinearAlgebra::distributed::Vector<double> &b)
const;
615 const ::Vector<double> &b);
626 const ::LinearAlgebra::distributed::Vector<double> &b);
636 Epetra_MultiVector &x,
637 const Epetra_MultiVector &b);
661 <<
"An error with error number " << arg1
662 <<
" occurred while calling a Trilinos function");
696 std::unique_ptr<Amesos_BaseSolver>
solver;
706# ifdef DEAL_II_TRILINOS_WITH_BELOS
713 template <
typename VectorType>
739 struct AdditionalData
744 AdditionalData(
const SolverName solver_name = SolverName::cg,
745 const bool right_preconditioning =
false)
746 : solver_name(solver_name)
747 , right_preconditioning(right_preconditioning)
753 SolverName solver_name;
758 bool right_preconditioning;
765 const AdditionalData &additional_data,
766 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters);
771 template <
typename OperatorType,
typename PreconditionerType>
773 solve(
const OperatorType &a,
776 const PreconditionerType &p);
780 const AdditionalData additional_data;
781 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters;
791# ifdef DEAL_II_TRILINOS_WITH_BELOS
802 template <
typename VectorType>
804 class MultiVecWrapper
805 :
public Belos::MultiVec<typename VectorType::value_type>
811 using value_type =
typename VectorType::value_type;
817 this_type_is_missing_a_specialization()
825 MultiVecWrapper(VectorType &vector)
827 this->vectors.resize(1);
828 this->vectors[0].reset(
836 MultiVecWrapper(
const VectorType &vector)
838 this->vectors.resize(1);
839 this->vectors[0].reset(
847 virtual ~MultiVecWrapper() =
default;
852 virtual Belos::MultiVec<value_type> *
853 Clone(
const int numvecs)
const
855 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
857 new_multi_vec->vectors.resize(numvecs);
859 for (
auto &vec : new_multi_vec->vectors)
861 vec = std::make_shared<VectorType>();
864 vec->reinit(*this->vectors[0]);
867 return new_multi_vec;
873 virtual Belos::MultiVec<value_type> *
884 virtual Belos::MultiVec<value_type> *
885 CloneCopy(
const std::vector<int> &index)
const
887 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
889 new_multi_vec->vectors.resize(
index.size());
891 for (
unsigned int i = 0; i <
index.size(); ++i)
894 this->vectors.size(),
897 new_multi_vec->vectors[i] = std::make_shared<VectorType>();
900 *new_multi_vec->vectors[i] = *this->vectors[
index[i]];
903 return new_multi_vec;
911 virtual Belos::MultiVec<value_type> *
912 CloneViewNonConst(
const std::vector<int> &index)
914 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
916 new_multi_vec->vectors.resize(
index.size());
918 for (
unsigned int i = 0; i <
index.size(); ++i)
921 this->vectors.size(),
924 new_multi_vec->vectors[i].reset(
925 this->vectors[index[i]].get(),
931 return new_multi_vec;
939 virtual const Belos::MultiVec<value_type> *
940 CloneView(
const std::vector<int> &index)
const
942 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
944 new_multi_vec->vectors.resize(
index.size());
946 for (
unsigned int i = 0; i <
index.size(); ++i)
949 this->vectors.size(),
952 new_multi_vec->vectors[i].reset(
953 this->vectors[index[i]].get(),
959 return new_multi_vec;
966 GetGlobalLength()
const
970 for (
unsigned int i = 1; i < this->vectors.size(); ++i)
973 return this->vectors[0]->size();
980 GetNumberVecs()
const
982 return vectors.size();
989 MvTimesMatAddMv(
const value_type alpha,
990 const Belos::MultiVec<value_type> &A_,
991 const Teuchos::SerialDenseMatrix<int, value_type> &B,
992 const value_type beta)
994 const auto &
A = try_to_get_underlying_vector(A_);
996 const unsigned int n_rows = B.numRows();
997 const unsigned int n_cols = B.numCols();
999 AssertThrow(n_rows ==
static_cast<unsigned int>(
A.GetNumberVecs()),
1001 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1004 for (
unsigned int i = 0; i < n_cols; ++i)
1005 (*this->vectors[i]) *= beta;
1007 for (
unsigned int i = 0; i < n_cols; ++i)
1008 for (
unsigned int j = 0; j < n_rows; ++j)
1009 this->vectors[i]->add(alpha * B(j, i), *
A.vectors[j]);
1016 MvAddMv(
const value_type alpha,
1017 const Belos::MultiVec<value_type> &A_,
1018 const value_type beta,
1019 const Belos::MultiVec<value_type> &B_)
1021 const auto &
A = try_to_get_underlying_vector(A_);
1022 const auto &B = try_to_get_underlying_vector(B_);
1026 AssertThrow(this->vectors.size() == B.vectors.size(),
1029 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1031 this->vectors[i]->equ(alpha, *
A.vectors[i]);
1032 this->vectors[i]->add(beta, *B.vectors[i]);
1040 MvScale(
const value_type alpha)
1042 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1043 (*this->vectors[i]) *= alpha;
1050 MvScale(
const std::vector<value_type> &alpha)
1061 MvTransMv(
const value_type alpha,
1062 const Belos::MultiVec<value_type> &A_,
1063 Teuchos::SerialDenseMatrix<int, value_type> &B)
const
1065 const auto &
A = try_to_get_underlying_vector(A_);
1067 const unsigned int n_rows = B.numRows();
1068 const unsigned int n_cols = B.numCols();
1070 AssertThrow(n_rows ==
static_cast<unsigned int>(
A.GetNumberVecs()),
1072 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1075 for (
unsigned int i = 0; i < n_rows; ++i)
1076 for (
unsigned int j = 0; j < n_cols; ++j)
1077 B(i, j) = alpha * ((*
A.vectors[i]) * (*this->vectors[j]));
1085 MvDot(
const Belos::MultiVec<value_type> &A_,
1086 std::vector<value_type> &b)
const
1088 const auto &
A = try_to_get_underlying_vector(A_);
1094 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1095 b[i] = (*this->vectors[i]) * (*
A.vectors[i]);
1103 std::vector<
typename Teuchos::ScalarTraits<value_type>::magnitudeType>
1105 Belos::NormType type = Belos::TwoNorm)
const
1110 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1111 normvec[i] = this->vectors[i]->
l2_norm();
1118 SetBlock(
const Belos::MultiVec<value_type> &A,
1119 const std::vector<int> &index)
1139 MvInit(
const value_type alpha)
1149 MvPrint(std::ostream &os)
const
1170 genericVector()
const
1180 static MultiVecWrapper<VectorType> &
1181 try_to_get_underlying_vector(Belos::MultiVec<value_type> &vec_in)
1183 auto vec =
dynamic_cast<MultiVecWrapper<VectorType> *
>(&vec_in);
1193 const static MultiVecWrapper<VectorType> &
1194 try_to_get_underlying_vector(
const Belos::MultiVec<value_type> &vec_in)
1196 auto vec =
dynamic_cast<const MultiVecWrapper<VectorType> *
>(&vec_in);
1204# ifdef HAVE_BELOS_TSQR
1206 factorExplicit(Belos::MultiVec<value_type> &,
1207 Teuchos::SerialDenseMatrix<int, value_type> &,
1215 Teuchos::SerialDenseMatrix<int, value_type> &,
1216 const typename Teuchos::ScalarTraits<value_type>::magnitudeType &)
1227 std::vector<std::shared_ptr<VectorType>> vectors;
1233 MultiVecWrapper() =
default;
1242 template <
typename OperatorType,
typename VectorType>
1244 class OperatorWrapper
1245 :
public Belos::Operator<typename VectorType::value_type>
1251 using value_type =
typename VectorType::value_type;
1257 this_type_is_missing_a_specialization()
1265 OperatorWrapper(
const OperatorType &op)
1272 virtual ~OperatorWrapper() =
default;
1278 Apply(
const Belos::MultiVec<value_type> &x,
1279 Belos::MultiVec<value_type> &y,
1280 Belos::ETrans trans = Belos::NOTRANS)
const
1285 op.vmult(MultiVecWrapper<VectorType>::try_to_get_underlying_vector(y)
1287 MultiVecWrapper<VectorType>::try_to_get_underlying_vector(x)
1295 HasApplyTranspose()
const
1305 const OperatorType &op;
1312 template <
typename VectorType>
1314 SolverBelos<VectorType>::SolverBelos(
1316 const AdditionalData &additional_data,
1317 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters)
1318 : solver_control(solver_control)
1319 , additional_data(additional_data)
1320 , belos_parameters(belos_parameters)
1325 template <
typename VectorType>
1327 template <
typename OperatorType,
typename PreconditionerType>
1328 void SolverBelos<VectorType>::solve(
const OperatorType &A_dealii,
1329 VectorType &x_dealii,
1330 const VectorType &b_dealii,
1331 const PreconditionerType &P_dealii)
1333 using value_type =
typename VectorType::value_type;
1335 using MV = Belos::MultiVec<value_type>;
1336 using OP = Belos::Operator<value_type>;
1338 Teuchos::RCP<OP>
A = Teuchos::rcp(
1339 new internal::OperatorWrapper<OperatorType, VectorType>(A_dealii));
1340 Teuchos::RCP<OP> P = Teuchos::rcp(
1341 new internal::OperatorWrapper<PreconditionerType, VectorType>(P_dealii));
1342 Teuchos::RCP<MV> X =
1343 Teuchos::rcp(
new internal::MultiVecWrapper<VectorType>(x_dealii));
1344 Teuchos::RCP<MV> B =
1345 Teuchos::rcp(
new internal::MultiVecWrapper<VectorType>(b_dealii));
1347 Teuchos::RCP<Belos::LinearProblem<value_type, MV, OP>> problem =
1348 Teuchos::rcp(
new Belos::LinearProblem<value_type, MV, OP>(A, X, B));
1350 if (additional_data.right_preconditioning ==
false)
1351 problem->setLeftPrec(P);
1353 problem->setRightPrec(P);
1355 const bool problem_set = problem->setProblem();
1360 r.reinit(x_dealii,
true);
1361 A_dealii.vmult(r, x_dealii);
1362 r.sadd(-1., 1., b_dealii);
1363 const auto norm_0 = r.l2_norm();
1368 double relative_tolerance_to_be_achieved =
1369 solver_control.tolerance() / norm_0;
1370 const unsigned int max_steps = solver_control.max_steps();
1372 if (
const auto *reduction_control =
1374 relative_tolerance_to_be_achieved =
1375 std::max(relative_tolerance_to_be_achieved,
1376 reduction_control->reduction());
1378 Teuchos::RCP<Teuchos::ParameterList> belos_parameters_copy(
1379 Teuchos::rcp(
new Teuchos::ParameterList(*belos_parameters)));
1381 belos_parameters_copy->set(
"Convergence Tolerance",
1382 relative_tolerance_to_be_achieved);
1383 belos_parameters_copy->set(
"Maximum Iterations",
1384 static_cast<int>(max_steps));
1386 Teuchos::RCP<Belos::SolverManager<value_type, MV, OP>> solver;
1388 if (additional_data.solver_name == SolverName::cg)
1389 solver = Teuchos::rcp(
1390 new Belos::BlockCGSolMgr<value_type, MV, OP>(problem,
1391 belos_parameters_copy));
1392 else if (additional_data.solver_name == SolverName::gmres)
1393 solver = Teuchos::rcp(
1394 new Belos::BlockGmresSolMgr<value_type, MV, OP>(problem,
1395 belos_parameters_copy));
1399 const auto flag = solver->solve();
1401 solver_control.check(solver->getNumIters(), solver->achievedTol() * norm_0);
1403 AssertThrow(flag == Belos::ReturnType::Converged ||
1406 (solver_control.last_step() == max_steps)),
1408 solver_control.last_value()));
@ iterate
Continue iteration.
virtual ~SolverBase()=default
SolverControl & control() const
void solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b, const PreconditionBase &preconditioner)
const AdditionalData additional_data
std::unique_ptr< AztecOO_StatusTest > status_test
void set_preconditioner(AztecOO &solver, const Preconditioner &preconditioner)
SolverControl & solver_control
void do_solve(const Preconditioner &preconditioner)
void solve(Epetra_Operator &A, ::LinearAlgebra::distributed::Vector< double > &x, const ::LinearAlgebra::distributed::Vector< double > &b, const PreconditionBase &preconditioner)
std::unique_ptr< Epetra_LinearProblem > linear_problem
enum TrilinosWrappers::SolverBase::SolverName solver_name
void solve(const SparseMatrix &A, ::LinearAlgebra::distributed::Vector< double > &x, const ::LinearAlgebra::distributed::Vector< double > &b, const PreconditionBase &preconditioner)
void solve(const SparseMatrix &A, ::LinearAlgebra::distributed::Vector< double > &x, const ::LinearAlgebra::distributed::Vector< double > &b)
virtual ~SolverDirect()=default
void initialize(const SparseMatrix &A)
void vmult(MPI::Vector &x, const MPI::Vector &b) const
std::unique_ptr< Amesos_BaseSolver > solver
SolverControl & solver_control
void solve(MPI::Vector &x, const MPI::Vector &b)
std::unique_ptr< Epetra_LinearProblem > linear_problem
AdditionalData additional_data
SolverControl solver_control_own
SolverControl & control() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
#define DeclException1(Exception1, type1, outsequence)
#define AssertThrow(cond, exc)
std::vector< index_type > data
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const unsigned int gmres_restart_parameter
const bool output_solver_details
bool output_solver_details