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/base/std_cxx14/memory.h> 23 # include <deal.II/lac/trilinos_precondition.h> 24 # include <deal.II/lac/trilinos_sparse_matrix.h> 25 # include <deal.II/lac/trilinos_vector.h> 27 # include <AztecOO_StatusTest.h> 28 # include <AztecOO_StatusTestCombo.h> 29 # include <AztecOO_StatusTestMaxIters.h> 30 # include <AztecOO_StatusTestResNorm.h> 31 # include <AztecOO_StatusType.h> 36 DEAL_II_NAMESPACE_OPEN
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)
104 const_cast<Epetra_Operator *
>(&A),
119 const Epetra_Operator &preconditioner)
124 const_cast<Epetra_Operator *
>(&A),
137 Epetra_MultiVector & x,
138 const Epetra_MultiVector &b,
144 const_cast<Epetra_Operator *
>(&A),
146 const_cast<Epetra_MultiVector *>(&b));
157 Epetra_MultiVector & x,
158 const Epetra_MultiVector &b,
159 const Epetra_Operator & preconditioner)
164 const_cast<Epetra_Operator *
>(&A),
166 const_cast<Epetra_MultiVector *>(&b));
176 const ::Vector<double> &b,
184 ExcMessage(
"Can only work in serial when using deal.II vectors."));
186 ExcMessage(
"Matrix is not compressed. Call compress() method."));
189 Epetra_Vector ep_b(View,
191 const_cast<double *
>(b.begin()));
206 const ::Vector<double> &b,
209 Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.
begin());
210 Epetra_Vector ep_b(View,
211 A.OperatorRangeMap(),
212 const_cast<double *
>(b.begin()));
217 std_cxx14::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
227 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,
262 A.OperatorDomainMap().NumMyElements());
265 A.OperatorRangeMap().NumMyElements());
267 Epetra_Vector ep_x(View, A.OperatorDomainMap(), x.
begin());
268 Epetra_Vector ep_b(View,
269 A.OperatorRangeMap(),
270 const_cast<double *
>(b.begin()));
275 std_cxx14::make_unique<Epetra_LinearProblem>(&A, &ep_x, &ep_b);
286 compute_residual(
const Epetra_MultiVector *
const residual_vector)
288 Assert(residual_vector->NumVectors() == 1,
289 ExcMessage(
"Residual multivector holds more than one vector"));
290 TrilinosScalar res_l2_norm = 0.0;
291 const int ierr = residual_vector->Norm2(&res_l2_norm);
296 class TrilinosReductionControl :
public AztecOO_StatusTest
299 TrilinosReductionControl(
const int & max_steps,
300 const double & tolerance,
301 const double & reduction,
304 virtual ~TrilinosReductionControl()
override =
default;
307 ResidualVectorRequired()
const override 309 return status_test_collection->ResidualVectorRequired();
312 virtual AztecOO_StatusType
313 CheckStatus(
int CurrentIter,
314 Epetra_MultiVector *CurrentResVector,
315 double CurrentResNormEst,
316 bool SolutionUpdated)
override 321 (CurrentResNormEst < 0.0 ? compute_residual(CurrentResVector) :
323 if (CurrentIter == 0)
324 initial_residual = current_residual;
326 return status_test_collection->CheckStatus(CurrentIter,
332 virtual AztecOO_StatusType
333 GetStatus()
const override 335 return status_test_collection->GetStatus();
338 virtual std::ostream &
339 Print(std::ostream &stream,
int indent = 0)
const override 341 return status_test_collection->Print(stream, indent);
345 get_initial_residual()
const 347 return initial_residual;
351 get_current_residual()
const 353 return current_residual;
357 double initial_residual;
358 double current_residual;
359 std::unique_ptr<AztecOO_StatusTestCombo> status_test_collection;
360 std::unique_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
361 std::unique_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
362 std::unique_ptr<AztecOO_StatusTestResNorm> status_test_rel_tol;
366 TrilinosReductionControl::TrilinosReductionControl(
367 const int & max_steps,
368 const double & tolerance,
369 const double & reduction,
371 : initial_residual(std::numeric_limits<double>::max())
372 , current_residual(std::numeric_limits<double>::max())
376 status_test_collection =
377 std_cxx14::make_unique<AztecOO_StatusTestCombo>(
378 AztecOO_StatusTestCombo::OR);
382 status_test_max_steps =
383 std_cxx14::make_unique<AztecOO_StatusTestMaxIters>(max_steps);
384 status_test_collection->AddStatusTest(*status_test_max_steps);
386 Assert(linear_problem.GetRHS()->NumVectors() == 1,
387 ExcMessage(
"RHS multivector holds more than one vector"));
390 status_test_abs_tol = std_cxx14::make_unique<AztecOO_StatusTestResNorm>(
391 *linear_problem.GetOperator(),
392 *(linear_problem.GetLHS()->operator()(0)),
393 *(linear_problem.GetRHS()->operator()(0)),
395 status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
396 AztecOO_StatusTestResNorm::TwoNorm);
397 status_test_abs_tol->DefineScaleForm(
398 AztecOO_StatusTestResNorm::None, AztecOO_StatusTestResNorm::TwoNorm);
399 status_test_collection->AddStatusTest(*status_test_abs_tol);
402 status_test_rel_tol = std_cxx14::make_unique<AztecOO_StatusTestResNorm>(
403 *linear_problem.GetOperator(),
404 *(linear_problem.GetLHS()->operator()(0)),
405 *(linear_problem.GetRHS()->operator()(0)),
407 status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit,
408 AztecOO_StatusTestResNorm::TwoNorm);
409 status_test_rel_tol->DefineScaleForm(
410 AztecOO_StatusTestResNorm::NormOfInitRes,
411 AztecOO_StatusTestResNorm::TwoNorm);
412 status_test_collection->AddStatusTest(*status_test_rel_tol);
419 template <
typename Preconditioner>
432 solver.SetAztecOption(AZ_solver, AZ_cg);
435 solver.SetAztecOption(AZ_solver, AZ_cgs);
438 solver.SetAztecOption(AZ_solver, AZ_gmres);
439 solver.SetAztecOption(AZ_kspace,
443 solver.SetAztecOption(AZ_solver, AZ_bicgstab);
446 solver.SetAztecOption(AZ_solver, AZ_tfqmr);
456 solver.SetAztecOption(AZ_output,
459 solver.SetAztecOption(AZ_conv, AZ_noscaled);
480 std_cxx14::make_unique<internal::TrilinosReductionControl>(
481 reduction_control->max_steps(),
482 reduction_control->tolerance(),
483 reduction_control->reduction(),
501 "option not implemented"));
506 "numerical breakdown"));
511 "loss of precision"));
516 "GMRES Hessenberg ill-conditioned"));
529 if (
const internal::TrilinosReductionControl
530 *
const reduction_control_status =
531 dynamic_cast<const internal::TrilinosReductionControl *const>(
540 if (
solver.NumIters() > 0)
546 0, reduction_control_status->get_initial_residual());
549 reduction_control_status->get_current_residual());
554 reduction_control_status->get_current_residual());
579 const int ierr = solver.SetPrecOperator(
580 const_cast<Epetra_Operator *>(preconditioner.
preconditioner.get()));
584 solver.SetAztecOption(AZ_precond, AZ_none);
591 const Epetra_Operator &preconditioner)
594 solver.SetPrecOperator(const_cast<Epetra_Operator *>(&preconditioner));
618 const bool output_solver_details,
619 const unsigned int restart_parameter)
636 const bool output_solver_details)
686 const std::string &solver_type)
687 : output_solver_details(output_solver_details)
688 , solver_type(solver_type)
734 std::string(
"You tried to select the solver type <") +
736 "> but this solver is not supported by Trilinos either " 737 "because it does not exist, or because Trilinos was not " 738 "configured for its use."));
743 verbose_cout <<
"Starting symbolic factorization" << std::endl;
744 ierr = solver->SymbolicFactorization();
747 verbose_cout <<
"Starting numeric factorization" << std::endl;
748 ierr = solver->NumericFactorization();
769 verbose_cout <<
"Starting solve" << std::endl;
772 int ierr = solver->Solve();
784 const ::LinearAlgebra::distributed::Vector<double> &b)
786 Epetra_Vector ep_x(View,
789 Epetra_Vector ep_b(View,
791 const_cast<double *
>(b.begin()));
804 verbose_cout <<
"Starting solve" << std::endl;
807 int ierr = solver->Solve();
834 std::string(
"You tried to select the solver type <") +
836 "> but this solver is not supported by Trilinos either " 837 "because it does not exist, or because Trilinos was not " 838 "configured for its use."));
843 verbose_cout <<
"Starting symbolic factorization" << std::endl;
844 ierr = solver->SymbolicFactorization();
847 verbose_cout <<
"Starting numeric factorization" << std::endl;
848 ierr = solver->NumericFactorization();
851 verbose_cout <<
"Starting solve" << std::endl;
852 ierr = solver->Solve();
887 const ::Vector<double> &b)
894 ExcMessage(
"Can only work in serial when using deal.II vectors."));
896 Epetra_Vector ep_b(View,
898 const_cast<double *
>(b.begin()));
914 const ::LinearAlgebra::distributed::Vector<double> &b)
923 Epetra_Vector ep_b(View,
925 const_cast<double *
>(b.begin()));
948 DEAL_II_NAMESPACE_CLOSE
950 #endif // DEAL_II_WITH_PETSC void set_preconditioner(AztecOO &solver, const Preconditioner &preconditioner)
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
const Epetra_Map & domain_partitioner() const
#define AssertDimension(dim1, dim2)
virtual State check(const unsigned int step, const double check_value)
const Epetra_CrsMatrix & trilinos_matrix() const
const bool output_solver_details
const AdditionalData additional_data
std::unique_ptr< AztecOO_StatusTest > status_test
const AdditionalData additional_data
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)
void solve(MPI::Vector &x, const MPI::Vector &b)
#define AssertThrow(cond, exc)
const AdditionalData additional_data
size_type local_size() const
void do_solve(const Preconditioner &preconditioner)
SolverControl & solver_control
static::ExceptionBase & ExcMessage(std::string arg1)
double last_value() const
SolverCGS(SolverControl &cn, const AdditionalData &data=AdditionalData())
Stop iteration, goal reached.
unsigned int max_steps() const
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
std::unique_ptr< Epetra_LinearProblem > linear_problem
void solve(const SparseMatrix &A, MPI::Vector &x, const MPI::Vector &b, const PreconditionBase &preconditioner)
std::unique_ptr< Epetra_LinearProblem > linear_problem
unsigned int last_step() const
SolverTFQMR(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData(const bool output_solver_details=false)
SolverBase(SolverControl &cn, const AdditionalData &data=AdditionalData())
const Epetra_Map & range_partitioner() const
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()
AdditionalData(const bool output_solver_details=false, const unsigned int gmres_restart_parameter=30)
std::shared_ptr< Epetra_Operator > preconditioner
SolverGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
const AdditionalData additional_data
std::pair< size_type, size_type > local_range() const
const Epetra_MultiVector & trilinos_vector() const
static::ExceptionBase & ExcInternalError()
static::ExceptionBase & ExcTrilinosError(int arg1)
AdditionalData(const bool output_solver_details=false)
const unsigned int gmres_restart_parameter