15#ifndef dealii_trilinos_solver_h
16#define dealii_trilinos_solver_h
21#ifdef DEAL_II_WITH_TRILINOS
35# include <Epetra_LinearProblem.h>
36# include <Epetra_Operator.h>
39# ifdef DEAL_II_TRILINOS_WITH_BELOS
40# include <BelosBlockCGSolMgr.hpp>
41# include <BelosBlockGmresSolMgr.hpp>
42# include <BelosEpetraAdapter.hpp>
43# include <BelosIteration.hpp>
44# include <BelosMultiVec.hpp>
45# include <BelosOperator.hpp>
46# include <BelosSolverManager.hpp>
62 class PreconditionBase;
216 Epetra_MultiVector &x,
217 const Epetra_MultiVector &b,
232 Epetra_MultiVector &x,
233 const Epetra_MultiVector &b,
251 const ::Vector<double> &b,
268 const ::Vector<double> &b,
280 const ::LinearAlgebra::distributed::Vector<double> &b,
293 const ::LinearAlgebra::distributed::Vector<double> &b,
308 <<
"An error with error number " << arg1
309 <<
" occurred while calling a Trilinos function");
325 template <
typename Preconditioner>
327 do_solve(
const Preconditioner &preconditioner);
332 template <
typename Preconditioner>
582 const ::LinearAlgebra::distributed::Vector<double> &b);
599 const ::LinearAlgebra::distributed::Vector<double> &b)
const;
620 const ::Vector<double> &b);
631 const ::LinearAlgebra::distributed::Vector<double> &b);
641 Epetra_MultiVector &x,
642 const Epetra_MultiVector &b);
666 <<
"An error with error number " << arg1
667 <<
" occurred while calling a Trilinos function");
701 std::unique_ptr<Amesos_BaseSolver>
solver;
711# ifdef DEAL_II_TRILINOS_WITH_BELOS
718 template <
typename VectorType>
744 struct AdditionalData
749 AdditionalData(
const SolverName solver_name = SolverName::cg,
750 const bool right_preconditioning =
false)
751 : solver_name(solver_name)
752 , right_preconditioning(right_preconditioning)
758 SolverName solver_name;
763 bool right_preconditioning;
770 const AdditionalData &additional_data,
771 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters);
776 template <
typename OperatorType,
typename PreconditionerType>
778 solve(
const OperatorType &a,
781 const PreconditionerType &p);
785 const AdditionalData additional_data;
786 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters;
796# ifdef DEAL_II_TRILINOS_WITH_BELOS
807 template <
typename VectorType>
809 class MultiVecWrapper
810 :
public Belos::MultiVec<typename VectorType::value_type>
816 using value_type =
typename VectorType::value_type;
822 this_type_is_missing_a_specialization()
830 MultiVecWrapper(VectorType &vector)
832 this->vectors.resize(1);
833 this->vectors[0].reset(
841 MultiVecWrapper(
const VectorType &vector)
843 this->vectors.resize(1);
844 this->vectors[0].reset(
852 virtual ~MultiVecWrapper() =
default;
857 virtual Belos::MultiVec<value_type> *
858 Clone(
const int numvecs)
const
860 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
862 new_multi_vec->vectors.resize(numvecs);
864 for (
auto &vec : new_multi_vec->vectors)
866 vec = std::make_shared<VectorType>();
869 vec->reinit(*this->vectors[0]);
872 return new_multi_vec;
878 virtual Belos::MultiVec<value_type> *
889 virtual Belos::MultiVec<value_type> *
890 CloneCopy(
const std::vector<int> &index)
const
892 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
894 new_multi_vec->vectors.resize(
index.size());
896 for (
unsigned int i = 0; i <
index.size(); ++i)
899 this->vectors.size(),
902 new_multi_vec->vectors[i] = std::make_shared<VectorType>();
905 *new_multi_vec->vectors[i] = *this->vectors[
index[i]];
908 return new_multi_vec;
916 virtual Belos::MultiVec<value_type> *
917 CloneViewNonConst(
const std::vector<int> &index)
919 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
921 new_multi_vec->vectors.resize(
index.size());
923 for (
unsigned int i = 0; i <
index.size(); ++i)
926 this->vectors.size(),
929 new_multi_vec->vectors[i].reset(
930 this->vectors[index[i]].get(),
936 return new_multi_vec;
944 virtual const Belos::MultiVec<value_type> *
945 CloneView(
const std::vector<int> &index)
const
947 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
949 new_multi_vec->vectors.resize(
index.size());
951 for (
unsigned int i = 0; i <
index.size(); ++i)
954 this->vectors.size(),
957 new_multi_vec->vectors[i].reset(
958 this->vectors[index[i]].get(),
964 return new_multi_vec;
970 virtual std::ptrdiff_t
971 GetGlobalLength()
const
975 for (
unsigned int i = 1; i < this->vectors.size(); ++i)
978 return this->vectors[0]->size();
985 GetNumberVecs()
const
987 return vectors.size();
994 MvTimesMatAddMv(
const value_type alpha,
995 const Belos::MultiVec<value_type> &A_,
996 const Teuchos::SerialDenseMatrix<int, value_type> &B,
997 const value_type beta)
999 const auto &
A = try_to_get_underlying_vector(A_);
1001 const unsigned int n_rows = B.numRows();
1002 const unsigned int n_cols = B.numCols();
1004 AssertThrow(n_rows ==
static_cast<unsigned int>(
A.GetNumberVecs()),
1006 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1009 for (
unsigned int i = 0; i < n_cols; ++i)
1010 (*this->vectors[i]) *= beta;
1012 for (
unsigned int i = 0; i < n_cols; ++i)
1013 for (
unsigned int j = 0; j < n_rows; ++j)
1014 this->vectors[i]->add(alpha * B(j, i), *
A.vectors[j]);
1021 MvAddMv(
const value_type alpha,
1022 const Belos::MultiVec<value_type> &A_,
1023 const value_type beta,
1024 const Belos::MultiVec<value_type> &B_)
1026 const auto &
A = try_to_get_underlying_vector(A_);
1027 const auto &B = try_to_get_underlying_vector(B_);
1031 AssertThrow(this->vectors.size() == B.vectors.size(),
1034 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1036 this->vectors[i]->equ(alpha, *
A.vectors[i]);
1037 this->vectors[i]->add(beta, *B.vectors[i]);
1045 MvScale(
const value_type alpha)
1047 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1048 (*this->vectors[i]) *= alpha;
1055 MvScale(
const std::vector<value_type> & )
1065 MvTransMv(
const value_type alpha,
1066 const Belos::MultiVec<value_type> &A_,
1067 Teuchos::SerialDenseMatrix<int, value_type> &B)
const
1069 const auto &
A = try_to_get_underlying_vector(A_);
1071 const unsigned int n_rows = B.numRows();
1072 const unsigned int n_cols = B.numCols();
1074 AssertThrow(n_rows ==
static_cast<unsigned int>(
A.GetNumberVecs()),
1076 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1079 for (
unsigned int i = 0; i < n_rows; ++i)
1080 for (
unsigned int j = 0; j < n_cols; ++j)
1081 B(i, j) = alpha * ((*
A.vectors[i]) * (*this->vectors[j]));
1089 MvDot(
const Belos::MultiVec<value_type> &A_,
1090 std::vector<value_type> &b)
const
1092 const auto &
A = try_to_get_underlying_vector(A_);
1098 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1099 b[i] = (*this->vectors[i]) * (*
A.vectors[i]);
1107 std::vector<
typename Teuchos::ScalarTraits<value_type>::magnitudeType>
1109 Belos::NormType type = Belos::TwoNorm)
const
1114 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1115 normvec[i] = this->vectors[i]->
l2_norm();
1122 SetBlock(
const Belos::MultiVec<value_type> & ,
1123 const std::vector<int> & )
1141 MvInit(
const value_type )
1150 MvPrint(std::ostream & )
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_DISABLE_EXTRA_DIAGNOSTICS
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
#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