15#ifndef dealii_trilinos_solver_h
16#define dealii_trilinos_solver_h
21#ifdef DEAL_II_WITH_TRILINOS
31# include <Epetra_LinearProblem.h>
32# include <Epetra_Operator.h>
35# ifdef DEAL_II_TRILINOS_WITH_BELOS
36# include <BelosBlockCGSolMgr.hpp>
37# include <BelosBlockGmresSolMgr.hpp>
38# include <BelosEpetraAdapter.hpp>
39# include <BelosIteration.hpp>
40# include <BelosMultiVec.hpp>
41# include <BelosOperator.hpp>
42# include <BelosSolverManager.hpp>
55 class PreconditionBase;
209 Epetra_MultiVector &x,
210 const Epetra_MultiVector &b,
225 Epetra_MultiVector &x,
226 const Epetra_MultiVector &b,
244 const ::Vector<double> &b,
261 const ::Vector<double> &b,
273 const ::LinearAlgebra::distributed::Vector<double> &b,
286 const ::LinearAlgebra::distributed::Vector<double> &b,
301 <<
"An error with error number " << arg1
302 <<
" occurred while calling a Trilinos function");
318 template <
typename Preconditioner>
320 do_solve(
const Preconditioner &preconditioner);
325 template <
typename Preconditioner>
575 const ::LinearAlgebra::distributed::Vector<double> &b);
592 const ::LinearAlgebra::distributed::Vector<double> &b)
const;
613 const ::Vector<double> &b);
624 const ::LinearAlgebra::distributed::Vector<double> &b);
637 <<
"An error with error number " << arg1
638 <<
" occurred while calling a Trilinos function");
672 std::unique_ptr<Amesos_BaseSolver>
solver;
682# ifdef DEAL_II_TRILINOS_WITH_BELOS
689 template <
typename VectorType>
714 struct AdditionalData
719 AdditionalData(
const SolverName solver_name = SolverName::cg,
720 const bool right_preconditioning =
false)
721 : solver_name(solver_name)
722 , right_preconditioning(right_preconditioning)
728 SolverName solver_name;
733 bool right_preconditioning;
740 const AdditionalData &additional_data,
741 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters);
746 template <
typename OperatorType,
typename PreconditionerType>
748 solve(
const OperatorType &a,
751 const PreconditionerType &p);
755 const AdditionalData additional_data;
756 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters;
766# ifdef DEAL_II_TRILINOS_WITH_BELOS
777 template <
typename VectorType>
778 class MultiVecWrapper
779 :
public Belos::MultiVec<typename VectorType::value_type>
785 using value_type =
typename VectorType::value_type;
791 this_type_is_missing_a_specialization()
799 MultiVecWrapper(VectorType &vector)
801 this->vectors.resize(1);
802 this->vectors[0].reset(
810 MultiVecWrapper(
const VectorType &vector)
812 this->vectors.resize(1);
813 this->vectors[0].reset(
814 &
const_cast<VectorType &
>(vector),
821 virtual ~MultiVecWrapper() =
default;
826 virtual Belos::MultiVec<value_type> *
827 Clone(
const int numvecs)
const
829 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
831 new_multi_vec->vectors.resize(numvecs);
833 for (
auto &vec : new_multi_vec->vectors)
835 vec = std::make_shared<VectorType>();
838 vec->reinit(*this->vectors[0]);
841 return new_multi_vec;
847 virtual Belos::MultiVec<value_type> *
858 virtual Belos::MultiVec<value_type> *
859 CloneCopy(
const std::vector<int> &index)
const
861 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
863 new_multi_vec->vectors.resize(
index.size());
865 for (
unsigned int i = 0; i <
index.size(); ++i)
868 this->vectors.size(),
871 new_multi_vec->vectors[i] = std::make_shared<VectorType>();
874 *new_multi_vec->vectors[i] = *this->vectors[
index[i]];
877 return new_multi_vec;
885 virtual Belos::MultiVec<value_type> *
886 CloneViewNonConst(
const std::vector<int> &index)
888 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
890 new_multi_vec->vectors.resize(
index.size());
892 for (
unsigned int i = 0; i <
index.size(); ++i)
895 this->vectors.size(),
898 new_multi_vec->vectors[i].reset(
899 this->vectors[index[i]].get(),
905 return new_multi_vec;
913 virtual const Belos::MultiVec<value_type> *
914 CloneView(
const std::vector<int> &index)
const
916 auto new_multi_vec =
new MultiVecWrapper<VectorType>;
918 new_multi_vec->vectors.resize(
index.size());
920 for (
unsigned int i = 0; i <
index.size(); ++i)
923 this->vectors.size(),
926 new_multi_vec->vectors[i].reset(
927 this->vectors[index[i]].get(),
933 return new_multi_vec;
940 GetGlobalLength()
const
944 for (
unsigned int i = 1; i < this->vectors.size(); ++i)
947 return this->vectors[0]->size();
954 GetNumberVecs()
const
956 return vectors.size();
963 MvTimesMatAddMv(
const value_type alpha,
964 const Belos::MultiVec<value_type> &A_,
965 const Teuchos::SerialDenseMatrix<int, value_type> &B,
966 const value_type beta)
968 const auto &A = try_to_get_underlying_vector(A_);
970 const unsigned int n_rows = B.numRows();
971 const unsigned int n_cols = B.numCols();
973 AssertThrow(n_rows ==
static_cast<unsigned int>(A.GetNumberVecs()),
975 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
978 for (
unsigned int i = 0; i < n_cols; ++i)
979 (*this->vectors[i]) *= beta;
981 for (
unsigned int i = 0; i < n_cols; ++i)
982 for (
unsigned int j = 0; j < n_rows; ++j)
983 this->vectors[i]->add(alpha * B(j, i), *A.vectors[j]);
990 MvAddMv(
const value_type alpha,
991 const Belos::MultiVec<value_type> &A_,
992 const value_type beta,
993 const Belos::MultiVec<value_type> &B_)
995 const auto &A = try_to_get_underlying_vector(A_);
996 const auto &B = try_to_get_underlying_vector(B_);
998 AssertThrow(this->vectors.size() == A.vectors.size(),
1000 AssertThrow(this->vectors.size() == B.vectors.size(),
1003 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1005 this->vectors[i]->equ(alpha, *A.vectors[i]);
1006 this->vectors[i]->add(beta, *B.vectors[i]);
1014 MvScale(
const value_type alpha)
1016 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1017 (*this->vectors[i]) *= alpha;
1024 MvScale(
const std::vector<value_type> &alpha)
1035 MvTransMv(
const value_type alpha,
1036 const Belos::MultiVec<value_type> &A_,
1037 Teuchos::SerialDenseMatrix<int, value_type> &B)
const
1039 const auto &A = try_to_get_underlying_vector(A_);
1041 const unsigned int n_rows = B.numRows();
1042 const unsigned int n_cols = B.numCols();
1044 AssertThrow(n_rows ==
static_cast<unsigned int>(A.GetNumberVecs()),
1046 AssertThrow(n_cols ==
static_cast<unsigned int>(this->GetNumberVecs()),
1049 for (
unsigned int i = 0; i < n_rows; ++i)
1050 for (
unsigned int j = 0; j < n_cols; ++j)
1051 B(i, j) = alpha * ((*A.vectors[i]) * (*this->vectors[j]));
1059 MvDot(
const Belos::MultiVec<value_type> &A_,
1060 std::vector<value_type> &b)
const
1062 const auto &A = try_to_get_underlying_vector(A_);
1064 AssertThrow(this->vectors.size() == A.vectors.size(),
1068 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1069 b[i] = (*this->vectors[i]) * (*A.vectors[i]);
1077 std::vector<
typename Teuchos::ScalarTraits<value_type>::magnitudeType>
1079 Belos::NormType type = Belos::TwoNorm)
const
1084 for (
unsigned int i = 0; i < this->vectors.size(); ++i)
1085 normvec[i] = this->vectors[i]->
l2_norm();
1092 SetBlock(
const Belos::MultiVec<value_type> &A,
1093 const std::vector<int> &index)
1113 MvInit(
const value_type alpha)
1123 MvPrint(std::ostream &os)
const
1144 genericVector()
const
1154 static MultiVecWrapper<VectorType> &
1155 try_to_get_underlying_vector(Belos::MultiVec<value_type> &vec_in)
1157 auto vec =
dynamic_cast<MultiVecWrapper<VectorType> *
>(&vec_in);
1167 const static MultiVecWrapper<VectorType> &
1168 try_to_get_underlying_vector(
const Belos::MultiVec<value_type> &vec_in)
1170 auto vec =
dynamic_cast<const MultiVecWrapper<VectorType> *
>(&vec_in);
1178# ifdef HAVE_BELOS_TSQR
1180 factorExplicit(Belos::MultiVec<value_type> &,
1181 Teuchos::SerialDenseMatrix<int, value_type> &,
1189 Teuchos::SerialDenseMatrix<int, value_type> &,
1190 const typename Teuchos::ScalarTraits<value_type>::magnitudeType &)
1201 std::vector<std::shared_ptr<VectorType>> vectors;
1207 MultiVecWrapper() =
default;
1216 template <
typename OperatorType,
typename VectorType>
1217 class OperatorWrapper
1218 :
public Belos::Operator<typename VectorType::value_type>
1224 using value_type =
typename VectorType::value_type;
1230 this_type_is_missing_a_specialization()
1238 OperatorWrapper(
const OperatorType &op)
1245 virtual ~OperatorWrapper() =
default;
1251 Apply(
const Belos::MultiVec<value_type> &x,
1252 Belos::MultiVec<value_type> &y,
1253 Belos::ETrans trans = Belos::NOTRANS)
const
1258 op.vmult(MultiVecWrapper<VectorType>::try_to_get_underlying_vector(y)
1260 MultiVecWrapper<VectorType>::try_to_get_underlying_vector(x)
1268 HasApplyTranspose()
const
1278 const OperatorType &op;
1285 template <
typename VectorType>
1286 SolverBelos<VectorType>::SolverBelos(
1288 const AdditionalData &additional_data,
1289 const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters)
1290 : solver_control(solver_control)
1291 , additional_data(additional_data)
1292 , belos_parameters(belos_parameters)
1297 template <
typename VectorType>
1298 template <
typename OperatorType,
typename PreconditionerType>
1300 SolverBelos<VectorType>::solve(
const OperatorType &A_dealii,
1301 VectorType &x_dealii,
1302 const VectorType &b_dealii,
1303 const PreconditionerType &P_dealii)
1305 using value_type =
typename VectorType::value_type;
1307 using MV = Belos::MultiVec<value_type>;
1308 using OP = Belos::Operator<value_type>;
1310 Teuchos::RCP<OP> A = Teuchos::rcp(
1311 new internal::OperatorWrapper<OperatorType, VectorType>(A_dealii));
1312 Teuchos::RCP<OP> P = Teuchos::rcp(
1313 new internal::OperatorWrapper<PreconditionerType, VectorType>(P_dealii));
1314 Teuchos::RCP<MV> X =
1315 Teuchos::rcp(
new internal::MultiVecWrapper<VectorType>(x_dealii));
1316 Teuchos::RCP<MV> B =
1317 Teuchos::rcp(
new internal::MultiVecWrapper<VectorType>(b_dealii));
1319 Teuchos::RCP<Belos::LinearProblem<value_type, MV, OP>> problem =
1320 Teuchos::rcp(
new Belos::LinearProblem<value_type, MV, OP>(A, X, B));
1322 if (additional_data.right_preconditioning ==
false)
1323 problem->setLeftPrec(P);
1325 problem->setRightPrec(P);
1327 const bool problem_set = problem->setProblem();
1332 r.reinit(x_dealii,
true);
1333 A_dealii.vmult(r, x_dealii);
1334 r.sadd(-1., 1., b_dealii);
1335 const auto norm_0 = r.l2_norm();
1340 double relative_tolerance_to_be_achieved =
1341 solver_control.tolerance() / norm_0;
1342 const unsigned int max_steps = solver_control.max_steps();
1344 if (
const auto *reduction_control =
1346 relative_tolerance_to_be_achieved =
1347 std::max(relative_tolerance_to_be_achieved,
1348 reduction_control->reduction());
1350 Teuchos::RCP<Teuchos::ParameterList> belos_parameters_copy(
1351 Teuchos::rcp(
new Teuchos::ParameterList(*belos_parameters)));
1353 belos_parameters_copy->set(
"Convergence Tolerance",
1354 relative_tolerance_to_be_achieved);
1355 belos_parameters_copy->set(
"Maximum Iterations",
1356 static_cast<int>(max_steps));
1358 Teuchos::RCP<Belos::SolverManager<value_type, MV, OP>> solver;
1360 if (additional_data.solver_name == SolverName::cg)
1361 solver = Teuchos::rcp(
1362 new Belos::BlockCGSolMgr<value_type, MV, OP>(problem,
1363 belos_parameters_copy));
1364 else if (additional_data.solver_name == SolverName::gmres)
1365 solver = Teuchos::rcp(
1366 new Belos::BlockGmresSolMgr<value_type, MV, OP>(problem,
1367 belos_parameters_copy));
1371 const auto flag = solver->solve();
1373 solver_control.check(solver->getNumIters(), solver->achievedTol() * norm_0);
1375 AssertThrow(flag == Belos::ReturnType::Converged ||
1378 (solver_control.last_step() == max_steps)),
1380 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_NAMESPACE_CLOSE
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)
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
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