mg_smoother.h

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

#ifndef dealii_mg_smoother_h
#define dealii_mg_smoother_h


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

#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/smartpointer.h>

#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/vector_memory.h>

#include <deal.II/multigrid/mg_base.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

/*
 * MGSmootherBase is defined in mg_base.h
 */


template <typename VectorType>
class MGSmoother : public MGSmootherBase<VectorType>
{
public:
  MGSmoother(const unsigned int steps     = 1,
             const bool         variable  = false,
             const bool         symmetric = false,
             const bool         transpose = false);

  void
  set_steps(const unsigned int);

  void
  set_variable(const bool);

  void
  set_symmetric(const bool);

  void
  set_transpose(const bool);

  void
  set_debug(const unsigned int level);

protected:
  mutable GrowingVectorMemory<VectorType> vector_memory;

  unsigned int steps;

  bool variable;

  bool symmetric;

  bool transpose;

  unsigned int debug;
};


template <typename VectorType>
class MGSmootherIdentity : public MGSmootherBase<VectorType>
{
public:
  virtual void
  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;

  virtual void
  clear();
};


namespace mg
{
  template <class RelaxationType, typename VectorType>
  class SmootherRelaxation : public MGLevelObject<RelaxationType>,
                             public MGSmoother<VectorType>
  {
  public:
    SmootherRelaxation(const unsigned int steps     = 1,
                       const bool         variable  = false,
                       const bool         symmetric = false,
                       const bool         transpose = false);

    template <typename MatrixType2>
    void
    initialize(const MGLevelObject<MatrixType2> &             matrices,
               const typename RelaxationType::AdditionalData &additional_data =
                 typename RelaxationType::AdditionalData());

    template <typename MatrixType2, class DATA>
    void
    initialize(const MGLevelObject<MatrixType2> &matrices,
               const MGLevelObject<DATA> &       additional_data);

    void
    clear() override;

    virtual void
    smooth(const unsigned int level,
           VectorType &       u,
           const VectorType & rhs) const override;

    virtual void
    apply(const unsigned int level,
          VectorType &       u,
          const VectorType & rhs) const override;

    std::size_t
    memory_consumption() const;
  };
} // namespace mg

template <typename MatrixType, class RelaxationType, typename VectorType>
class MGSmootherRelaxation : public MGSmoother<VectorType>
{
public:
  MGSmootherRelaxation(const unsigned int steps     = 1,
                       const bool         variable  = false,
                       const bool         symmetric = false,
                       const bool         transpose = false);

  template <typename MatrixType2>
  void
  initialize(const MGLevelObject<MatrixType2> &             matrices,
             const typename RelaxationType::AdditionalData &additional_data =
               typename RelaxationType::AdditionalData());

  template <typename MatrixType2, class DATA>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const MGLevelObject<DATA> &       additional_data);

  template <typename MatrixType2, class DATA>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const DATA &                      additional_data,
             const unsigned int                block_row,
             const unsigned int                block_col);

  template <typename MatrixType2, class DATA>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const MGLevelObject<DATA> &       additional_data,
             const unsigned int                block_row,
             const unsigned int                block_col);

  void
  clear();

  virtual void
  smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;

  virtual void
  apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;

  MGLevelObject<RelaxationType> smoothers;

  std::size_t
  memory_consumption() const;


private:
  MGLevelObject<LinearOperator<VectorType>> matrices;
};



template <typename MatrixType, typename PreconditionerType, typename VectorType>
class MGSmootherPrecondition : public MGSmoother<VectorType>
{
public:
  MGSmootherPrecondition(const unsigned int steps     = 1,
                         const bool         variable  = false,
                         const bool         symmetric = false,
                         const bool         transpose = false);

