13#ifndef dealii_trilinos_solver_h
14#define dealii_trilinos_solver_h
19#ifndef DEAL_II_TRILINOS_WITH_EPETRA
24#ifdef DEAL_II_TRILINOS_WITH_EPETRA
37# include <Epetra_LinearProblem.h>
38# include <Epetra_Operator.h>
41# ifdef DEAL_II_TRILINOS_WITH_BELOS
42# include <BelosBlockCGSolMgr.hpp>
43# include <BelosBlockGmresSolMgr.hpp>
44# include <BelosEpetraAdapter.hpp>
45# include <BelosIteration.hpp>
46# include <BelosMultiVec.hpp>
47# include <BelosOperator.hpp>
48# include <BelosSolverManager.hpp>
60#ifdef DEAL_II_TRILINOS_WITH_EPETRA
215 Epetra_MultiVector &x,
216 const Epetra_MultiVector &b,
231 Epetra_MultiVector &x,
232 const Epetra_MultiVector &b,
250 const ::Vector<double> &b,
267 const ::Vector<double> &b,
279 const ::LinearAlgebra::distributed::Vector<double> &b,
292 const ::LinearAlgebra::distributed::Vector<double> &b,
307 <<
"An error with error number " << arg1
308 <<
" occurred while calling a Trilinos function");
324 template <
typename Preconditioner>
326 do_solve(
const Preconditioner &preconditioner);
331 template <
typename Preconditioner>
581 const ::LinearAlgebra::distributed::Vector<double> &b);
598 const ::LinearAlgebra::distributed::Vector<double> &b)
const;
619 const ::Vector<double> &b);
630 const ::LinearAlgebra::distributed::Vector<double> &b);
640 Epetra_MultiVector &x,
641 const Epetra_MultiVector &b);
665 <<
"An error with error number " << arg1
666 <<
" occurred while calling a Trilinos function");
700 std::unique_ptr<Amesos_BaseSolver>
solver;
710# ifdef DEAL_II_TRILINOS_WITH_BELOS
717 template <
typename VectorType>
775 template <
typename OperatorType,
typename PreconditionerType>
780 const PreconditionerType &p);
795# ifdef DEAL_II_TRILINOS_WITH_BELOS
806 template <
typename VectorType>
808 class MultiVecWrapper
809 :
public Belos::MultiVec<typename VectorType::value_type>
815 using value_type =
typename VectorType::value_type;
821 this_type_is_missing_a_specialization()
829 MultiVecWrapper(VectorType &vector)
831 this->vectors.resize(1);
832 this->vectors[0].reset(
840 MultiVecWrapper(
const VectorType &vector)
842 this->vectors.resize(1);
843 this->vectors[0].reset(
844 &
const_cast<VectorType &
>(vector),
851 virtual ~MultiVecWrapper() =
default;
856 virtual Belos::MultiVec<value_type> *
857 Clone(
const int numvecs)
const
859 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
861 new_multi_vec->vectors.resize(numvecs);
863 for (
auto &vec : new_multi_vec->vectors)
865 vec = std::make_shared<VectorType>();
868 vec->reinit(*this->vectors[0]);
871 return new_multi_vec;
877 virtual Belos::MultiVec<value_type> *
888 virtual Belos::MultiVec<value_type> *
889 CloneCopy(
const std::vector<int> &index)
const
891 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
893 new_multi_vec->vectors.resize(index.size());
895 for (
unsigned int i = 0; i < index.size(); ++i)
898 this->vectors.size(),
901 new_multi_vec->vectors[i] = std::make_shared<VectorType>();
904 *new_multi_vec->vectors[i] = *this->vectors[index[i]];
907 return new_multi_vec;
915 virtual Belos::MultiVec<value_type> *
916 CloneViewNonConst(
const std::vector<int> &index)
918 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
920 new_multi_vec->vectors.resize(index.size());
922 for (
unsigned int i = 0; i < index.size(); ++i)
925 this->vectors.size(),
928 new_multi_vec->vectors[i].reset(
929 this->vectors[index[i]].get(),
935 return new_multi_vec;
943 virtual const Belos::MultiVec<value_type> *
944 CloneView(
const std::vector<int> &index)
const
946 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
948 new_multi_vec->vectors.resize(index.size());
950 for (
unsigned int i = 0; i < index.size(); ++i)
953 this->vectors.size(),
956 new_multi_vec->vectors[i].reset(
957 this->vectors[index[i]].get(),
963 return new_multi_vec;
969 virtual std::ptrdiff_t
970 GetGlobalLength()
const
974 for (
unsigned int i = 1; i < this->vectors.size(); ++i)
977 return this->vectors[0]->size();
984 GetNumberVecs()
const
986 return vectors.size();
994 const Belos::MultiVec<value_type> &A_,
995 const Teuchos::SerialDenseMatrix<int, value_type> &B,
998 const auto &A = try_to_get_underlying_vector(A_);
1000 const unsigned int n_rows = B.numRows();
1001 const unsigned int n_cols = B.numCols();
1003 AssertThrow(n_rows ==
static_cast<unsigned int>(A.GetNumberVecs()),
1005 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1008 for (
unsigned int i = 0; i < n_cols; ++i)
1009 (*this->vectors[i]) *= beta;
1011 for (
unsigned int i = 0; i < n_cols; ++i)
1012 for (
unsigned int j = 0; j < n_rows; ++j)
1013 this->vectors[i]->add(alpha * B(j, i), *A.vectors[j]);
1021 const Belos::MultiVec<value_type> &A_,
1023 const Belos::MultiVec<value_type> &B_)
1025 const auto &A = try_to_get_underlying_vector(A_);
1026 const auto &B = try_to_get_underlying_vector(B_);
1028 AssertThrow(this->vectors.size() == A.vectors.size(),
1030 AssertThrow(this->vectors.size() == B.vectors.size(),
1033 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1035 this->vectors[i]->equ(alpha, *A.vectors[i]);
1036 this->vectors[i]->add(beta, *B.vectors[i]);
1046 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1047 (*this->vectors[i]) *= alpha;
1054 MvScale(
const std::vector<value_type> & )
1065 const Belos::MultiVec<value_type> &A_,
1066 Teuchos::SerialDenseMatrix<int, value_type> &B)
const
1068 const auto &A = try_to_get_underlying_vector(A_);
1070 const unsigned int n_rows = B.numRows();
1071 const unsigned int n_cols = B.numCols();
1073 AssertThrow(n_rows ==
static_cast<unsigned int>(A.GetNumberVecs()),
1075 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1078 for (
unsigned int i = 0; i < n_rows; ++i)
1079 for (
unsigned int j = 0; j < n_cols; ++j)
1080 B(i, j) = alpha * ((*A.vectors[i]) * (*this->vectors[j]));
1088 MvDot(
const Belos::MultiVec<value_type> &A_,
1089 std::vector<value_type> &b)
const
1091 const auto &A = try_to_get_underlying_vector(A_);
1093 AssertThrow(this->vectors.size() == A.vectors.size(),
1097 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1098 b[i] = (*this->vectors[i]) * (*A.vectors[i]);
1106 std::vector<
typename Teuchos::ScalarTraits<value_type>::magnitudeType>
1108 Belos::NormType type = Belos::TwoNorm)
const
1113 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1114 normvec[i] = this->vectors[i]->l2_norm();
1121 SetBlock(
const Belos::MultiVec<value_type> & ,
1122 const std::vector<int> & )
1149 MvPrint(std::ostream & )
const
1169 genericVector()
const
1179 static MultiVecWrapper<VectorType> &
1180 try_to_get_underlying_vector(Belos::MultiVec<value_type> &vec_in)
1182 auto vec =
dynamic_cast<MultiVecWrapper<VectorType> *
>(&vec_in);
1192 const static MultiVecWrapper<VectorType> &
1193 try_to_get_underlying_vector(
const Belos::MultiVec<value_type> &vec_in)
1195 auto vec =
dynamic_cast<const MultiVecWrapper<VectorType> *
>(&vec_in);
1203# ifdef HAVE_BELOS_TSQR
1205 factorExplicit(Belos::MultiVec<value_type> &,
1206 Teuchos::SerialDenseMatrix<int, value_type> &,
1214 Teuchos::SerialDenseMatrix<int, value_type> &,
1215 const typename Teuchos::ScalarTraits<value_type>::magnitudeType &)
1226 std::vector<std::shared_ptr<VectorType>> vectors;
1232 MultiVecWrapper() =
default;
1241 template <
typename OperatorType,
typename VectorType>
1243 class OperatorWrapper
1244 :
public Belos::Operator<typename VectorType::value_type>
1250 using value_type =
typename VectorType::value_type;
1256 this_type_is_missing_a_specialization()
1264 OperatorWrapper(
const OperatorType &op)
1271 virtual ~OperatorWrapper() =
default;
1277 Apply(
const Belos::MultiVec<value_type> &x,
1278 Belos::MultiVec<value_type> &y,
1279 Belos::ETrans trans = Belos::NOTRANS)
const
1284 op.vmult(MultiVecWrapper<VectorType>::try_to_get_underlying_vector(y)
1286 MultiVecWrapper<VectorType>::try_to_get_underlying_vector(x)
1294 HasApplyTranspose()
const
1304 const OperatorType &op;
1311 template <
typename VectorType>
1314 SolverControl &solver_control,
1315 const AdditionalData &additional_data,
1316 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters)
1317 : solver_control(solver_control)
1318 , additional_data(additional_data)
1319 , belos_parameters(belos_parameters)
1324 template <
typename VectorType>
1326 template <
typename OperatorType,
typename PreconditionerType>
1328 VectorType &x_dealii,
1329 const VectorType &b_dealii,
1330 const PreconditionerType &P_dealii)
1332 using value_type =
typename VectorType::value_type;
1334 using MV = Belos::MultiVec<value_type>;
1335 using OP = Belos::Operator<value_type>;
1337 Teuchos::RCP<OP>
A = Teuchos::rcp(
1338 new internal::OperatorWrapper<OperatorType, VectorType>(A_dealii));
1339 Teuchos::RCP<OP> P = Teuchos::rcp(
1340 new internal::OperatorWrapper<PreconditionerType, VectorType>(P_dealii));
1341 Teuchos::RCP<MV> X =
1342 Teuchos::rcp(
new internal::MultiVecWrapper<VectorType>(x_dealii));
1343 Teuchos::RCP<MV> B =
1344 Teuchos::rcp(
new internal::MultiVecWrapper<VectorType>(b_dealii));
1346 Teuchos::RCP<Belos::LinearProblem<value_type, MV, OP>> problem =
1347 Teuchos::rcp(
new Belos::LinearProblem<value_type, MV, OP>(A, X, B));
1349 if (additional_data.right_preconditioning ==
false)
1350 problem->setLeftPrec(P);
1352 problem->setRightPrec(P);
1354 const bool problem_set = problem->setProblem();
1359 r.reinit(x_dealii,
true);
1360 A_dealii.vmult(r, x_dealii);
1361 r.sadd(-1., 1., b_dealii);
1362 const auto norm_0 = r.l2_norm();
1367 double relative_tolerance_to_be_achieved =
1368 solver_control.tolerance() / norm_0;
1369 const unsigned int max_steps = solver_control.max_steps();
1371 if (
const auto *reduction_control =
1372 dynamic_cast<ReductionControl *
>(&solver_control))
1373 relative_tolerance_to_be_achieved =
1374 std::max(relative_tolerance_to_be_achieved,
1375 reduction_control->reduction());
1377 Teuchos::RCP<Teuchos::ParameterList> belos_parameters_copy(
1378 Teuchos::rcp(
new Teuchos::ParameterList(*belos_parameters)));
1380 belos_parameters_copy->set(
"Convergence Tolerance",
1381 relative_tolerance_to_be_achieved);
1382 belos_parameters_copy->set(
"Maximum Iterations",
1383 static_cast<int>(max_steps));
1385 Teuchos::RCP<Belos::SolverManager<value_type, MV, OP>> solver;
1387 if (additional_data.solver_name == SolverName::cg)
1388 solver = Teuchos::rcp(
1389 new Belos::BlockCGSolMgr<value_type, MV, OP>(problem,
1390 belos_parameters_copy));
1391 else if (additional_data.solver_name == SolverName::gmres)
1392 solver = Teuchos::rcp(
1393 new Belos::BlockGmresSolMgr<value_type, MV, OP>(problem,
1394 belos_parameters_copy));
1398 const auto flag = solver->solve();
1400 solver_control.check(solver->getNumIters(), solver->achievedTol() * norm_0);
1402 AssertThrow(flag == Belos::ReturnType::Converged ||
1403 ((
dynamic_cast<IterationNumberControl *
>(&solver_control) !=
1405 (solver_control.last_step() == max_steps)),
1406 SolverControl::NoConvergence(solver_control.last_step(),
1407 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)
std::unique_ptr< Epetra_LinearProblem > linear_problem
SolverBase(SolverControl &cn, const AdditionalData &data=AdditionalData())
enum TrilinosWrappers::SolverBase::SolverName solver_name
SolverBelos(SolverControl &solver_control, const AdditionalData &additional_data, const Teuchos::RCP< Teuchos::ParameterList > &belos_parameters)
const AdditionalData additional_data
void solve(const OperatorType &a, VectorType &x, const VectorType &b, const PreconditionerType &p)
SolverControl & solver_control
const Teuchos::RCP< Teuchos::ParameterList > & belos_parameters
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
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
SolverDirect(const AdditionalData &data=AdditionalData())
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
#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)
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
const unsigned int gmres_restart_parameter
const bool output_solver_details
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
AdditionalData(const SolverName solver_name=SolverName::cg, const bool right_preconditioning=false)
bool right_preconditioning
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
bool output_solver_details