From 54941c1c64ba2a754f94bb6f036e2741f04b9006 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Sat, 28 Jan 2017 18:34:56 +0100 Subject: [PATCH] Implement Trilinos AztecOO_StatusTest for ReductionControl. This patch ensures that Trilinos solvers now respect the convergence criterion specified by ReductionControl. To do so, we create an AztecOO_StatusTest that monitors the heuristics that ReductionControl uses to test for convergence. Fixes #3843 --- .../changes/minor/20170128Jean-PaulPelteret | 4 + include/deal.II/lac/trilinos_solver.h | 6 + source/lac/trilinos_solver.cc | 203 +++++++++++++++++- 3 files changed, 212 insertions(+), 1 deletion(-) create mode 100644 doc/news/changes/minor/20170128Jean-PaulPelteret diff --git a/doc/news/changes/minor/20170128Jean-PaulPelteret b/doc/news/changes/minor/20170128Jean-PaulPelteret new file mode 100644 index 0000000000..b985ce1a2b --- /dev/null +++ b/doc/news/changes/minor/20170128Jean-PaulPelteret @@ -0,0 +1,4 @@ +Fixed: Trilinos solvers now respect the convergence criterion specified by +ReductionControl. +
+(Jean-Paul Pelteret, 2017/01/28) diff --git a/include/deal.II/lac/trilinos_solver.h b/include/deal.II/lac/trilinos_solver.h index 2ec95a4a85..6bf091828d 100644 --- a/include/deal.II/lac/trilinos_solver.h +++ b/include/deal.II/lac/trilinos_solver.h @@ -323,6 +323,12 @@ namespace TrilinosWrappers */ std_cxx11::shared_ptr linear_problem; + /** + * A structure that contains a Trilinos object that can query the linear + * solver and determine whether the convergence criterion have been met. + */ + std_cxx11::shared_ptr status_test; + /** * A structure that contains the Trilinos solver and preconditioner * objects. diff --git a/source/lac/trilinos_solver.cc b/source/lac/trilinos_solver.cc index 400b7df1ee..bc4bafa3ce 100644 --- a/source/lac/trilinos_solver.cc +++ b/source/lac/trilinos_solver.cc @@ -22,7 +22,16 @@ # include # include +DEAL_II_DISABLE_EXTRA_DIAGNOSTICS +# include +# include +# include +# include +# include +DEAL_II_ENABLE_EXTRA_DIAGNOSTICS + # include +# include DEAL_II_NAMESPACE_OPEN @@ -282,6 +291,144 @@ namespace TrilinosWrappers } + namespace internal + { + namespace + { + double + compute_residual (const Epetra_MultiVector *const residual_vector) + { + Assert(residual_vector->NumVectors() == 1, + ExcMessage("Residual multivector holds more than one vector")); + TrilinosScalar res_l2_norm = 0.0; + const int ierr = residual_vector->Norm2 (&res_l2_norm); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + return res_l2_norm; + } + + class TrilinosReductionControl : public AztecOO_StatusTest + { + public: + TrilinosReductionControl (const double &max_steps, + const double &tolerance, + const double &reduction, + const Epetra_LinearProblem &linear_problem); + + virtual ~TrilinosReductionControl() {} + + virtual bool + ResidualVectorRequired () const + { + return status_test_collection->ResidualVectorRequired(); + } + + virtual AztecOO_StatusType + CheckStatus (int CurrentIter, + Epetra_MultiVector *CurrentResVector, + double CurrentResNormEst, + bool SolutionUpdated) + { + // Note: CurrentResNormEst is set to -1.0 if no estimate of the + // residual value is available + current_residual = (CurrentResNormEst < 0.0 ? + compute_residual(CurrentResVector) : + CurrentResNormEst); + if (CurrentIter == 0) + initial_residual = current_residual; + + return status_test_collection->CheckStatus(CurrentIter, + CurrentResVector, + CurrentResNormEst, + SolutionUpdated); + + } + + virtual AztecOO_StatusType + GetStatus () const + { + return status_test_collection->GetStatus(); + } + + virtual std::ostream & + Print (std::ostream &stream, + int indent = 0) const + { + return status_test_collection->Print(stream,indent); + } + + double + get_initial_residual() const + { + return initial_residual; + } + + double + get_current_residual() const + { + return current_residual; + } + + private: + double initial_residual; + double current_residual; + std_cxx11::shared_ptr status_test_collection; + std_cxx11::shared_ptr status_test_max_steps; + std_cxx11::shared_ptr status_test_abs_tol; + std_cxx11::shared_ptr status_test_rel_tol; + }; + + + TrilinosReductionControl::TrilinosReductionControl( + const double &max_steps, + const double &tolerance, + const double &reduction, + const Epetra_LinearProblem &linear_problem ) + : + initial_residual (std::numeric_limits::max()), + current_residual(std::numeric_limits::max()) + { + // Consider linear problem converged if any of the collection + // of criterion are met + status_test_collection.reset( + new AztecOO_StatusTestCombo (AztecOO_StatusTestCombo::OR) ); + + // Maximum number of iterations + status_test_max_steps.reset( + new AztecOO_StatusTestMaxIters(max_steps) ); + status_test_collection->AddStatusTest(*status_test_max_steps); + + Assert(linear_problem.GetRHS()->NumVectors() == 1, + ExcMessage("RHS multivector holds more than one vector")); + + // Residual norm is below some absolute value + status_test_abs_tol.reset( + new AztecOO_StatusTestResNorm(*linear_problem.GetOperator(), + *(linear_problem.GetLHS()->operator()(0)), + *(linear_problem.GetRHS()->operator()(0)), + tolerance) ); + status_test_abs_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit, + AztecOO_StatusTestResNorm::TwoNorm); + status_test_abs_tol->DefineScaleForm(AztecOO_StatusTestResNorm::None, + AztecOO_StatusTestResNorm::TwoNorm); + status_test_collection->AddStatusTest(*status_test_abs_tol); + + // Residual norm, scaled by some initial value, is below some threshold + status_test_rel_tol.reset( + new AztecOO_StatusTestResNorm(*linear_problem.GetOperator(), + *(linear_problem.GetLHS()->operator()(0)), + *(linear_problem.GetRHS()->operator()(0)), + reduction) ); + status_test_rel_tol->DefineResForm(AztecOO_StatusTestResNorm::Explicit, + AztecOO_StatusTestResNorm::TwoNorm); + status_test_rel_tol->DefineScaleForm(AztecOO_StatusTestResNorm::NormOfInitRes, + AztecOO_StatusTestResNorm::TwoNorm); + status_test_collection->AddStatusTest(*status_test_rel_tol); + } + + } + } + + template void SolverBase::do_solve(const Preconditioner &preconditioner) @@ -322,6 +469,33 @@ namespace TrilinosWrappers AZ_all : AZ_none); solver.SetAztecOption (AZ_conv, AZ_noscaled); + // By default, the Trilinos solver chooses convergence criterion based on + // the number of iterations made and an absolute tolerance. + // This implies that the use of the standard Trilinos convergence test + // actually coincides with dealii::IterationNumberControl because the + // solver, unless explicitly told otherwise, will Iterate() until a number + // of max_steps() are taken or an absolute tolerance() is attained. + // It is therefore suitable for use with both SolverControl or + // IterationNumberControl. The final check at the end will determine whether + // failure to converge to the defined residual norm constitutes failure + // (SolverControl) or is alright (IterationNumberControl). + // In the case that the SolverControl wants to perform ReductionControl, + // then we have to do a little extra something by prescribing a custom + // status test. + if (!status_test) + { + if (const ReductionControl* const reduction_control + = dynamic_cast(&solver_control)) + { + status_test.reset(new internal::TrilinosReductionControl( + reduction_control->max_steps(), + reduction_control->tolerance(), + reduction_control->reduction(), + *linear_problem) ); + solver.SetStatusTest(status_test.get()); + } + } + // ... and then solve! ierr = solver.Iterate (solver_control.max_steps(), solver_control.tolerance()); @@ -350,7 +524,34 @@ namespace TrilinosWrappers // Finally, let the deal.II SolverControl object know what has // happened. If the solve succeeded, the status of the solver control will // turn into SolverControl::success. - solver_control.check (solver.NumIters(), solver.TrueResidual()); + // If the residual is not computed/stored by the solver, as can happen for + // certain choices of solver or if a custom status test is set, then the + // result returned by TrueResidual() is equal to -1. In this case we must + // compute it ourself. + if (const internal::TrilinosReductionControl* const reduction_control_status + = dynamic_cast(status_test.get())) + { + Assert(dynamic_cast(&solver_control), ExcInternalError()); + + // Check to see if solver converged in one step + // This can happen if the matrix is diagonal and a non-trivial + // preconditioner is used. + if (solver.NumIters() > 0) + { + // For ReductionControl, we must first register the initial residual + // value. This is the basis from which it will determine whether the + // current residual corresponds to a converged state. + solver_control.check (0, reduction_control_status->get_initial_residual()); + solver_control.check (solver.NumIters(), reduction_control_status->get_current_residual()); + } + else + solver_control.check (solver.NumIters(), reduction_control_status->get_current_residual()); + } + else + { + Assert(solver.TrueResidual() >= 0.0, ExcInternalError()); + solver_control.check (solver.NumIters(), solver.TrueResidual()); + } if (solver_control.last_check() != SolverControl::success) AssertThrow(false, SolverControl::NoConvergence (solver_control.last_step(), -- 2.39.5