  template <typename MatrixType2>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const typename PreconditionerType::AdditionalData &
               additional_data = typename PreconditionerType::AdditionalData());

  template <typename MatrixType2, class DATA>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const MGLevelObject<DATA> &       additional_data);

  template <typename MatrixType2, class DATA>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const DATA &                      additional_data,
             const unsigned int                block_row,
             const unsigned int                block_col);

  template <typename MatrixType2, class DATA>
  void
  initialize(const MGLevelObject<MatrixType2> &matrices,
             const MGLevelObject<DATA> &       additional_data,
             const unsigned int                block_row,
             const unsigned int                block_col);

  void
  clear() override;

  virtual void
  smooth(const unsigned int level,
         VectorType &       u,
         const VectorType & rhs) const override;

  virtual void
  apply(const unsigned int level,
        VectorType &       u,
        const VectorType & rhs) const override;

  MGLevelObject<PreconditionerType> smoothers;

  std::size_t
  memory_consumption() const;


private:
  MGLevelObject<LinearOperator<VectorType>> matrices;
};

/* ------------------------------- Inline functions --------------------------
 */

#ifndef DOXYGEN

template <typename VectorType>
inline void
MGSmootherIdentity<VectorType>::smooth(const unsigned int,
                                       VectorType &,
                                       const VectorType &) const
{}

template <typename VectorType>
inline void
MGSmootherIdentity<VectorType>::clear()
{}

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

template <typename VectorType>
inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
                                          const bool         variable,
                                          const bool         symmetric,
                                          const bool         transpose)
  : steps(steps)
  , variable(variable)
  , symmetric(symmetric)
  , transpose(transpose)
  , debug(0)
{}


template <typename VectorType>
inline void
MGSmoother<VectorType>::set_steps(const unsigned int s)
{
  steps = s;
}


template <typename VectorType>
inline void
MGSmoother<VectorType>::set_debug(const unsigned int s)
{
  debug = s;
}


template <typename VectorType>
inline void
MGSmoother<VectorType>::set_variable(const bool flag)
{
  variable = flag;
}


template <typename VectorType>
inline void
MGSmoother<VectorType>::set_symmetric(const bool flag)
{
  symmetric = flag;
}


template <typename VectorType>
inline void
MGSmoother<VectorType>::set_transpose(const bool flag)
{
  transpose = flag;
}

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

namespace mg
{
  template <class RelaxationType, typename VectorType>
  inline SmootherRelaxation<RelaxationType, VectorType>::SmootherRelaxation(
    const unsigned int steps,
    const bool         variable,
    const bool         symmetric,
    const bool         transpose)
    : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
  {}


  template <class RelaxationType, typename VectorType>
  inline void
  SmootherRelaxation<RelaxationType, VectorType>::clear()
  {
    MGLevelObject<RelaxationType>::clear_elements();
  }


  template <class RelaxationType, typename VectorType>
  template <typename MatrixType2>
  inline void
  SmootherRelaxation<RelaxationType, VectorType>::initialize(
    const MGLevelObject<MatrixType2> &             m,
    const typename RelaxationType::AdditionalData &data)
  {
    const unsigned int min = m.min_level();
    const unsigned int max = m.max_level();

    this->resize(min, max);

    for (unsigned int i = min; i <= max; ++i)
      (*this)[i].initialize(m[i], data);
  }


  template <class RelaxationType, typename VectorType>
  template <typename MatrixType2, class DATA>
  inline void
  SmootherRelaxation<RelaxationType, VectorType>::initialize(
    const MGLevelObject<MatrixType2> &m,
    const MGLevelObject<DATA> &       data)
  {
    const unsigned int min = std::max(m.min_level(), data.min_level());
    const unsigned int max = std::min(m.max_level(), data.max_level());

    this->resize(min, max);

    for (unsigned int i = min; i <= max; ++i)
      (*this)[i].initialize(m[i], data[i]);
  }


  template <class RelaxationType, typename VectorType>
  inline void
  SmootherRelaxation<RelaxationType, VectorType>::smooth(
    const unsigned int level,
    VectorType &       u,
    const VectorType & rhs) const
  {
    unsigned int maxlevel = this->max_level();
    unsigned int steps2   = this->steps;

    if (this->variable)
      steps2 *= (1 << (maxlevel - level));

    bool T = this->transpose;
    if (this->symmetric && (steps2 % 2 == 0))
      T = false;
    if (this->debug > 0)
      deallog << 'S' << level << ' ';

    for (unsigned int i = 0; i < steps2; ++i)
      {
        if (T)
          (*this)[level].Tstep(u, rhs);
        else
          (*this)[level].step(u, rhs);
        if (this->symmetric)
          T = !T;
      }
  }


