18 #ifdef DEAL_II_WITH_TRILINOS
26 # include <AztecOO_StatusTest.h>
27 # include <AztecOO_StatusTestCombo.h>
28 # include <AztecOO_StatusTestMaxIters.h>
29 # include <AztecOO_StatusTestResNorm.h>
30 # include <AztecOO_StatusType.h>
41 const bool output_solver_details,
42 const unsigned int gmres_restart_parameter)
43 : output_solver_details(output_solver_details)
44 , gmres_restart_parameter(gmres_restart_parameter)
60 : solver_name(solver_name)
62 , additional_data(data)
84 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()),
86 const_cast<Epetra_MultiVector *
>(&
b.trilinos_vector()));
106 const_cast<Epetra_MultiVector *
>(
107 &
b.trilinos_vector()));
127 const_cast<Epetra_MultiVector *
>(
128 &
b.trilinos_vector()));
139 Epetra_MultiVector &x,
140 const Epetra_MultiVector &
b,
148 const_cast<Epetra_MultiVector *
>(
160 Epetra_MultiVector &x,
161 const Epetra_MultiVector &
b,
169 const_cast<Epetra_MultiVector *
>(
180 const ::Vector<double> &
b,
187 Assert(
A.local_range().second ==
A.m(),
188 ExcMessage(
"Can only work in serial when using deal.II vectors."));
189 Assert(
A.trilinos_matrix().Filled(),
190 ExcMessage(
"Matrix is not compressed. Call compress() method."));
192 Epetra_Vector ep_x(
View,
A.trilinos_matrix().DomainMap(), x.
begin());
193 Epetra_Vector ep_b(
View,
194 A.trilinos_matrix().RangeMap(),
195 const_cast<double *
>(
b.begin()));
200 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()), &ep_x, &ep_b);
210 const ::Vector<double> &
b,
213 Epetra_Vector ep_x(
View,
A.OperatorDomainMap(), x.
begin());
214 Epetra_Vector ep_b(
View,
215 A.OperatorRangeMap(),
216 const_cast<double *
>(
b.begin()));
220 linear_problem = std::make_unique<Epetra_LinearProblem>(&
A, &ep_x, &ep_b);
230 const ::LinearAlgebra::distributed::Vector<double> &
b,
236 A.trilinos_matrix().DomainMap().NumMyElements());
238 A.trilinos_matrix().RangeMap().NumMyElements());
240 Epetra_Vector ep_x(
View,
A.trilinos_matrix().DomainMap(), x.
begin());
241 Epetra_Vector ep_b(
View,
242 A.trilinos_matrix().RangeMap(),
243 const_cast<double *
>(
b.begin()));
248 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()), &ep_x, &ep_b);
258 const ::LinearAlgebra::distributed::Vector<double> &
b,
262 A.OperatorDomainMap().NumMyElements());
264 A.OperatorRangeMap().NumMyElements());
266 Epetra_Vector ep_x(
View,
A.OperatorDomainMap(), x.
begin());
267 Epetra_Vector ep_b(
View,
268 A.OperatorRangeMap(),
269 const_cast<double *
>(
b.begin()));
273 linear_problem = std::make_unique<Epetra_LinearProblem>(&
A, &ep_x, &ep_b);
284 compute_residual(
const Epetra_MultiVector *
const residual_vector)
286 Assert(residual_vector->NumVectors() == 1,
287 ExcMessage(
"Residual multivector holds more than one vector"));
289 const int ierr = residual_vector->Norm2(&res_l2_norm);
294 class TrilinosReductionControl :
public AztecOO_StatusTest
297 TrilinosReductionControl(
const int max_steps,
298 const double tolerance,
300 const Epetra_LinearProblem &linear_problem);
302 virtual ~TrilinosReductionControl()
override =
default;
305 ResidualVectorRequired()
const override
310 virtual AztecOO_StatusType
311 CheckStatus(
int CurrentIter,
312 Epetra_MultiVector *CurrentResVector,
313 double CurrentResNormEst,
314 bool SolutionUpdated)
override
319 (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
321 if (CurrentIter == 0)
330 virtual AztecOO_StatusType
331 GetStatus()
const override
336 virtual std::ostream &
337 Print(std::ostream &stream,
int indent = 0)
const override
343 get_initial_residual()
const
349 get_current_residual()
const
364 TrilinosReductionControl::TrilinosReductionControl(
366 const double tolerance,
368 const Epetra_LinearProblem &linear_problem)
374 AztecOO_StatusTestCombo::OR))
379 std::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
382 Assert(linear_problem.GetRHS()->NumVectors() == 1,
383 ExcMessage(
"RHS multivector holds more than one vector"));
387 *linear_problem.GetOperator(),
388 *(linear_problem.GetLHS()->operator()(0)),
389 *(linear_problem.GetRHS()->operator()(0)),
392 AztecOO_StatusTestResNorm::TwoNorm);
394 AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
399 *linear_problem.GetOperator(),
400 *(linear_problem.GetLHS()->operator()(0)),
401 *(linear_problem.GetRHS()->operator()(0)),
404 AztecOO_StatusTestResNorm::TwoNorm);
406 AztecOO_StatusTestResNorm::NormOfInitRes,
407 AztecOO_StatusTestResNorm::TwoNorm);
415 template <
typename Preconditioner>
417 SolverBase::do_solve(
const Preconditioner &preconditioner)
422 solver.SetProblem(*linear_problem);
428 solver.SetAztecOption(AZ_solver, AZ_cg);
431 solver.SetAztecOption(AZ_solver, AZ_cgs);
434 solver.SetAztecOption(AZ_solver, AZ_gmres);
435 solver.SetAztecOption(AZ_kspace,
436 additional_data.gmres_restart_parameter);
439 solver.SetAztecOption(AZ_solver, AZ_bicgstab);
442 solver.SetAztecOption(AZ_solver, AZ_tfqmr);
449 set_preconditioner(solver, preconditioner);
452 solver.SetAztecOption(AZ_output,
453 additional_data.output_solver_details ? AZ_all :
455 solver.SetAztecOption(AZ_conv, AZ_noscaled);
475 status_test = std::make_unique<internal::TrilinosReductionControl>(
476 reduction_control->max_steps(),
477 reduction_control->tolerance(),
478 reduction_control->reduction(),
480 solver.SetStatusTest(status_test.get());
486 solver.Iterate(solver_control.max_steps(), solver_control.tolerance());
496 "option not implemented"));
501 "numerical breakdown"));
506 "loss of precision"));
511 "GMRES Hessenberg ill-conditioned"));
524 if (
const internal::TrilinosReductionControl
525 *
const reduction_control_status =
526 dynamic_cast<const internal::TrilinosReductionControl *const
>(
535 if (solver.NumIters() > 0)
540 solver_control.check(
541 0, reduction_control_status->get_initial_residual());
542 solver_control.check(
544 reduction_control_status->get_current_residual());
547 solver_control.check(
549 reduction_control_status->get_current_residual());
554 solver_control.check(solver.NumIters(), solver.TrueResidual());
560 solver_control.last_value()));
567 SolverBase::set_preconditioner(AztecOO &solver,
574 const int ierr = solver.SetPrecOperator(
579 solver.SetAztecOption(AZ_precond, AZ_none);
585 SolverBase::set_preconditioner(AztecOO &solver,
589 solver.SetPrecOperator(
const_cast<Epetra_Operator *
>(&preconditioner));
633 const std::string &solver_type)
634 : output_solver_details(output_solver_details)
635 , solver_type(solver_type)
664 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()));
681 std::string(
"You tried to select the solver type <") +
683 "> but this solver is not supported by Trilinos either "
684 "because it does not exist, or because Trilinos was not "
685 "configured for its use."));
690 verbose_cout <<
"Starting symbolic factorization" << std::endl;
691 ierr =
solver->SymbolicFactorization();
694 verbose_cout <<
"Starting numeric factorization" << std::endl;
695 ierr =
solver->NumericFactorization();
708 const_cast<Epetra_MultiVector *
>(&
b.trilinos_vector()));
716 verbose_cout <<
"Starting solve" << std::endl;
719 int ierr =
solver->Solve();
731 const ::LinearAlgebra::distributed::Vector<double> &
b)
733 Epetra_Vector ep_x(
View,
736 Epetra_Vector ep_b(
View,
738 const_cast<double *
>(
b.begin()));
751 verbose_cout <<
"Starting solve" << std::endl;
754 int ierr =
solver->Solve();
781 std::string(
"You tried to select the solver type <") +
783 "> but this solver is not supported by Trilinos either "
784 "because it does not exist, or because Trilinos was not "
785 "configured for its use."));
790 verbose_cout <<
"Starting symbolic factorization" << std::endl;
791 ierr =
solver->SymbolicFactorization();
794 verbose_cout <<
"Starting numeric factorization" << std::endl;
795 ierr =
solver->NumericFactorization();
798 verbose_cout <<
"Starting solve" << std::endl;
822 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()),
824 const_cast<Epetra_MultiVector *
>(&
b.trilinos_vector()));
834 const ::Vector<double> &
b)
840 Assert(
A.local_range().second ==
A.m(),
841 ExcMessage(
"Can only work in serial when using deal.II vectors."));
842 Epetra_Vector ep_x(
View,
A.trilinos_matrix().DomainMap(), x.
begin());
843 Epetra_Vector ep_b(
View,
844 A.trilinos_matrix().RangeMap(),
845 const_cast<double *
>(
b.begin()));
850 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()), &ep_x, &ep_b);
861 const ::LinearAlgebra::distributed::Vector<double> &
b)
864 A.trilinos_matrix().DomainMap().NumMyElements());
866 A.trilinos_matrix().RangeMap().NumMyElements());
867 Epetra_Vector ep_x(
View,
A.trilinos_matrix().DomainMap(), x.
begin());
868 Epetra_Vector ep_b(
View,
869 A.trilinos_matrix().RangeMap(),
870 const_cast<double *
>(
b.begin()));
875 const_cast<Epetra_CrsMatrix *
>(&
A.trilinos_matrix()), &ep_x, &ep_b);
size_type locally_owned_size() const
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
unsigned int last_step() const
double last_value() const
virtual State check(const unsigned int step, const double check_value)
@ 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
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())
const AdditionalData additional_data
SolverDirect(SolverControl &cn, const AdditionalData &data=AdditionalData())
void initialize(const SparseMatrix &A)
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
SolverControl & control() const
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
__global__ void reduction(Number *result, const Number *v, const size_type N)
static ::ExceptionBase & ExcTrilinosError(int arg1)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcTrilinosError(int arg1)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
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, const std::string &solver_type="Amesos_Klu")
bool output_solver_details
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_abs_tol
std::unique_ptr< AztecOO_StatusTestResNorm > status_test_rel_tol
std::unique_ptr< AztecOO_StatusTestMaxIters > status_test_max_steps
std::unique_ptr< AztecOO_StatusTestCombo > status_test_collection