Clean up implementation in SolverBicgstab

diff --git a/include/deal.II/lac/solver_bicgstab.h b/include/deal.II/lac/solver_bicgstab.h
index f47e8cfe9460fabf0610f7547ade5f6d413cf37b..681910c7ae35bb362c9ca547c312cd2c1c6ebba2 100644 (file)
--- 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 <typename VectorType = Vector<double>>
-class SolverBicgstab : public SolverBase<VectorType>,
-                       protected internal::SolverBicgstabData
+class SolverBicgstab : public SolverBase<VectorType>
 {
 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<VectorType>::Pointer Vr;
-
-  /**
-   * Auxiliary vector.
-   */
-  typename VectorMemory<VectorType>::Pointer Vrbar;
-
-  /**
-   * Auxiliary vector.
-   */
-  typename VectorMemory<VectorType>::Pointer Vp;
-
-  /**
-   * Auxiliary vector.
-   */
-  typename VectorMemory<VectorType>::Pointer Vy;
-
-  /**
-   * Auxiliary vector.
-   */
-  typename VectorMemory<VectorType>::Pointer Vz;
-
-  /**
-   * Auxiliary vector.
-   */
-  typename VectorMemory<VectorType>::Pointer Vt;
-
-  /**
-   * Auxiliary vector.
-   */
-  typename VectorMemory<VectorType>::Pointer Vv;
-
-  /**
-   * A pointer to the right hand side vector passed to solve().
-   */
-  const VectorType *Vb;
-
   /**
    * Computation of the stopping criterion.
    */
   template <typename MatrixType>
   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 <typename MatrixType, typename PreconditionerType>
   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<VectorType>::SolverBicgstab(SolverControl &           cn,
                                            VectorMemory<VectorType> &mem,
                                            const AdditionalData &    data)
   : SolverBase<VectorType>(cn, mem)
-  , Vx(nullptr)
-  , Vb(nullptr)
   , additional_data(data)
 {}
 
@@ -325,8 +237,6 @@ template <typename VectorType>
 SolverBicgstab<VectorType>::SolverBicgstab(SolverControl &       cn,
                                            const AdditionalData &data)
   : SolverBase<VectorType>(cn)
-  , Vx(nullptr)
-  , Vb(nullptr)
   , additional_data(data)
 {}
 
@@ -337,13 +247,11 @@ template <typename MatrixType>
 double
 SolverBicgstab<VectorType>::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 <typename VectorType>
 template <typename MatrixType, typename PreconditionerType>
 typename SolverBicgstab<VectorType>::IterationResult
 SolverBicgstab<VectorType>::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<VectorType>::Pointer Vr(this->memory);
+  typename VectorMemory<VectorType>::Pointer Vrbar(this->memory);
+  typename VectorMemory<VectorType>::Pointer Vp(this->memory);
+  typename VectorMemory<VectorType>::Pointer Vy(this->memory);
+  typename VectorMemory<VectorType>::Pointer Vz(this->memory);
+  typename VectorMemory<VectorType>::Pointer Vt(this->memory);
+  typename VectorMemory<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::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<VectorType>::solve(const MatrixType &        A,
 {
   LogStream::Prefix prefix("Bicgstab");
 
-  // Allocate temporary memory.
-  Vr    = typename VectorMemory<VectorType>::Pointer(this->memory);
-  Vrbar = typename VectorMemory<VectorType>::Pointer(this->memory);
-  Vp    = typename VectorMemory<VectorType>::Pointer(this->memory);
-  Vy    = typename VectorMemory<VectorType>::Pointer(this->memory);
-  Vz    = typename VectorMemory<VectorType>::Pointer(this->memory);
-  Vt    = typename VectorMemory<VectorType>::Pointer(this->memory);
-  Vv    = typename VectorMemory<VectorType>::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 cf383afbce8962b06b468a455652010a3f326da5..fda0205828aa2c80a37d8f069337a0c66ba9f869 100644 (file)
--- 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 (file)
index 16dd36e..0000000
--- 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/lac/solver_bicgstab.h>
-
-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<double>())
-{}
-
-
-DEAL_II_NAMESPACE_CLOSE