From 1e48462f8c7824347d83863236960b4ab71fcb77 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Wed, 9 Dec 2020 00:34:10 +0000 Subject: [PATCH] Rework breakdown condition and use input parameter. Fix infinite loop in case of breakdown --- include/deal.II/lac/solver_bicgstab.h | 81 +++++++++++---------------- 1 file changed, 33 insertions(+), 48 deletions(-) diff --git a/include/deal.II/lac/solver_bicgstab.h b/include/deal.II/lac/solver_bicgstab.h index a8705964f8..dd5daa3b71 100644 --- a/include/deal.II/lac/solver_bicgstab.h +++ b/include/deal.II/lac/solver_bicgstab.h @@ -140,10 +140,13 @@ public: * Constructor. * * The default is to perform an exact residual computation and breakdown - * parameter 1e-10. + * parameter is the minimum finite value representable by the value_type of + * VectorType. */ - explicit AdditionalData(const bool exact_residual = true, - const double breakdown = 1.e-10) + explicit AdditionalData( + const bool exact_residual = true, + const double breakdown = + std::numeric_limits::min()) : exact_residual(exact_residual) , breakdown(breakdown) {} @@ -256,13 +259,6 @@ protected: AdditionalData additional_data; private: - /** - * Everything before the iteration loop. - */ - template - SolverControl::State - start(const MatrixType &A); - /** * A structure returned by the iterate() function representing what it found * is happening during the iteration. @@ -348,20 +344,6 @@ SolverBicgstab::criterion(const MatrixType &A, -template -template -SolverControl::State -SolverBicgstab::start(const MatrixType &A) -{ - A.vmult(*Vr, *Vx); - Vr->sadd(-1., 1., *Vb); - res = Vr->l2_norm(); - - return this->iteration_status(step, res, *Vx); -} - - - template void SolverBicgstab::print_vectors(const unsigned int, @@ -378,9 +360,14 @@ typename SolverBicgstab::IterationResult SolverBicgstab::iterate(const MatrixType & A, const PreconditionerType &preconditioner) { - // TODO:[GK] Implement "use the length of the computed orthogonal residual" in - // the BiCGStab method. - SolverControl::State state = SolverControl::iterate; + A.vmult(*Vr, *Vx); + Vr->sadd(-1., 1., *Vb); + res = Vr->l2_norm(); + + SolverControl::State state = this->iteration_status(step, res, *Vx); + if (state == SolverControl::State::success) + return IterationResult(false, state, step, res); + alpha = omega = rho = 1.; VectorType &r = *Vr; @@ -399,8 +386,12 @@ SolverBicgstab::iterate(const MatrixType & A, ++step; rhobar = r * rbar; - beta = rhobar * alpha / (rho * omega); - rho = rhobar; + if (std::fabs(rhobar) < additional_data.breakdown) + { + return IterationResult(true, state, step, res); + } + beta = rhobar * alpha / (rho * omega); + rho = rhobar; if (startup == true) { p = r; @@ -415,14 +406,13 @@ SolverBicgstab::iterate(const MatrixType & A, preconditioner.vmult(y, p); A.vmult(v, y); rhobar = rbar * v; + if (std::fabs(rhobar) < additional_data.breakdown) + { + return IterationResult(true, state, step, res); + } alpha = rho / rhobar; - // TODO:[?] Find better breakdown criterion - - if (std::fabs(alpha) > 1.e10) - return IterationResult(true, state, step, res); - res = std::sqrt(r.add_and_dot(-alpha, v, r)); // check for early success, see the lac/bicgstab_early testcase as to @@ -440,8 +430,13 @@ SolverBicgstab::iterate(const MatrixType & A, preconditioner.vmult(z, r); A.vmult(t, z); - rhobar = t * r; - omega = rhobar / (t * t); + rhobar = t * r; + auto t_squared = t * t; + if (t_squared < additional_data.breakdown) + { + return IterationResult(true, state, step, res); + } + omega = rhobar / (t * t); Vx->add(alpha, y, omega, z); if (additional_data.exact_residual) @@ -495,22 +490,12 @@ SolverBicgstab::solve(const MatrixType & A, step = 0; IterationResult state(false, SolverControl::failure, 0, 0); - - // Iterate while the inner iteration returns a breakdown, i.e., try and try - // until we succeed. do { - if (step != 0) - deallog << "Restart step " << step << std::endl; - if (start(A) == SolverControl::success) - { - state.state = SolverControl::success; - break; - } state = iterate(A, preconditioner); - ++step; } - while (state.breakdown == true); + while (state.state == SolverControl::iterate); + // Release the temporary memory again. Vr.reset(); -- 2.39.5