mg_smoother.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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 ();

    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;

    std::size_t memory_consumption () const;
  };
}

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 ();

  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<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();
  }
}


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

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