  template <class RelaxationType, typename VectorType>
  inline void
  SmootherRelaxation<RelaxationType, VectorType>::apply(
    const unsigned int level,
    VectorType &       u,
    const VectorType & rhs) const
  {
    unsigned int maxlevel = this->max_level();
    unsigned int steps2   = this->steps;

    if (this->variable)
      steps2 *= (1 << (maxlevel - level));

    bool T = this->transpose;
    if (this->symmetric && (steps2 % 2 == 0))
      T = false;
    if (this->debug > 0)
      deallog << 'S' << level << ' ';

    if (T)
      (*this)[level].Tvmult(u, rhs);
    else
      (*this)[level].vmult(u, rhs);
    if (this->symmetric)
      T = !T;
    for (unsigned int i = 1; i < steps2; ++i)
      {
        if (T)
          (*this)[level].Tstep(u, rhs);
        else
          (*this)[level].step(u, rhs);
        if (this->symmetric)
          T = !T;
      }
  }


  template <class RelaxationType, typename VectorType>
  inline std::size_t
  SmootherRelaxation<RelaxationType, VectorType>::memory_consumption() const
  {
    return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
           MGLevelObject<RelaxationType>::memory_consumption() +
           this->vector_memory.memory_consumption();
  }
} // namespace mg


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

template <typename MatrixType, class RelaxationType, typename VectorType>
inline MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
  MGSmootherRelaxation(const unsigned int steps,
                       const bool         variable,
                       const bool         symmetric,
                       const bool         transpose)
  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
{}



template <typename MatrixType, class RelaxationType, typename VectorType>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::clear()
{
  smoothers.clear_elements();

  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
  for (; i <= max_level; ++i)
    matrices[i] = LinearOperator<VectorType>();
}


template <typename MatrixType, class RelaxationType, typename VectorType>
template <typename MatrixType2>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &             m,
  const typename RelaxationType::AdditionalData &data)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      // Workaround: Unfortunately, not every "m[i]" object has a rich
      // enough interface to populate reinit_(domain|range)_vector. Thus,
      // apply an empty LinearOperator exemplar.
      matrices[i] =
        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
      smoothers[i].initialize(m[i], data);
    }
}

template <typename MatrixType, class RelaxationType, typename VectorType>
template <typename MatrixType2, class DATA>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &m,
  const MGLevelObject<DATA> &       data)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      // Workaround: Unfortunately, not every "m[i]" object has a rich
      // enough interface to populate reinit_(domain|range)_vector. Thus,
      // apply an empty LinearOperator exemplar.
      matrices[i] =
        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
      smoothers[i].initialize(m[i], data[i]);
    }
}

template <typename MatrixType, class RelaxationType, typename VectorType>
template <typename MatrixType2, class DATA>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &m,
  const DATA &                      data,
  const unsigned int                row,
  const unsigned int                col)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      // Workaround: Unfortunately, not every "m[i]" object has a rich
      // enough interface to populate reinit_(domain|range)_vector. Thus,
      // apply an empty LinearOperator exemplar.
      matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
                                                m[i].block(row, col));
      smoothers[i].initialize(m[i].block(row, col), data);
    }
}

template <typename MatrixType, class RelaxationType, typename VectorType>
template <typename MatrixType2, class DATA>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &m,
  const MGLevelObject<DATA> &       data,
  const unsigned int                row,
  const unsigned int                col)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      // Workaround: Unfortunately, not every "m[i]" object has a rich
      // enough interface to populate reinit_(domain|range)_vector. Thus,
      // apply an empty LinearOperator exemplar.
      matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
                                                m[i].block(row, col));
      smoothers[i].initialize(m[i].block(row, col), data[i]);
    }
}


