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>
40 const bool output_solver_details,
41 const unsigned int gmres_restart_parameter)
42 : output_solver_details(output_solver_details)
43 , gmres_restart_parameter(gmres_restart_parameter)
59 : solver_name(solver_name)
61 , additional_data(
data)
83 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()),
85 const_cast<Epetra_MultiVector *
>(&b.trilinos_vector()));
103 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
104 &
x.trilinos_vector(),
105 const_cast<Epetra_MultiVector *
>(
106 &b.trilinos_vector()));
124 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
125 &
x.trilinos_vector(),
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,
186 Assert(A.local_range().second == A.m(),
187 ExcMessage(
"Can only work in serial when using deal.II vectors."));
188 Assert(A.trilinos_matrix().Filled(),
189 ExcMessage(
"Matrix is not compressed. Call compress() method."));
193 A.trilinos_matrix().RangeMap(),
194 const_cast<double *
>(b.begin()));
199 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()), &
ep_x, &
ep_b);
209 const ::Vector<double> &b,
214 A.OperatorRangeMap(),
215 const_cast<double *
>(b.begin()));
229 const ::LinearAlgebra::distributed::Vector<double> &b,
235 A.trilinos_matrix().DomainMap().NumMyElements());
237 A.trilinos_matrix().RangeMap().NumMyElements());
241 A.trilinos_matrix().RangeMap(),
242 const_cast<double *
>(b.begin()));
247 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()), &
ep_x, &
ep_b);
257 const ::LinearAlgebra::distributed::Vector<double> &b,
258 const PreconditionBase &preconditioner)
261 A.OperatorDomainMap().NumMyElements());
263 A.OperatorRangeMap().NumMyElements());
267 A.OperatorRangeMap(),
268 const_cast<double *
>(b.begin()));
286 ExcMessage(
"Residual multivector holds more than one vector"));
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;
335 virtual std::ostream &
336 Print(std::ostream &stream,
int indent = 0)
const override
363 TrilinosReductionControl::TrilinosReductionControl(
365 const double tolerance,
366 const double reduction,
367 const Epetra_LinearProblem &linear_problem)
378 std::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
381 Assert(linear_problem.GetRHS()->NumVectors() == 1,
382 ExcMessage(
"RHS multivector holds more than one vector"));
386 *linear_problem.GetOperator(),
387 *(linear_problem.GetLHS()->operator()(0)),
388 *(linear_problem.GetRHS()->operator()(0)),
398 *linear_problem.GetOperator(),
399 *(linear_problem.GetLHS()->operator()(0)),
400 *(linear_problem.GetRHS()->operator()(0)),
414 template <
typename Preconditioner>
421 solver.SetProblem(*linear_problem);
435 additional_data.gmres_restart_parameter);
448 set_preconditioner(solver, preconditioner);
452 additional_data.output_solver_details ?
AZ_all :
474 status_test = std::make_unique<internal::TrilinosReductionControl>(
479 solver.SetStatusTest(status_test.get());
485 solver.Iterate(solver_control.max_steps(), solver_control.tolerance());
495 "option not implemented"));
500 "numerical breakdown"));
505 "loss of precision"));
510 "GMRES Hessenberg ill-conditioned"));
523 if (
const internal::TrilinosReductionControl
525 dynamic_cast<const internal::TrilinosReductionControl *
>(
534 if (solver.NumIters() > 0)
539 solver_control.check(
541 solver_control.check(
546 solver_control.check(
553 solver_control.check(solver.NumIters(), solver.TrueResidual());
559 solver_control.last_value()));
566 SolverBase::set_preconditioner(AztecOO &solver,
573 const int ierr = solver.SetPrecOperator(
584 SolverBase::set_preconditioner(AztecOO &solver,
588 solver.SetPrecOperator(
const_cast<Epetra_Operator *
>(&preconditioner));
632 const std::string &solver_type)
633 : output_solver_details(output_solver_details)
634 , solver_type(solver_type)
648 , additional_data(
data.output_solver_details,
data.solver_type)
670 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()));
687 "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;
700 verbose_cout <<
"Starting numeric factorization" << std::endl;
727 const ::LinearAlgebra::distributed::Vector<double> &b)
741 const_cast<Epetra_MultiVector *
>(&b.trilinos_vector()));
764 const ::LinearAlgebra::distributed::Vector<double> &b)
const
771 const_cast<double *
>(b.begin()));
814 "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;
827 verbose_cout <<
"Starting numeric factorization" << std::endl;
849 Epetra_MultiVector &
x,
850 const Epetra_MultiVector &b)
855 std::make_unique<Epetra_LinearProblem>(
const_cast<Epetra_Operator *
>(&A),
857 const_cast<Epetra_MultiVector *
>(
874 const unsigned int m =
rhs.m();
875 const unsigned int n =
rhs.n();
886 for (
unsigned int i = 0; i < n; ++i)
895 mat.OperatorRangeMap(),
899 mat.OperatorDomainMap(),
918 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()),
919 &
x.trilinos_vector(),
920 const_cast<Epetra_MultiVector *
>(&b.trilinos_vector()));
930 const ::Vector<double> &b)
936 Assert(A.local_range().second == A.m(),
937 ExcMessage(
"Can only work in serial when using deal.II vectors."));
940 A.trilinos_matrix().RangeMap(),
941 const_cast<double *
>(b.begin()));
946 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()), &
ep_x, &
ep_b);
957 const ::LinearAlgebra::distributed::Vector<double> &b)
960 A.trilinos_matrix().DomainMap().NumMyElements());
962 A.trilinos_matrix().RangeMap().NumMyElements());
965 A.trilinos_matrix().RangeMap(),
966 const_cast<double *
>(b.begin()));
971 const_cast<Epetra_CrsMatrix *
>(&A.trilinos_matrix()), &
ep_x, &
ep_b);
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.
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())
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_NAMESPACE_CLOSE
#define DEAL_II_NOT_IMPLEMENTED()
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)
std::vector< index_type > data
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