# include <deal.II/lac/trilinos_vector_base.h>
# include <deal.II/lac/trilinos_precondition.h>
+DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
+# include <AztecOO_StatusTest.h>
+# include <AztecOO_StatusTestMaxIters.h>
+# include <AztecOO_StatusTestResNorm.h>
+# include <AztecOO_StatusTestCombo.h>
+# include <AztecOO_StatusType.h>
+DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
+
# include <cmath>
+# include <limits>
DEAL_II_NAMESPACE_OPEN
}
+ 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<AztecOO_StatusTestCombo> status_test_collection;
+ std_cxx11::shared_ptr<AztecOO_StatusTestMaxIters> status_test_max_steps;
+ std_cxx11::shared_ptr<AztecOO_StatusTestResNorm> status_test_abs_tol;
+ std_cxx11::shared_ptr<AztecOO_StatusTestResNorm> 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<double>::max()),
+ current_residual(std::numeric_limits<double>::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<typename Preconditioner>
void
SolverBase::do_solve(const Preconditioner &preconditioner)
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<const ReductionControl *const>(&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());
// 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<const internal::TrilinosReductionControl *const>(status_test.get()))
+ {
+ Assert(dynamic_cast<const ReductionControl *const>(&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(),
projects / dealii.git / commitdiff