template <typename MatrixType, class RelaxationType, typename VectorType>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::smooth(
  const unsigned int level,
  VectorType &       u,
  const VectorType & rhs) const
{
  unsigned int maxlevel = smoothers.max_level();
  unsigned int steps2   = this->steps;

  if (this->variable)
    steps2 *= (1 << (maxlevel - level));

  bool T = this->transpose;
  if (this->symmetric && (steps2 % 2 == 0))
    T = false;
  if (this->debug > 0)
    deallog << 'S' << level << ' ';

  for (unsigned int i = 0; i < steps2; ++i)
    {
      if (T)
        smoothers[level].Tstep(u, rhs);
      else
        smoothers[level].step(u, rhs);
      if (this->symmetric)
        T = !T;
    }
}


template <typename MatrixType, class RelaxationType, typename VectorType>
inline void
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::apply(
  const unsigned int level,
  VectorType &       u,
  const VectorType & rhs) const
{
  unsigned int maxlevel = smoothers.max_level();
  unsigned int steps2   = this->steps;

  if (this->variable)
    steps2 *= (1 << (maxlevel - level));

  bool T = this->transpose;
  if (this->symmetric && (steps2 % 2 == 0))
    T = false;
  if (this->debug > 0)
    deallog << 'S' << level << ' ';

  if (T)
    smoothers[level].Tvmult(u, rhs);
  else
    smoothers[level].vmult(u, rhs);
  if (this->symmetric)
    T = !T;
  for (unsigned int i = 1; i < steps2; ++i)
    {
      if (T)
        smoothers[level].Tstep(u, rhs);
      else
        smoothers[level].step(u, rhs);
      if (this->symmetric)
        T = !T;
    }
}



template <typename MatrixType, class RelaxationType, typename VectorType>
inline std::size_t
MGSmootherRelaxation<MatrixType, RelaxationType, VectorType>::
  memory_consumption() const
{
  return sizeof(*this) + matrices.memory_consumption() +
         smoothers.memory_consumption() +
         this->vector_memory.memory_consumption();
}


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

