solver_relaxation.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2010 - 2017 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.
//
// ---------------------------------------------------------------------

#ifndef dealii_solver_relaxation_h
#define dealii_solver_relaxation_h


#include <deal.II/base/config.h>

#include <deal.II/base/logstream.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_control.h>

DEAL_II_NAMESPACE_OPEN

template <typename VectorType = Vector<double>>
class SolverRelaxation : public Solver<VectorType>
{
public:
  struct AdditionalData
  {};

  SolverRelaxation(SolverControl &       cn,
                   const AdditionalData &data = AdditionalData());

  virtual ~SolverRelaxation();

  template <typename MatrixType, class RelaxationType>
  void
  solve(const MatrixType &    A,
        VectorType &          x,
        const VectorType &    b,
        const RelaxationType &R);
};

//----------------------------------------------------------------------//

template <class VectorType>
SolverRelaxation<VectorType>::SolverRelaxation(SolverControl &cn,
                                               const AdditionalData &)
  : Solver<VectorType>(cn)
{}



template <class VectorType>
SolverRelaxation<VectorType>::~SolverRelaxation()
{}


template <class VectorType>
template <typename MatrixType, class RelaxationType>
void
SolverRelaxation<VectorType>::solve(const MatrixType &    A,
                                    VectorType &          x,
                                    const VectorType &    b,
                                    const RelaxationType &R)
{
  GrowingVectorMemory<VectorType> mem;
  SolverControl::State            conv = SolverControl::iterate;

  // Memory allocation
  typename VectorMemory<VectorType>::Pointer Vr(mem);
  VectorType &                               r = *Vr;
  r.reinit(x);
  typename VectorMemory<VectorType>::Pointer Vd(mem);
  VectorType &                               d = *Vd;
  d.reinit(x);

  LogStream::Prefix prefix("Relaxation");

  int iter = 0;
  // Main loop
  for (; conv == SolverControl::iterate; iter++)
    {
      // Compute residual
      A.vmult(r, x);
      r.sadd(-1., 1., b);

      // The required norm of the
      // (preconditioned)
      // residual is computed in
      // criterion() and stored
      // in res.
      conv = this->iteration_status(iter, r.l2_norm(), x);
      if (conv != SolverControl::iterate)
        break;
      R.step(x, b);
    }

  // in case of failure: throw exception
  AssertThrow(conv == SolverControl::success,
              SolverControl::NoConvergence(iter, r.l2_norm()));
  // otherwise exit as normal
}


DEAL_II_NAMESPACE_CLOSE

#endif