17#ifdef DEAL_II_WITH_TRILINOS
25# include <AztecOO_StatusTest.h>
26# include <AztecOO_StatusTestCombo.h>
27# include <AztecOO_StatusTestMaxIters.h>
28# include <AztecOO_StatusTestResNorm.h>
29# include <AztecOO_StatusType.h>
85 const_cast<Epetra_MultiVector *
>(&b.trilinos_vector()));
103 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
105 const_cast<Epetra_MultiVector *
>(
106 &b.trilinos_vector()));
124 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
126 const_cast<Epetra_MultiVector *
>(
127 &b.trilinos_vector()));
138 Epetra_MultiVector &x,
139 const Epetra_MultiVector &b,
145 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
147 const_cast<Epetra_MultiVector *
>(
159 Epetra_MultiVector &x,
160 const Epetra_MultiVector &b,
166 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
168 const_cast<Epetra_MultiVector *
>(
179 const ::Vector<double> &b,
187 ExcMessage(
"Can only work in serial when using deal.II vectors."));
189 ExcMessage(
"Matrix is not compressed. Call compress() method."));
192 Epetra_Vector ep_b(View,
194 const_cast<double *
>(b.begin()));
209 const ::Vector<double> &b,
212 Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.
begin());
213 Epetra_Vector ep_b(View,
214 A.OperatorRangeMap(),
215 const_cast<double *
>(b.begin()));
219 linear_problem = std::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
229 const ::LinearAlgebra::distributed::Vector<double> &b,
240 Epetra_Vector ep_b(View,
242 const_cast<double *
>(b.begin()));
257 const ::LinearAlgebra::distributed::Vector<double> &b,
261 A.OperatorDomainMap().NumMyElements());
263 A.OperatorRangeMap().NumMyElements());
265 Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.
begin());
266 Epetra_Vector ep_b(View,
267 A.OperatorRangeMap(),
268 const_cast<double *
>(b.begin()));
272 linear_problem = std::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
283 compute_residual(
const Epetra_MultiVector *
const residual_vector)
285 Assert(residual_vector->NumVectors() == 1,
286 ExcMessage(
"Residual multivector holds more than one vector"));
288 const int ierr = residual_vector->Norm2(&res_l2_norm);
293 class TrilinosReductionControl :
public AztecOO_StatusTest
296 TrilinosReductionControl(
const int max_steps,
297 const double tolerance,
298 const double reduction,
299 const Epetra_LinearProblem &linear_problem);
301 virtual ~TrilinosReductionControl()
override =
default;
304 ResidualVectorRequired()
const override
306 return status_test_collection->ResidualVectorRequired();
309 virtual AztecOO_StatusType
310 CheckStatus(
int CurrentIter,
311 Epetra_MultiVector *CurrentResVector,
312 double CurrentResNormEst,
313 bool SolutionUpdated)
override
318 (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
320 if (CurrentIter == 0)
321 initial_residual = current_residual;
323 return status_test_collection->CheckStatus(CurrentIter,
329 virtual AztecOO_StatusType
330 GetStatus()
const override
332 return status_test_collection->GetStatus();
335 virtual std::ostream &
336 Print(std::ostream &stream,
int indent = 0)
const override
338 return status_test_collection->Print(stream, indent);
342 get_initial_residual()
const
344 return initial_residual;
348 get_current_residual()
const
350 return current_residual;
354 double initial_residual;
355 double current_residual;
356 std::unique_ptr<AztecOO_StatusTestCombo> status_test_collection;
357 std::unique_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
358 std::unique_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
359 std::unique_ptr<AztecOO_StatusTestResNorm> status_test_rel_tol;
363 TrilinosReductionControl::TrilinosReductionControl(
365 const double tolerance,
366 const double reduction,
367 const Epetra_LinearProblem &linear_problem)
368 : initial_residual(std::numeric_limits<double>::
max())
369 , current_residual(std::numeric_limits<double>::
max())
372 , status_test_collection(std::make_unique<AztecOO_StatusTestCombo>(
373 AztecOO_StatusTestCombo::OR))
377 status_test_max_steps =
378 std::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
379 status_test_collection->AddStatusTest(*status_test_max_steps);
382 ExcMessage(
"RHS multivector holds more than one vector"));
385 status_test_abs_tol = std::make_unique<AztecOO_StatusTestResNorm>(
390 status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
391 AztecOO_StatusTestResNorm::TwoNorm);
392 status_test_abs_tol->DefineScaleForm(
393 AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
394 status_test_collection->AddStatusTest(*status_test_abs_tol);
397 status_test_rel_tol = std::make_unique<AztecOO_StatusTestResNorm>(
402 status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
403 AztecOO_StatusTestResNorm::TwoNorm);
404 status_test_rel_tol->DefineScaleForm(
405 AztecOO_StatusTestResNorm::NormOfInitRes,
406 AztecOO_StatusTestResNorm::TwoNorm);
407 status_test_collection->AddStatusTest(*status_test_rel_tol);
414 template <
typename Preconditioner>
427 solver.SetAztecOption(AZ_solver, AZ_cg);
430 solver.SetAztecOption(AZ_solver, AZ_cgs);
433 solver.SetAztecOption(AZ_solver, AZ_gmres);
434 solver.SetAztecOption(AZ_kspace,
438 solver.SetAztecOption(AZ_solver, AZ_bicgstab);
441 solver.SetAztecOption(AZ_solver, AZ_tfqmr);
451 solver.SetAztecOption(AZ_output,
454 solver.SetAztecOption(AZ_conv, AZ_noscaled);
474 status_test = std::make_unique<internal::TrilinosReductionControl>(
475 reduction_control->max_steps(),
476 reduction_control->tolerance(),
477 reduction_control->reduction(),
495 "option not implemented"));
500 "numerical breakdown"));
505 "loss of precision"));
510 "GMRES Hessenberg ill-conditioned"));
523 if (
const internal::TrilinosReductionControl
524 *
const reduction_control_status =
525 dynamic_cast<const internal::TrilinosReductionControl *const
>(
534 if (
solver.NumIters() > 0)
540 0, reduction_control_status->get_initial_residual());
543 reduction_control_status->get_current_residual());
548 reduction_control_status->get_current_residual());
573 const int ierr =
solver.SetPrecOperator(
578 solver.SetAztecOption(AZ_precond, AZ_none);
687 std::string(
"You tried to select the solver type <") +
689 "> but this solver is not supported by Trilinos either "
690 "because it does not exist, or because Trilinos was not "
691 "configured for its use."));
696 verbose_cout <<
"Starting symbolic factorization" << std::endl;
697 ierr =
solver->SymbolicFactorization();
700 verbose_cout <<
"Starting numeric factorization" << std::endl;
701 ierr =
solver->NumericFactorization();
727 const ::LinearAlgebra::distributed::Vector<double> &b)
741 const_cast<Epetra_MultiVector *
>(&b.trilinos_vector()));
749 verbose_cout <<
"Starting solve" << std::endl;
752 int ierr =
solver->Solve();
764 const ::LinearAlgebra::distributed::Vector<double> &b)
const
766 Epetra_Vector ep_x(View,
769 Epetra_Vector ep_b(View,
771 const_cast<double *
>(b.begin()));
784 verbose_cout <<
"Starting solve" << std::endl;
787 int ierr =
solver->Solve();
814 std::string(
"You tried to select the solver type <") +
816 "> but this solver is not supported by Trilinos either "
817 "because it does not exist, or because Trilinos was not "
818 "configured for its use."));
823 verbose_cout <<
"Starting symbolic factorization" << std::endl;
824 ierr =
solver->SymbolicFactorization();
827 verbose_cout <<
"Starting numeric factorization" << std::endl;
828 ierr =
solver->NumericFactorization();
831 verbose_cout <<
"Starting solve" << std::endl;
857 const_cast<Epetra_MultiVector *
>(&b.trilinos_vector()));
867 const ::Vector<double> &b)
874 ExcMessage(
"Can only work in serial when using deal.II vectors."));
876 Epetra_Vector ep_b(View,
878 const_cast<double *
>(b.begin()));
894 const ::LinearAlgebra::distributed::Vector<double> &b)
901 Epetra_Vector ep_b(View,
903 const_cast<double *
>(b.begin()));
size_type locally_owned_size() const
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
@ success
Stop iteration, goal reached.
const Epetra_MultiVector & trilinos_vector() const
Teuchos::RCP< Epetra_Operator > preconditioner
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
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
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())
const Epetra_CrsMatrix & trilinos_matrix() const
std::pair< size_type, size_type > local_range() const
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
#define DEAL_II_NOT_IMPLEMENTED()
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 bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
bool output_solver_details