template <typename MatrixType, typename PreconditionerType, typename VectorType>
inline MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
  MGSmootherPrecondition(const unsigned int steps,
                         const bool         variable,
                         const bool         symmetric,
                         const bool         transpose)
  : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
{}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::clear()
{
  smoothers.clear_elements();

  unsigned int i = matrices.min_level(), max_level = matrices.max_level();
  for (; i <= max_level; ++i)
    matrices[i] = LinearOperator<VectorType>();
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
template <typename MatrixType2>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &                 m,
  const typename PreconditionerType::AdditionalData &data)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      // Workaround: Unfortunately, not every "m[i]" object has a rich
      // enough interface to populate reinit_(domain|range)_vector. Thus,
      // apply an empty LinearOperator exemplar.
      matrices[i] =
        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
      smoothers[i].initialize(m[i], data);
    }
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
template <typename MatrixType2, class DATA>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &m,
  const MGLevelObject<DATA> &       data)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      // Workaround: Unfortunately, not every "m[i]" object has a rich
      // enough interface to populate reinit_(domain|range)_vector. Thus,
      // apply an empty LinearOperator exemplar.
      matrices[i] =
        linear_operator<VectorType>(LinearOperator<VectorType>(), m[i]);
      smoothers[i].initialize(m[i], data[i]);
    }
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
template <typename MatrixType2, class DATA>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &m,
  const DATA &                      data,
  const unsigned int                row,
  const unsigned int                col)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      matrices[i] = &(m[i].block(row, col));
      smoothers[i].initialize(m[i].block(row, col), data);
    }
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
template <typename MatrixType2, class DATA>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::initialize(
  const MGLevelObject<MatrixType2> &m,
  const MGLevelObject<DATA> &       data,
  const unsigned int                row,
  const unsigned int                col)
{
  const unsigned int min = m.min_level();
  const unsigned int max = m.max_level();

  Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
  Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));

  matrices.resize(min, max);
  smoothers.resize(min, max);

  for (unsigned int i = min; i <= max; ++i)
    {
      matrices[i] = &(m[i].block(row, col));
      smoothers[i].initialize(m[i].block(row, col), data[i]);
    }
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::smooth(
  const unsigned int level,
  VectorType &       u,
  const VectorType & rhs) const
{
  unsigned int maxlevel = matrices.max_level();
  unsigned int steps2   = this->steps;

  if (this->variable)
    steps2 *= (1 << (maxlevel - level));

  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);

  r->reinit(u, true);
  d->reinit(u, true);

  bool T = this->transpose;
  if (this->symmetric && (steps2 % 2 == 0))
    T = false;
  if (this->debug > 0)
    deallog << 'S' << level << ' ';

  for (unsigned int i = 0; i < steps2; ++i)
    {
      if (T)
        {
          if (this->debug > 0)
            deallog << 'T';
          matrices[level].Tvmult(*r, u);
          r->sadd(-1., 1., rhs);
          if (this->debug > 2)
            deallog << ' ' << r->l2_norm() << ' ';
          smoothers[level].Tvmult(*d, *r);
          if (this->debug > 1)
            deallog << ' ' << d->l2_norm() << ' ';
        }
      else
        {
          if (this->debug > 0)
            deallog << 'N';
          matrices[level].vmult(*r, u);
          r->sadd(-1., rhs);
          if (this->debug > 2)
            deallog << ' ' << r->l2_norm() << ' ';
          smoothers[level].vmult(*d, *r);
          if (this->debug > 1)
            deallog << ' ' << d->l2_norm() << ' ';
        }
      u += *d;
      if (this->symmetric)
        T = !T;
    }
  if (this->debug > 0)
    deallog << std::endl;
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
inline void
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::apply(
  const unsigned int level,
  VectorType &       u,
  const VectorType & rhs) const
{
  unsigned int maxlevel = matrices.max_level();
  unsigned int steps2   = this->steps;

  if (this->variable)
    steps2 *= (1 << (maxlevel - level));

  bool T = this->transpose;
  if (this->symmetric && (steps2 % 2 == 0))
    T = false;
  if (this->debug > 0)
    deallog << 'S' << level << ' ';

  // first step where we overwrite the result
  if (this->debug > 2)
    deallog << ' ' << rhs.l2_norm() << ' ';
  if (this->debug > 0)
    deallog << (T ? 'T' : 'N');
  if (T)
    smoothers[level].Tvmult(u, rhs);
  else
    smoothers[level].vmult(u, rhs);
  if (this->debug > 1)
    deallog << ' ' << u.l2_norm() << ' ';
  if (this->symmetric)
    T = !T;

  typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
  typename VectorMemory<VectorType>::Pointer d(this->vector_memory);

  if (steps2 > 1)
    {
      r->reinit(u, true);
      d->reinit(u, true);
    }

  for (unsigned int i = 1; i < steps2; ++i)
    {
      if (T)
        {
          if (this->debug > 0)
            deallog << 'T';
          matrices[level].Tvmult(*r, u);
          r->sadd(-1., 1., rhs);
          if (this->debug > 2)
            deallog << ' ' << r->l2_norm() << ' ';
          smoothers[level].Tvmult(*d, *r);
          if (this->debug > 1)
            deallog << ' ' << d->l2_norm() << ' ';
        }
      else
        {
          if (this->debug > 0)
            deallog << 'N';
          matrices[level].vmult(*r, u);
          r->sadd(-1., rhs);
          if (this->debug > 2)
            deallog << ' ' << r->l2_norm() << ' ';
          smoothers[level].vmult(*d, *r);
          if (this->debug > 1)
            deallog << ' ' << d->l2_norm() << ' ';
        }
      u += *d;
      if (this->symmetric)
        T = !T;
    }
  if (this->debug > 0)
    deallog << std::endl;
}



template <typename MatrixType, typename PreconditionerType, typename VectorType>
inline std::size_t
MGSmootherPrecondition<MatrixType, PreconditionerType, VectorType>::
  memory_consumption() const
{
  return sizeof(*this) + matrices.memory_consumption() +
         smoothers.memory_consumption() +
         this->vector_memory.memory_consumption();
}


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif