From 74ddd199f3e3932aa4df8ad3916644d32f3b08de Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 12 Sep 2022 17:04:05 +0200 Subject: [PATCH] Clean up implementation in SolverBicgstab --- include/deal.II/lac/solver_bicgstab.h | 245 ++++++++------------------ source/lac/CMakeLists.txt | 1 - source/lac/solver_bicgstab.cc | 33 ---- 3 files changed, 73 insertions(+), 206 deletions(-) delete mode 100644 source/lac/solver_bicgstab.cc diff --git a/include/deal.II/lac/solver_bicgstab.h b/include/deal.II/lac/solver_bicgstab.h index f47e8cfe94..681910c7ae 100644 --- a/include/deal.II/lac/solver_bicgstab.h +++ b/include/deal.II/lac/solver_bicgstab.h @@ -36,54 +36,6 @@ DEAL_II_NAMESPACE_OPEN * @{ */ -namespace internal -{ - /** - * Class containing the non-parameter non-template values used by the - * SolverBicgstab class. - */ - class SolverBicgstabData - { - protected: - /** - * Auxiliary value. - */ - double alpha; - /** - * Auxiliary value. - */ - double beta; - /** - * Auxiliary value. - */ - double omega; - /** - * Auxiliary value. - */ - double rho; - /** - * Auxiliary value. - */ - double rhobar; - - /** - * Current iteration step. - */ - unsigned int step; - - /** - * Residual. - */ - double res; - - /** - * Default constructor. This is protected so that only SolverBicgstab can - * create instances. - */ - SolverBicgstabData(); - }; -} // namespace internal - /** * Bicgstab algorithm by van der Vorst. * @@ -124,8 +76,7 @@ namespace internal * to observe the progress of the iteration. */ template > -class SolverBicgstab : public SolverBase, - protected internal::SolverBicgstabData +class SolverBicgstab : public SolverBase { public: /** @@ -194,57 +145,15 @@ public: const PreconditionerType &preconditioner); protected: - /** - * A pointer to the solution vector passed to solve(). - */ - VectorType *Vx; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vr; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vrbar; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vp; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vy; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vz; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vt; - - /** - * Auxiliary vector. - */ - typename VectorMemory::Pointer Vv; - - /** - * A pointer to the right hand side vector passed to solve(). - */ - const VectorType *Vb; - /** * Computation of the stopping criterion. */ template double - criterion(const MatrixType &A, const VectorType &x, const VectorType &b); + criterion(const MatrixType &A, + const VectorType &x, + const VectorType &b, + VectorType & t); /** * Interface for derived class. This function gets the current iteration @@ -286,9 +195,14 @@ private: */ template IterationResult - iterate(const MatrixType &A, const PreconditionerType &preconditioner); + iterate(const MatrixType & A, + VectorType & x, + const VectorType & b, + const PreconditionerType &preconditioner, + const unsigned int step); }; + /** @} */ /*-------------------------Inline functions -------------------------------*/ @@ -314,8 +228,6 @@ SolverBicgstab::SolverBicgstab(SolverControl & cn, VectorMemory &mem, const AdditionalData & data) : SolverBase(cn, mem) - , Vx(nullptr) - , Vb(nullptr) , additional_data(data) {} @@ -325,8 +237,6 @@ template SolverBicgstab::SolverBicgstab(SolverControl & cn, const AdditionalData &data) : SolverBase(cn) - , Vx(nullptr) - , Vb(nullptr) , additional_data(data) {} @@ -337,13 +247,11 @@ template double SolverBicgstab::criterion(const MatrixType &A, const VectorType &x, - const VectorType &b) + const VectorType &b, + VectorType & t) { - A.vmult(*Vt, x); - Vt->add(-1., b); - res = Vt->l2_norm(); - - return res; + A.vmult(t, x); + return std::sqrt(t.add_and_dot(-1.0, b, t)); } @@ -362,18 +270,21 @@ template template typename SolverBicgstab::IterationResult SolverBicgstab::iterate(const MatrixType & A, - const PreconditionerType &preconditioner) + VectorType & x, + const VectorType & b, + const PreconditionerType &preconditioner, + const unsigned int last_step) { - 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.; - + // Allocate temporary memory. + typename VectorMemory::Pointer Vr(this->memory); + typename VectorMemory::Pointer Vrbar(this->memory); + typename VectorMemory::Pointer Vp(this->memory); + typename VectorMemory::Pointer Vy(this->memory); + typename VectorMemory::Pointer Vz(this->memory); + typename VectorMemory::Pointer Vt(this->memory); + typename VectorMemory::Pointer Vv(this->memory); + + // Define a few aliases for simpler use of the vectors VectorType &r = *Vr; VectorType &rbar = *Vrbar; VectorType &p = *Vp; @@ -382,24 +293,46 @@ SolverBicgstab::iterate(const MatrixType & A, VectorType &t = *Vt; VectorType &v = *Vv; - rbar = r; - bool startup = true; + r.reinit(x, true); + rbar.reinit(x, true); + p.reinit(x, true); + y.reinit(x, true); + z.reinit(x, true); + t.reinit(x, true); + v.reinit(x, true); + + A.vmult(r, x); + r.sadd(-1., 1., b); + double res = r.l2_norm(); + + unsigned int step = last_step; + + SolverControl::State state = this->iteration_status(step, res, x); + if (state == SolverControl::State::success) + return IterationResult(false, state, step, res); + + rbar = r; + + double alpha = 1.; + double rho = 1.; + double omega = 1.; do { ++step; - rhobar = r * rbar; + const double rhobar = (step == 1 + last_step) ? res * res : r * rbar; + if (std::fabs(rhobar) < additional_data.breakdown) { return IterationResult(true, state, step, res); } - beta = rhobar * alpha / (rho * omega); - rho = rhobar; - if (startup == true) + + const double beta = rhobar * alpha / (rho * omega); + rho = rhobar; + if (step == last_step + 1) { - p = r; - startup = false; + p = r; } else { @@ -409,13 +342,13 @@ SolverBicgstab::iterate(const MatrixType & A, preconditioner.vmult(y, p); A.vmult(v, y); - rhobar = rbar * v; - if (std::fabs(rhobar) < additional_data.breakdown) + const double rbar_dot_v = rbar * v; + if (std::fabs(rbar_dot_v) < additional_data.breakdown) { return IterationResult(true, state, step, res); } - alpha = rho / rhobar; + alpha = rho / rbar_dot_v; res = std::sqrt(r.add_and_dot(-alpha, v, r)); @@ -425,34 +358,34 @@ SolverBicgstab::iterate(const MatrixType & A, // note: the vector *Vx we pass to the iteration_status signal here is // only the current approximation, not the one we will return with, which // will be x=*Vx + alpha*y - if (this->iteration_status(step, res, *Vx) == SolverControl::success) + if (this->iteration_status(step, res, x) == SolverControl::success) { - Vx->add(alpha, y); - print_vectors(step, *Vx, r, y); + x.add(alpha, y); + print_vectors(step, x, r, y); return IterationResult(false, SolverControl::success, step, res); } preconditioner.vmult(z, r); A.vmult(t, z); - rhobar = t * r; - auto t_squared = t * t; + const double t_dot_r = t * r; + const double 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); + omega = t_dot_r / t_squared; + x.add(alpha, y, omega, z); if (additional_data.exact_residual) { r.add(-omega, t); - res = criterion(A, *Vx, *Vb); + res = criterion(A, x, b, t); } else res = std::sqrt(r.add_and_dot(-omega, t, r)); - state = this->iteration_status(step, res, *Vx); - print_vectors(step, *Vx, r, y); + state = this->iteration_status(step, res, x); + print_vectors(step, x, r, y); } while (state == SolverControl::iterate); @@ -471,45 +404,13 @@ SolverBicgstab::solve(const MatrixType & A, { LogStream::Prefix prefix("Bicgstab"); - // Allocate temporary memory. - Vr = typename VectorMemory::Pointer(this->memory); - Vrbar = typename VectorMemory::Pointer(this->memory); - Vp = typename VectorMemory::Pointer(this->memory); - Vy = typename VectorMemory::Pointer(this->memory); - Vz = typename VectorMemory::Pointer(this->memory); - Vt = typename VectorMemory::Pointer(this->memory); - Vv = typename VectorMemory::Pointer(this->memory); - - Vr->reinit(x, true); - Vrbar->reinit(x, true); - Vp->reinit(x, true); - Vy->reinit(x, true); - Vz->reinit(x, true); - Vt->reinit(x, true); - Vv->reinit(x, true); - - Vx = &x; - Vb = &b; - - step = 0; - IterationResult state(false, SolverControl::failure, 0, 0); do { - state = iterate(A, preconditioner); + state = iterate(A, x, b, preconditioner, state.last_step); } while (state.state == SolverControl::iterate); - - // Release the temporary memory again. - Vr.reset(); - Vrbar.reset(); - Vp.reset(); - Vy.reset(); - Vz.reset(); - Vt.reset(); - Vv.reset(); - // In case of failure: throw exception AssertThrow(state.state == SolverControl::success, SolverControl::NoConvergence(state.last_step, diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index cf383afbce..fda0205828 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -35,7 +35,6 @@ SET(_unity_include_src relaxation_block.cc read_write_vector.cc solver.cc - solver_bicgstab.cc solver_control.cc sparse_decomposition.cc sparse_direct.cc diff --git a/source/lac/solver_bicgstab.cc b/source/lac/solver_bicgstab.cc deleted file mode 100644 index 16dd36ee9f..0000000000 --- a/source/lac/solver_bicgstab.cc +++ /dev/null @@ -1,33 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1998 - 2018 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE.md at -// the top level directory of deal.II. -// -// --------------------------------------------------------------------- - - -#include - -DEAL_II_NAMESPACE_OPEN - - -internal::SolverBicgstabData::SolverBicgstabData() - : alpha(0.) - , beta(0.) - , omega(0.) - , rho(0.) - , rhobar(0.) - , step(numbers::invalid_unsigned_int) - , res(numbers::signaling_nan()) -{} - - -DEAL_II_NAMESPACE_CLOSE -- 2.39.5