Reference documentation for deal.II version GIT a189bc2bdf 2022-12-07 02:45:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
trilinos_solver.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_trilinos_solver_h
17 #define dealii_trilinos_solver_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #ifdef DEAL_II_WITH_TRILINOS
23 
24 # include <deal.II/lac/exceptions.h>
27 # include <deal.II/lac/vector.h>
28 
29 // for AztecOO solvers
30 # include <Amesos.h>
31 # include <AztecOO.h>
32 # include <Epetra_LinearProblem.h>
33 # include <Epetra_Operator.h>
34 
35 // for Belos solvers
36 # ifdef DEAL_II_TRILINOS_WITH_BELOS
37 # include <BelosBlockCGSolMgr.hpp>
38 # include <BelosBlockGmresSolMgr.hpp>
39 # include <BelosEpetraAdapter.hpp>
40 # include <BelosIteration.hpp>
41 # include <BelosMultiVec.hpp>
42 # include <BelosOperator.hpp>
43 # include <BelosSolverManager.hpp>
44 # endif
45 
46 # include <memory>
47 
48 
50 
51 namespace TrilinosWrappers
52 {
53  // forward declarations
54 # ifndef DOXYGEN
55  class SparseMatrix;
56  class PreconditionBase;
57 # endif
58 
59 
78  class SolverBase
79  {
80  public:
89  {
93  cg,
97  cgs,
109  tfqmr
111 
117  {
126  explicit AdditionalData(const bool output_solver_details = false,
127  const unsigned int gmres_restart_parameter = 30);
128 
134 
138  const unsigned int gmres_restart_parameter;
139  };
140 
145  const AdditionalData &data = AdditionalData());
146 
151  SolverBase(const enum SolverName solver_name,
152  SolverControl & cn,
153  const AdditionalData &data = AdditionalData());
154 
158  virtual ~SolverBase() = default;
159 
165  void
166  solve(const SparseMatrix & A,
167  MPI::Vector & x,
168  const MPI::Vector & b,
169  const PreconditionBase &preconditioner);
170 
178  void
179  solve(const Epetra_Operator & A,
180  MPI::Vector & x,
181  const MPI::Vector & b,
182  const PreconditionBase &preconditioner);
183 
193  void
194  solve(const Epetra_Operator &A,
195  MPI::Vector & x,
196  const MPI::Vector & b,
197  const Epetra_Operator &preconditioner);
198 
208  void
209  solve(const Epetra_Operator & A,
210  Epetra_MultiVector & x,
211  const Epetra_MultiVector &b,
212  const PreconditionBase & preconditioner);
213 
224  void
225  solve(const Epetra_Operator & A,
226  Epetra_MultiVector & x,
227  const Epetra_MultiVector &b,
228  const Epetra_Operator & preconditioner);
229 
230 
231 
242  void
243  solve(const SparseMatrix & A,
244  ::Vector<double> & x,
245  const ::Vector<double> &b,
246  const PreconditionBase & preconditioner);
247 
259  void
261  ::Vector<double> & x,
262  const ::Vector<double> &b,
263  const PreconditionBase & preconditioner);
264 
271  void
272  solve(const SparseMatrix & A,
274  const ::LinearAlgebra::distributed::Vector<double> &b,
275  const PreconditionBase &preconditioner);
276 
284  void
287  const ::LinearAlgebra::distributed::Vector<double> &b,
288  const PreconditionBase &preconditioner);
289 
290 
294  SolverControl &
295  control() const;
296 
301  int,
302  << "An error with error number " << arg1
303  << " occurred while calling a Trilinos function");
304 
305  protected:
313 
314  private:
319  template <typename Preconditioner>
320  void
321  do_solve(const Preconditioner &preconditioner);
322 
326  template <typename Preconditioner>
327  void
328  set_preconditioner(AztecOO &solver, const Preconditioner &preconditioner);
329 
335  std::unique_ptr<Epetra_LinearProblem> linear_problem;
336 
341  std::unique_ptr<AztecOO_StatusTest> status_test;
342 
347  AztecOO solver;
348 
353  };
354 
355 
356  // provide a declaration for two explicit specializations
357  template <>
358  void
359  SolverBase::set_preconditioner(AztecOO & solver,
360  const PreconditionBase &preconditioner);
361 
362  template <>
363  void
364  SolverBase::set_preconditioner(AztecOO & solver,
365  const Epetra_Operator &preconditioner);
366 
367 
368 
374  class SolverCG : public SolverBase
375  {
376  public:
382  {
386  explicit AdditionalData(const bool output_solver_details = false);
387  };
388 
396  SolverCG(SolverControl &cn, const AdditionalData &data = AdditionalData());
397 
398  protected:
403  };
404 
405 
406 
412  class SolverCGS : public SolverBase
413  {
414  public:
419  {
423  explicit AdditionalData(const bool output_solver_details = false);
424  };
425 
434 
435  protected:
440  };
441 
442 
443 
448  class SolverGMRES : public SolverBase
449  {
450  public:
455  {
460  explicit AdditionalData(const bool output_solver_details = false,
461  const unsigned int restart_parameter = 30);
462  };
463 
472  const AdditionalData &data = AdditionalData());
473 
474  protected:
479  };
480 
481 
482 
489  class SolverBicgstab : public SolverBase
490  {
491  public:
496  {
500  explicit AdditionalData(const bool output_solver_details = false);
501  };
502 
511  const AdditionalData &data = AdditionalData());
512 
513  protected:
518  };
519 
520 
521 
528  class SolverTFQMR : public SolverBase
529  {
530  public:
535  {
539  explicit AdditionalData(const bool output_solver_details = false);
540  };
541 
550  const AdditionalData &data = AdditionalData());
551 
552  protected:
557  };
558 
559 
560 
574  {
575  public:
581  {
585  explicit AdditionalData(const bool output_solver_details = false,
586  const std::string &solver_type = "Amesos_Klu");
587 
593 
612  std::string solver_type;
613  };
614 
619  const AdditionalData &data = AdditionalData());
620 
624  virtual ~SolverDirect() = default;
625 
632  void
633  initialize(const SparseMatrix &A);
634 
640  void
641  solve(MPI::Vector &x, const MPI::Vector &b);
642 
648  void
650  const ::LinearAlgebra::distributed::Vector<double> &b);
651 
658  void
659  solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b);
660 
668  void
669  solve(const SparseMatrix & A,
670  ::Vector<double> & x,
671  const ::Vector<double> &b);
672 
679  void
680  solve(const SparseMatrix & A,
682  const ::LinearAlgebra::distributed::Vector<double> &b);
683 
687  SolverControl &
688  control() const;
689 
694  int,
695  << "An error with error number " << arg1
696  << " occurred while calling a Trilinos function");
697 
698  private:
703  void
704  do_solve();
705 
713 
719  std::unique_ptr<Epetra_LinearProblem> linear_problem;
720 
725  std::unique_ptr<Amesos_BaseSolver> solver;
726 
731  };
732 
733 
734 
735 # ifdef DEAL_II_TRILINOS_WITH_BELOS
742  template <typename VectorType>
743  class SolverBelos
744  {
745  public:
750  enum SolverName
751  {
755  cg,
759  gmres,
760  } solver_name;
761 
767  struct AdditionalData
768  {
772  AdditionalData(const SolverName solver_name = SolverName::cg,
773  const bool right_preconditioning = false)
774  : solver_name(solver_name)
775  , right_preconditioning(right_preconditioning)
776  {}
777 
781  SolverName solver_name;
782 
786  bool right_preconditioning;
787  };
788 
792  SolverBelos(SolverControl & solver_control,
793  const AdditionalData & additional_data,
794  const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters);
795 
799  template <typename OperatorType, typename PreconditionerType>
800  void
801  solve(const OperatorType & a,
802  VectorType & x,
803  const VectorType & b,
804  const PreconditionerType &p);
805 
806  private:
807  SolverControl & solver_control;
808  const AdditionalData additional_data;
809  const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters;
810  };
811 # endif
812 
813 } // namespace TrilinosWrappers
814 
815 
816 
817 # ifndef DOXYGEN
818 
819 # ifdef DEAL_II_TRILINOS_WITH_BELOS
820 namespace TrilinosWrappers
821 {
822  namespace internal
823  {
830  template <class VectorType>
831  class MultiVecWrapper
832  : public Belos::MultiVec<typename VectorType::value_type>
833  {
834  public:
838  using value_type = typename VectorType::value_type;
839 
843  static bool
844  this_type_is_missing_a_specialization()
845  {
846  return false;
847  }
848 
852  MultiVecWrapper(VectorType &vector)
853  {
854  this->vectors.resize(1);
855  this->vectors[0].reset(
856  &vector,
857  [](auto *) { /*Nothing to do, since vector is owned outside.*/ });
858  }
859 
863  MultiVecWrapper(const VectorType &vector)
864  {
865  this->vectors.resize(1);
866  this->vectors[0].reset(
867  &const_cast<VectorType &>(vector),
868  [](auto *) { /*Nothing to do, since vector is owned outside.*/ });
869  }
870 
874  virtual ~MultiVecWrapper() = default;
875 
879  virtual Belos::MultiVec<value_type> *
880  Clone(const int numvecs) const
881  {
882  auto new_multi_vec = new MultiVecWrapper<VectorType>;
883 
884  new_multi_vec->vectors.resize(numvecs);
885 
886  for (auto &vec : new_multi_vec->vectors)
887  {
888  vec = std::make_shared<VectorType>();
889 
890  AssertThrow(this->vectors.size() > 0, ExcInternalError());
891  vec->reinit(*this->vectors[0]);
892  }
893 
894  return new_multi_vec;
895  }
896 
900  virtual Belos::MultiVec<value_type> *
901  CloneCopy() const
902  {
903  AssertThrow(false, ExcNotImplemented());
904  }
905 
911  virtual Belos::MultiVec<value_type> *
912  CloneCopy(const std::vector<int> &index) const
913  {
914  auto new_multi_vec = new MultiVecWrapper<VectorType>;
915 
916  new_multi_vec->vectors.resize(index.size());
917 
918  for (unsigned int i = 0; i < index.size(); ++i)
919  {
920  AssertThrow(static_cast<unsigned int>(index[i]) <
921  this->vectors.size(),
922  ExcInternalError());
923 
924  new_multi_vec->vectors[i] = std::make_shared<VectorType>();
925 
926  AssertIndexRange(index[i], this->vectors.size());
927  *new_multi_vec->vectors[i] = *this->vectors[index[i]];
928  }
929 
930  return new_multi_vec;
931  }
932 
938  virtual Belos::MultiVec<value_type> *
939  CloneViewNonConst(const std::vector<int> &index)
940  {
941  auto new_multi_vec = new MultiVecWrapper<VectorType>;
942 
943  new_multi_vec->vectors.resize(index.size());
944 
945  for (unsigned int i = 0; i < index.size(); ++i)
946  {
947  AssertThrow(static_cast<unsigned int>(index[i]) <
948  this->vectors.size(),
949  ExcInternalError());
950 
951  new_multi_vec->vectors[i].reset(
952  this->vectors[index[i]].get(),
953  [](
954  auto
955  *) { /*Nothing to do, since we are creating only a view.*/ });
956  }
957 
958  return new_multi_vec;
959  }
960 
966  virtual const Belos::MultiVec<value_type> *
967  CloneView(const std::vector<int> &index) const
968  {
969  auto new_multi_vec = new MultiVecWrapper<VectorType>;
970 
971  new_multi_vec->vectors.resize(index.size());
972 
973  for (unsigned int i = 0; i < index.size(); ++i)
974  {
975  AssertThrow(static_cast<unsigned int>(index[i]) <
976  this->vectors.size(),
977  ExcInternalError());
978 
979  new_multi_vec->vectors[i].reset(
980  this->vectors[index[i]].get(),
981  [](
982  auto
983  *) { /*Nothing to do, since we are creating only a view.*/ });
984  }
985 
986  return new_multi_vec;
987  }
988 
992  virtual ptrdiff_t
993  GetGlobalLength() const
994  {
995  AssertThrow(this->vectors.size() > 0, ExcInternalError());
996 
997  for (unsigned int i = 1; i < this->vectors.size(); ++i)
998  AssertDimension(this->vectors[0]->size(), this->vectors[i]->size());
999 
1000  return this->vectors[0]->size();
1001  }
1002 
1006  virtual int
1007  GetNumberVecs() const
1008  {
1009  return vectors.size();
1010  }
1011 
1015  virtual void
1016  MvTimesMatAddMv(const value_type alpha,
1017  const Belos::MultiVec<value_type> & A_,
1018  const Teuchos::SerialDenseMatrix<int, value_type> &B,
1019  const value_type beta)
1020  {
1021  const auto &A = try_to_get_underlying_vector(A_);
1022 
1023  const unsigned int n_rows = B.numRows();
1024  const unsigned int n_cols = B.numCols();
1025 
1026  AssertThrow(n_rows == static_cast<unsigned int>(A.GetNumberVecs()),
1027  ExcInternalError());
1028  AssertThrow(n_cols == static_cast<unsigned int>(this->GetNumberVecs()),
1029  ExcInternalError());
1030 
1031  for (unsigned int i = 0; i < n_cols; ++i)
1032  (*this->vectors[i]) *= beta;
1033 
1034  for (unsigned int i = 0; i < n_cols; ++i)
1035  for (unsigned int j = 0; j < n_rows; ++j)
1036  this->vectors[i]->add(alpha * B(j, i), *A.vectors[j]);
1037  }
1038 
1042  virtual void
1043  MvAddMv(const value_type alpha,
1044  const Belos::MultiVec<value_type> &A_,
1045  const value_type beta,
1046  const Belos::MultiVec<value_type> &B_)
1047  {
1048  const auto &A = try_to_get_underlying_vector(A_);
1049  const auto &B = try_to_get_underlying_vector(B_);
1050 
1051  AssertThrow(this->vectors.size() == A.vectors.size(),
1052  ExcInternalError());
1053  AssertThrow(this->vectors.size() == B.vectors.size(),
1054  ExcInternalError());
1055 
1056  for (unsigned int i = 0; i < this->vectors.size(); ++i)
1057  {
1058  this->vectors[i]->equ(alpha, *A.vectors[i]);
1059  this->vectors[i]->add(beta, *B.vectors[i]);
1060  }
1061  }
1062 
1066  virtual void
1067  MvScale(const value_type alpha)
1068  {
1069  for (unsigned int i = 0; i < this->vectors.size(); ++i)
1070  (*this->vectors[i]) *= alpha;
1071  }
1072 
1076  virtual void
1077  MvScale(const std::vector<value_type> &alpha)
1078  {
1079  AssertThrow(false, ExcNotImplemented());
1080  (void)alpha;
1081  }
1082 
1087  virtual void
1088  MvTransMv(const value_type alpha,
1089  const Belos::MultiVec<value_type> & A_,
1090  Teuchos::SerialDenseMatrix<int, value_type> &B) const
1091  {
1092  const auto &A = try_to_get_underlying_vector(A_);
1093 
1094  const unsigned int n_rows = B.numRows();
1095  const unsigned int n_cols = B.numCols();
1096 
1097  AssertThrow(n_rows == static_cast<unsigned int>(A.GetNumberVecs()),
1098  ExcInternalError());
1099  AssertThrow(n_cols == static_cast<unsigned int>(this->GetNumberVecs()),
1100  ExcInternalError());
1101 
1102  for (unsigned int i = 0; i < n_rows; ++i)
1103  for (unsigned int j = 0; j < n_cols; ++j)
1104  B(i, j) = alpha * ((*A.vectors[i]) * (*this->vectors[j]));
1105  }
1106 
1111  virtual void
1112  MvDot(const Belos::MultiVec<value_type> &A_,
1113  std::vector<value_type> & b) const
1114  {
1115  const auto &A = try_to_get_underlying_vector(A_);
1116 
1117  AssertThrow(this->vectors.size() == A.vectors.size(),
1118  ExcInternalError());
1119  AssertThrow(this->vectors.size() == b.size(), ExcInternalError());
1120 
1121  for (unsigned int i = 0; i < this->vectors.size(); ++i)
1122  b[i] = (*this->vectors[i]) * (*A.vectors[i]);
1123  }
1124 
1128  virtual void
1129  MvNorm(
1130  std::vector<typename Teuchos::ScalarTraits<value_type>::magnitudeType>
1131  & normvec,
1132  Belos::NormType type = Belos::TwoNorm) const
1133  {
1134  AssertThrow(type == Belos::TwoNorm, ExcNotImplemented());
1135  AssertThrow(this->vectors.size() == normvec.size(), ExcInternalError());
1136 
1137  for (unsigned int i = 0; i < this->vectors.size(); ++i)
1138  normvec[i] = this->vectors[i]->l2_norm();
1139  }
1140 
1144  virtual void
1145  SetBlock(const Belos::MultiVec<value_type> &A,
1146  const std::vector<int> & index)
1147  {
1148  AssertThrow(false, ExcNotImplemented());
1149  (void)A;
1150  (void)index;
1151  }
1152 
1156  virtual void
1157  MvRandom()
1158  {
1159  AssertThrow(false, ExcNotImplemented());
1160  }
1161 
1165  virtual void
1166  MvInit(const value_type alpha)
1167  {
1168  AssertThrow(false, ExcNotImplemented());
1169  (void)alpha;
1170  }
1171 
1175  virtual void
1176  MvPrint(std::ostream &os) const
1177  {
1178  AssertThrow(false, ExcNotImplemented());
1179  (void)os;
1180  }
1181 
1185  VectorType &
1186  genericVector()
1187  {
1188  AssertThrow(GetNumberVecs() == 1, ExcNotImplemented());
1189 
1190  return *vectors[0];
1191  }
1192 
1196  const VectorType &
1197  genericVector() const
1198  {
1199  AssertThrow(GetNumberVecs() == 1, ExcNotImplemented());
1200 
1201  return *vectors[0];
1202  }
1203 
1207  static MultiVecWrapper<VectorType> &
1208  try_to_get_underlying_vector(Belos::MultiVec<value_type> &vec_in)
1209  {
1210  auto vec = dynamic_cast<MultiVecWrapper<VectorType> *>(&vec_in);
1211 
1212  AssertThrow(vec, ExcInternalError());
1213 
1214  return *vec;
1215  }
1216 
1220  const static MultiVecWrapper<VectorType> &
1221  try_to_get_underlying_vector(const Belos::MultiVec<value_type> &vec_in)
1222  {
1223  auto vec = dynamic_cast<const MultiVecWrapper<VectorType> *>(&vec_in);
1224 
1225  AssertThrow(vec, ExcInternalError());
1226 
1227  return *vec;
1228  }
1229 
1230 
1231 # ifdef HAVE_BELOS_TSQR
1232  virtual void
1233  factorExplicit(Belos::MultiVec<value_type> &,
1234  Teuchos::SerialDenseMatrix<int, value_type> &,
1235  const bool = false)
1236  {
1237  Assert(false, ExcNotImplemented());
1238  }
1239 
1240  virtual int
1241  revealRank(
1242  Teuchos::SerialDenseMatrix<int, value_type> &,
1243  const typename Teuchos::ScalarTraits<value_type>::magnitudeType &)
1244  {
1245  Assert(false, ExcNotImplemented());
1246  }
1247 
1248 # endif
1249 
1250  private:
1254  std::vector<std::shared_ptr<VectorType>> vectors;
1255 
1260  MultiVecWrapper() = default;
1261  };
1262 
1269  template <class OperatorType, class VectorType>
1270  class OperatorWrapper
1271  : public Belos::Operator<typename VectorType::value_type>
1272  {
1273  public:
1277  using value_type = typename VectorType::value_type;
1278 
1282  static bool
1283  this_type_is_missing_a_specialization()
1284  {
1285  return false;
1286  }
1287 
1291  OperatorWrapper(const OperatorType &op)
1292  : op(op)
1293  {}
1294 
1298  virtual ~OperatorWrapper() = default;
1299 
1303  virtual void
1304  Apply(const Belos::MultiVec<value_type> &x,
1305  Belos::MultiVec<value_type> & y,
1306  Belos::ETrans trans = Belos::NOTRANS) const
1307  {
1308  // TODO: check for Tvmult
1309  AssertThrow(trans == Belos::NOTRANS, ExcNotImplemented());
1310 
1311  op.vmult(MultiVecWrapper<VectorType>::try_to_get_underlying_vector(y)
1312  .genericVector(),
1313  MultiVecWrapper<VectorType>::try_to_get_underlying_vector(x)
1314  .genericVector());
1315  }
1316 
1320  virtual bool
1321  HasApplyTranspose() const
1322  {
1323  // TODO: check for Tvmult
1324  return false;
1325  }
1326 
1327  private:
1331  const OperatorType &op;
1332  };
1333 
1334  } // namespace internal
1335 
1336 
1337 
1338  template <typename VectorType>
1339  SolverBelos<VectorType>::SolverBelos(
1340  SolverControl & solver_control,
1341  const AdditionalData & additional_data,
1342  const Teuchos::RCP<Teuchos::ParameterList> &belos_parameters)
1343  : solver_control(solver_control)
1344  , additional_data(additional_data)
1345  , belos_parameters(belos_parameters)
1346  {}
1347 
1348 
1349 
1350  template <typename VectorType>
1351  template <typename OperatorType, typename PreconditionerType>
1352  void
1353  SolverBelos<VectorType>::solve(const OperatorType & A_dealii,
1354  VectorType & x_dealii,
1355  const VectorType & b_dealii,
1356  const PreconditionerType &P_dealii)
1357  {
1358  using value_type = typename VectorType::value_type;
1359 
1360  using MV = Belos::MultiVec<value_type>;
1361  using OP = Belos::Operator<value_type>;
1362 
1363  Teuchos::RCP<OP> A = Teuchos::rcp(
1364  new internal::OperatorWrapper<OperatorType, VectorType>(A_dealii));
1365  Teuchos::RCP<OP> P = Teuchos::rcp(
1366  new internal::OperatorWrapper<PreconditionerType, VectorType>(P_dealii));
1367  Teuchos::RCP<MV> X =
1368  Teuchos::rcp(new internal::MultiVecWrapper<VectorType>(x_dealii));
1369  Teuchos::RCP<MV> B =
1370  Teuchos::rcp(new internal::MultiVecWrapper<VectorType>(b_dealii));
1371 
1372  Teuchos::RCP<Belos::LinearProblem<value_type, MV, OP>> problem =
1373  Teuchos::rcp(new Belos::LinearProblem<value_type, MV, OP>(A, X, B));
1374 
1375  if (additional_data.right_preconditioning == false)
1376  problem->setLeftPrec(P);
1377  else
1378  problem->setRightPrec(P);
1379 
1380  const bool problem_set = problem->setProblem();
1381  AssertThrow(problem_set, ExcInternalError());
1382 
1383  // compute initial residal
1384  VectorType r;
1385  r.reinit(x_dealii, true);
1386  A_dealii.vmult(r, x_dealii);
1387  r.sadd(-1., 1., b_dealii);
1388  const auto norm_0 = r.l2_norm();
1389 
1390  if (solver_control.check(0, norm_0) != SolverControl::iterate)
1391  return;
1392 
1393  double relative_tolerance_to_be_achieved =
1394  solver_control.tolerance() / norm_0;
1395  const unsigned int max_steps = solver_control.max_steps();
1396 
1397  if (const auto *reduction_control =
1398  dynamic_cast<ReductionControl *>(&solver_control))
1399  relative_tolerance_to_be_achieved =
1400  std::max(relative_tolerance_to_be_achieved,
1401  reduction_control->reduction());
1402 
1403  Teuchos::RCP<Teuchos::ParameterList> belos_parameters_copy(
1404  Teuchos::rcp(new Teuchos::ParameterList(*belos_parameters)));
1405 
1406  belos_parameters_copy->set("Convergence Tolerance",
1407  relative_tolerance_to_be_achieved);
1408  belos_parameters_copy->set("Maximum Iterations",
1409  static_cast<int>(max_steps));
1410 
1411  Teuchos::RCP<Belos::SolverManager<value_type, MV, OP>> solver;
1412 
1413  if (additional_data.solver_name == SolverName::cg)
1414  solver = Teuchos::rcp(
1415  new Belos::BlockCGSolMgr<value_type, MV, OP>(problem,
1416  belos_parameters_copy));
1417  else if (additional_data.solver_name == SolverName::gmres)
1418  solver = Teuchos::rcp(
1419  new Belos::BlockGmresSolMgr<value_type, MV, OP>(problem,
1420  belos_parameters_copy));
1421  else
1422  AssertThrow(false, ExcNotImplemented());
1423 
1424  const auto flag = solver->solve();
1425 
1426  solver_control.check(solver->getNumIters(), solver->achievedTol() * norm_0);
1427 
1428  AssertThrow(flag == Belos::ReturnType::Converged ||
1429  ((dynamic_cast<IterationNumberControl *>(&solver_control) !=
1430  nullptr) &&
1431  (solver_control.last_step() == max_steps)),
1432  SolverControl::NoConvergence(solver_control.last_step(),
1433  solver_control.last_value()));
1434  }
1435 
1436 } // namespace TrilinosWrappers
1437 # endif
1438 
1439 # endif
1440 
1442 
1443 #endif // DEAL_II_WITH_TRILINOS
1444 
1445 #endif
@ 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)
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
const AdditionalData additional_data
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
virtual ~SolverDirect()=default
SolverDirect(SolverControl &cn, const AdditionalData &data=AdditionalData())
void initialize(const SparseMatrix &A)
std::unique_ptr< Amesos_BaseSolver > solver
void solve(MPI::Vector &x, const MPI::Vector &b)
std::unique_ptr< Epetra_LinearProblem > linear_problem
SolverControl & control() const
const AdditionalData additional_data
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1695
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:509
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
std::vector< value_type > l2_norm(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
static const char A
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
AdditionalData(const bool output_solver_details=false)
AdditionalData(const bool output_solver_details=false)
AdditionalData(const bool output_solver_details=false)
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
AdditionalData(const bool output_solver_details=false, const unsigned int restart_parameter=30)
AdditionalData(const bool output_solver_details=false)