16 #include <deal.II/lac/trilinos_solver.h> 18 #ifdef DEAL_II_WITH_TRILINOS 20 # include <deal.II/base/conditional_ostream.h> 21 # include <deal.II/lac/trilinos_sparse_matrix.h> 22 # include <deal.II/lac/trilinos_vector_base.h> 23 # include <deal.II/lac/trilinos_precondition.h> 25 DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
26 # include <AztecOO_StatusTest.h> 27 # include <AztecOO_StatusTestMaxIters.h> 28 # include <AztecOO_StatusTestResNorm.h> 29 # include <AztecOO_StatusTestCombo.h> 30 # include <AztecOO_StatusType.h> 31 DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
36 DEAL_II_NAMESPACE_OPEN
42 const unsigned int gmres_restart_parameter)
44 output_solver_details (output_solver_details),
45 gmres_restart_parameter (gmres_restart_parameter)
61 solver_name (solver_name),
91 (
new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.
trilinos_matrix()),
113 (
new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
128 const Epetra_Operator &preconditioner)
135 (
new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
148 Epetra_MultiVector &x,
149 const Epetra_MultiVector &b,
157 (
new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
159 const_cast<Epetra_MultiVector *>(&b)));
170 Epetra_MultiVector &x,
171 const Epetra_MultiVector &b,
172 const Epetra_Operator &preconditioner)
179 (
new Epetra_LinearProblem(const_cast<Epetra_Operator *>(&A),
181 const_cast<Epetra_MultiVector *>(&b)));
191 const ::Vector<double> &b,
203 ExcMessage (
"Can only work in serial when using deal.II vectors."));
205 ExcMessage (
"Matrix is not compressed. Call compress() method."));
208 Epetra_Vector ep_b (View, A.
range_partitioner(),
const_cast<double *
>(b.begin()));
224 const ::Vector<double> &b,
229 Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.
begin());
230 Epetra_Vector ep_b (View, A.OperatorRangeMap(),
const_cast<double *
>(b.begin()));
244 const ::LinearAlgebra::distributed::Vector<double> &b,
253 AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
257 Epetra_Vector ep_b (View, A.
range_partitioner(),
const_cast<double *
>(b.begin()));
273 const ::LinearAlgebra::distributed::Vector<double> &b,
279 A.OperatorDomainMap().NumMyElements());
280 AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
281 A.OperatorRangeMap().NumMyElements());
283 Epetra_Vector ep_x (View, A.OperatorDomainMap(), x.
begin());
284 Epetra_Vector ep_b (View, A.OperatorRangeMap(),
const_cast<double *
>(b.begin()));
299 compute_residual (
const Epetra_MultiVector *
const residual_vector)
301 Assert(residual_vector->NumVectors() == 1,
302 ExcMessage(
"Residual multivector holds more than one vector"));
303 TrilinosScalar res_l2_norm = 0.0;
304 const int ierr = residual_vector->Norm2 (&res_l2_norm);
309 class TrilinosReductionControl :
public AztecOO_StatusTest
312 TrilinosReductionControl (
const double &max_steps,
313 const double &tolerance,
314 const double &reduction,
317 virtual ~TrilinosReductionControl() {}
320 ResidualVectorRequired ()
const 322 return status_test_collection->ResidualVectorRequired();
325 virtual AztecOO_StatusType
326 CheckStatus (
int CurrentIter,
327 Epetra_MultiVector *CurrentResVector,
328 double CurrentResNormEst,
329 bool SolutionUpdated)
333 current_residual = (CurrentResNormEst < 0.0 ?
334 compute_residual(CurrentResVector) :
336 if (CurrentIter == 0)
337 initial_residual = current_residual;
339 return status_test_collection->CheckStatus(CurrentIter,
346 virtual AztecOO_StatusType
349 return status_test_collection->GetStatus();
352 virtual std::ostream &
353 Print (std::ostream &stream,
354 int indent = 0)
const 356 return status_test_collection->Print(stream,indent);
360 get_initial_residual()
const 362 return initial_residual;
366 get_current_residual()
const 368 return current_residual;
372 double initial_residual;
373 double current_residual;
374 std_cxx11::shared_ptr<AztecOO_StatusTestCombo> status_test_collection;
375 std_cxx11::shared_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
376 std_cxx11::shared_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
377 std_cxx11::shared_ptr<AztecOO_StatusTestResNorm> status_test_rel_tol;
381 TrilinosReductionControl::TrilinosReductionControl(
382 const double &max_steps,
383 const double &tolerance,
384 const double &reduction,
387 initial_residual (std::numeric_limits<double>::max()),
388 current_residual(std::numeric_limits<double>::max())
392 status_test_collection.reset(
393 new AztecOO_StatusTestCombo (AztecOO_StatusTestCombo::OR) );
396 status_test_max_steps.reset(
397 new AztecOO_StatusTestMaxIters(max_steps) );
398 status_test_collection->AddStatusTest(*status_test_max_steps);
400 Assert(linear_problem.GetRHS()->NumVectors() == 1,
401 ExcMessage(
"RHS multivector holds more than one vector"));
404 status_test_abs_tol.reset(
405 new AztecOO_StatusTestResNorm(*linear_problem.GetOperator(),
406 *(linear_problem.GetLHS()->operator()(0)),
407 *(linear_problem.GetRHS()->operator()(0)),
409 status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
410 AztecOO_StatusTestResNorm::TwoNorm);
411 status_test_abs_tol->DefineScaleForm(AztecOO_StatusTestResNorm::None,
412 AztecOO_StatusTestResNorm::TwoNorm);
413 status_test_collection->AddStatusTest(*status_test_abs_tol);
416 status_test_rel_tol.reset(
417 new AztecOO_StatusTestResNorm(*linear_problem.GetOperator(),
418 *(linear_problem.GetLHS()->operator()(0)),
419 *(linear_problem.GetRHS()->operator()(0)),
421 status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
422 AztecOO_StatusTestResNorm::TwoNorm);
423 status_test_rel_tol->DefineScaleForm(AztecOO_StatusTestResNorm::NormOfInitRes,
424 AztecOO_StatusTestResNorm::TwoNorm);
425 status_test_collection->AddStatusTest(*status_test_rel_tol);
432 template<
typename Preconditioner>
445 solver.SetAztecOption(AZ_solver, AZ_cg);
448 solver.SetAztecOption(AZ_solver, AZ_cgs);
451 solver.SetAztecOption(AZ_solver, AZ_gmres);
455 solver.SetAztecOption(AZ_solver, AZ_bicgstab);
458 solver.SetAztecOption(AZ_solver, AZ_tfqmr);
470 solver.SetAztecOption (AZ_conv, AZ_noscaled);
490 status_test.reset(
new internal::TrilinosReductionControl(
491 reduction_control->max_steps(),
492 reduction_control->tolerance(),
493 reduction_control->reduction(),
510 "option not implemented"));
513 "numerical breakdown"));
516 "loss of precision"));
519 "GMRES Hessenberg ill-conditioned"));
531 if (
const internal::TrilinosReductionControl*
const reduction_control_status
532 = dynamic_cast<const internal::TrilinosReductionControl *const>(
status_test.get()))
539 if (
solver.NumIters() > 0)
572 const int ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>
577 solver.SetAztecOption(AZ_precond,AZ_none);
584 const Epetra_Operator &preconditioner)
586 const int ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>(&preconditioner));
596 output_solver_details (output_solver_details)
615 const unsigned int restart_parameter)
617 output_solver_details (output_solver_details),
618 restart_parameter (restart_parameter)
639 output_solver_details (output_solver_details)
660 output_solver_details (output_solver_details)
681 output_solver_details (output_solver_details)
701 const std::string &solver_type)
703 output_solver_details (output_solver_details),
704 solver_type(solver_type)
758 ExcMessage (std::string (
"You tried to select the solver type <") +
760 "> but this solver is not supported by Trilinos either " 761 "because it does not exist, or because Trilinos was not " 762 "configured for its use.")
769 verbose_cout <<
"Starting symbolic factorization" << std::endl;
770 ierr = solver->SymbolicFactorization();
773 verbose_cout <<
"Starting numeric factorization" << std::endl;
774 ierr = solver->NumericFactorization();
796 verbose_cout <<
"Starting solve" << std::endl;
797 ierr = solver->Solve ();
826 string (
"You tried to select the solver type <")
828 "> but this solver is not supported by Trilinos either " 829 "because it does not exist, or because Trilinos was not " 830 "configured for its use."));
832 solver.reset (Factory.
836 verbose_cout <<
"Starting symbolic factorization" << std::endl;
837 ierr = solver->SymbolicFactorization ();
840 verbose_cout <<
"Starting numeric factorization" << std::endl;
841 ierr = solver->NumericFactorization ();
844 verbose_cout <<
"Starting solve" << std::endl;
845 ierr = solver->Solve();
867 (
new Epetra_LinearProblem(const_cast<Epetra_CrsMatrix *>(&A.
trilinos_matrix()),
879 const ::Vector<double> &b)
889 ExcMessage (
"Can only work in serial when using deal.II vectors."));
891 Epetra_Vector ep_b (View, A.
range_partitioner(),
const_cast<double *
>(b.begin()));
907 const ::LinearAlgebra::distributed::Vector<double> &b)
911 AssertDimension (static_cast<TrilinosWrappers::types::int_type>(b.local_size()),
914 Epetra_Vector ep_b (View, A.
range_partitioner(),
const_cast<double *
>(b.begin()));
939 DEAL_II_NAMESPACE_CLOSE
941 #endif // DEAL_II_WITH_PETSC void set_preconditioner(AztecOO &solver, const Preconditioner &preconditioner)
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
bool output_solver_details
#define AssertDimension(dim1, dim2)
virtual State check(const unsigned int step, const double check_value)
std_cxx11::shared_ptr< Epetra_LinearProblem > linear_problem
const Epetra_CrsMatrix & trilinos_matrix() const
const bool output_solver_details
const AdditionalData additional_data
const AdditionalData additional_data
std_cxx11::shared_ptr< Epetra_LinearProblem > linear_problem
AdditionalData(const bool output_solver_details=false, const unsigned int restart_parameter=30)
const AdditionalData additional_data
void initialize(const SparseMatrix &A)
static::ExceptionBase & ExcTrilinosError(int arg1)
bool output_solver_details
#define AssertThrow(cond, exc)
const AdditionalData additional_data
void solve(const SparseMatrix &A, VectorBase &x, const VectorBase &b, const PreconditionBase &preconditioner)
size_type local_size() const
void do_solve(const Preconditioner &preconditioner)
SolverControl & solver_control
static::ExceptionBase & ExcMessage(std::string arg1)
bool output_solver_details
std_cxx11::shared_ptr< Epetra_Operator > preconditioner
bool output_solver_details
double last_value() const
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.
unsigned int max_steps() const
SolverBase(SolverControl &cn)
SolverControl & control() const
#define Assert(cond, exc)
AdditionalData(const bool output_solver_details=false)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
SolverBicgstab(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverControl & solver_control
bool output_solver_details
std_cxx11::shared_ptr< AztecOO_StatusTest > status_test
unsigned int last_step() const
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData(const bool output_solver_details=false)
const Epetra_Map & range_partitioner() const 1
void solve(VectorBase &x, const VectorBase &b)
const Epetra_Map & domain_partitioner() const 1
const AdditionalData additional_data
AdditionalData(const bool output_solver_details=false, const std::string &solver_type="Amesos_Klu")
SolverDirect(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
bool output_solver_details
SolverControl & control() const
AdditionalData(const bool output_solver_details=false)
static::ExceptionBase & ExcNotImplemented()
unsigned int restart_parameter
const Epetra_MultiVector & trilinos_vector() const
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
std::pair< size_type, size_type > local_range() const
static::ExceptionBase & ExcInternalError()
static::ExceptionBase & ExcTrilinosError(int arg1)
AdditionalData(const bool output_solver_details=false)
const unsigned int gmres_restart_parameter