multigrid.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/solver/base.hpp>
 
// includes, system
#include <deque>
#include <vector>
 
namespace FEAT
{
  namespace Solver
  {
    enum class MultiGridCycle
    {
      V,
      F,
      W
    };
 
    enum class MultiGridAdaptCGC
    {
      Fixed,
      MinEnergy,
      MinDefect
    };
 
    inline std::ostream& operator<<(std::ostream& os, MultiGridCycle cycle)
    {
      switch(cycle)
      {
      case MultiGridCycle::V:
        return os << "V";
      case MultiGridCycle::F:
        return os << "F";
      case MultiGridCycle::W:
        return os << "W";
      default:
        return os << "?";
      }
    }
 
    inline std::istream& operator>>(std::istream& is, MultiGridCycle& cycle)
    {
      // read character
      char c(' ');
      if((is >> c).fail())
        return is;
 
      switch(c)
      {
      case 'v':
      case 'V':
        cycle = MultiGridCycle::V;
        break;
 
      case 'f':
      case 'F':
        cycle = MultiGridCycle::F;
        break;
 
      case 'w':
      case 'W':
        cycle = MultiGridCycle::W;
        break;
 
      default:
        // invalid character
        is.putback(c);
        is.setstate(std::ios_base::failbit);
        break;
      }
 
      return is;
    }
 
    template<
      typename SystemMatrix_,
      typename SystemFilter_,
      typename TransferOperator_>
    class MultiGridLevelBase
    {
    public:
      typedef SystemMatrix_ SystemMatrixType;
      typedef SystemFilter_ SystemFilterType;
      typedef TransferOperator_ TransferOperatorType;
      typedef typename SystemMatrix_::VectorTypeR SystemVectorType;
      typedef SolverBase<SystemVectorType> SolverType;
 
    public:
      virtual ~MultiGridLevelBase()
      {
      }
 
      virtual const SystemMatrixType& get_system_matrix() const = 0;
 
      virtual const SystemFilterType& get_system_filter() const = 0;
 
      virtual const TransferOperatorType* get_transfer_operator() const = 0;
 
      virtual std::shared_ptr<SolverType> get_coarse_solver() = 0;
 
      virtual std::shared_ptr<SolverType> get_smoother_pre() = 0;
 
      virtual std::shared_ptr<SolverType> get_smoother_post() = 0;
 
      virtual std::shared_ptr<SolverType> get_smoother_peak() = 0;
    }; // class MultiGridLevelBase<...>
 
    template<
      typename SystemMatrix_,
      typename SystemFilter_,
      typename TransferOperator_>
    class MultiGridLevelStd :
      public MultiGridLevelBase<SystemMatrix_, SystemFilter_, TransferOperator_>
    {
    public:
      typedef MultiGridLevelBase<SystemMatrix_, SystemFilter_, TransferOperator_> BaseClass;
      typedef typename BaseClass::SystemMatrixType SystemMatrixType;
      typedef typename BaseClass::SystemFilterType SystemFilterType;
      typedef typename BaseClass::SystemVectorType SystemVectorType;
      typedef typename BaseClass::TransferOperatorType TransferOperatorType;
      typedef typename BaseClass::SolverType SolverType;
 
    protected:
      const SystemMatrixType& system_matrix;
      const SystemFilterType& system_filter;
      const TransferOperatorType* transfer_operator;
      std::shared_ptr<SolverType> coarse_solver;
      std::shared_ptr<SolverType> smoother_pre;
      std::shared_ptr<SolverType> smoother_post;
      std::shared_ptr<SolverType> smoother_peak;
 
    public:
      explicit MultiGridLevelStd(
        const SystemMatrixType& sys_matrix,
        const SystemFilterType& sys_filter,
        std::shared_ptr<SolverType> crs_solver)
         :
        system_matrix(sys_matrix),
        system_filter(sys_filter),
        transfer_operator(nullptr),
        coarse_solver(crs_solver),
        smoother_pre(),
        smoother_post(),
        smoother_peak()
      {
      }
 
      explicit MultiGridLevelStd(
        const SystemMatrixType& sys_matrix,
        const SystemFilterType& sys_filter,
        const TransferOperatorType& trans_operat,
        std::shared_ptr<SolverType> smooth_pre,
        std::shared_ptr<SolverType> smooth_post,
        std::shared_ptr<SolverType> smooth_peak,
        std::shared_ptr<SolverType> crs_solver = nullptr)
         :
        system_matrix(sys_matrix),
        system_filter(sys_filter),
        transfer_operator(&trans_operat),
        coarse_solver(crs_solver),
        smoother_pre(smooth_pre),
        smoother_post(smooth_post),
        smoother_peak(smooth_peak)
      {
      }
 
      virtual ~MultiGridLevelStd()
      {
      }
 
      virtual const SystemMatrixType& get_system_matrix() const override
      {
        return system_matrix;
      }
 
      virtual const SystemFilterType& get_system_filter() const override
      {
        return system_filter;
      }
 
      virtual const TransferOperatorType* get_transfer_operator() const override
      {
        return transfer_operator;
      }
 
      virtual std::shared_ptr<SolverType> get_coarse_solver() override
      {
        return coarse_solver;
      }
 
      virtual std::shared_ptr<SolverType> get_smoother_pre() override
      {
        return smoother_pre;
      }
 
      virtual std::shared_ptr<SolverType> get_smoother_post() override
      {
        return smoother_post;
      }
 
      virtual std::shared_ptr<SolverType> get_smoother_peak() override
      {
        return smoother_peak;
      }
    }; // class MultiGridLevelStd<...>
 
    // forward declaration
    template<
      typename SystemMatrix_,
      typename SystemFilter_,
      typename TransferOperator_>
    class MultiGrid;
 
    template<
      typename SystemMatrix_,
      typename SystemFilter_,
      typename TransferOperator_>
    class MultiGridHierarchy
    {
    public:
      friend class MultiGrid<SystemMatrix_, SystemFilter_, TransferOperator_>;
 
      typedef MultiGridLevelBase<SystemMatrix_, SystemFilter_, TransferOperator_> LevelType;
      typedef MultiGridLevelStd<SystemMatrix_, SystemFilter_, TransferOperator_> StdLevelType;
 
      typedef typename LevelType::SystemMatrixType SystemMatrixType;
      typedef typename LevelType::SystemFilterType SystemFilterType;
      typedef typename LevelType::SystemVectorType SystemVectorType;
      typedef typename LevelType::TransferOperatorType TransferOperatorType;
      typedef typename LevelType::SolverType SolverType;
 
    protected:
      class LevelInfo
      {
      public:
        std::shared_ptr<LevelType> level;
        std::deque<SolverType*> unique_solvers;
        SystemVectorType vec_rhs, vec_sol, vec_def, vec_cor, vec_tmp;
 
        double time_smooth;
 
        double time_coarse;
 
        double time_defect;
 
        double time_transfer;
 
        explicit LevelInfo(std::shared_ptr<LevelType> lvl) :
          level(lvl),
          unique_solvers()
        {
          reset_timings();
        }
 
        void reset_timings()
        {
          time_smooth = time_coarse = time_defect = time_transfer = 0.0;
        }
 
        void init_symbolic()
        {
          // build unique solver list
          unique_solvers.clear();
          _push_solver(level->get_coarse_solver().get());
          _push_solver(level->get_smoother_pre().get());
          _push_solver(level->get_smoother_post().get());
          _push_solver(level->get_smoother_peak().get());
 
          // initialize all unique solvers in forward order
          for(auto it = unique_solvers.begin(); it != unique_solvers.end(); ++it)
            (*it)->init_symbolic();
 
          // get level system matrix
          const SystemMatrixType& matrix = level->get_system_matrix();
 
          // create our vectors
          vec_rhs = matrix.create_vector_r();
          vec_sol = matrix.create_vector_r();
          vec_def = matrix.create_vector_r();
          vec_cor = matrix.create_vector_r();
          vec_tmp = matrix.create_vector_r();
        }
 
        void done_symbolic()
        {
          // release our vectors
          vec_tmp.clear();
          vec_cor.clear();
          vec_def.clear();
          vec_sol.clear();
          vec_rhs.clear();
 
          // release all unique solvers in reverse order
          for(auto it = unique_solvers.rbegin(); it != unique_solvers.rend(); ++it)
            (*it)->done_symbolic();
        }
 
        void init_numeric()
        {
          // initialize all unique solvers in forward order
          for(auto it = unique_solvers.begin(); it != unique_solvers.end(); ++it)
            (*it)->init_numeric();
        }
 
        void done_numeric()
        {
          // release all unique solvers in reverse order
          for(auto it = unique_solvers.rbegin(); it != unique_solvers.rend(); ++it)
            (*it)->done_numeric();
        }
 
      protected:
        void _push_solver(SolverType* solver)
        {
          if(solver == nullptr)
            return;
 
          // loop over all existing solvers to determine uniqueness
          for(auto it = unique_solvers.begin(); it != unique_solvers.end(); ++it)
          {
            // do we already have this solver object?
            if(*it == solver)
              return;
          }
 
          // alright, that's a new unique solver object, so add it
          unique_solvers.push_back(solver);
        }
      }; // class LevelInfo
 
      LevelInfo& _get_level_info(Index lvl)
      {
        XASSERTM(lvl < size_physical(), "invalid level index");
        return _levels.at(std::size_t(lvl));
      }
 
    protected:
      std::deque<LevelInfo> _levels;
      std::size_t _virt_size;
      bool _have_init_symbolic;
      bool _have_init_numeric;
 
    public:
      explicit MultiGridHierarchy(std::size_t size_virt) :
        _levels(),
        _virt_size(size_virt),
        _have_init_symbolic(false),
        _have_init_numeric(false)
      {
      }
 
      virtual ~MultiGridHierarchy()
      {
        // delete levels in reverse order
        while(!_levels.empty())
          _levels.pop_back();
      }
 
      void push_level(std::shared_ptr<LevelType> level)
      {
        XASSERTM(!_have_init_symbolic, "cannot push new level, init_symbolic() was already called");
        XASSERTM(!_have_init_numeric,  "cannot push new level, init_numeric() was already called");
        XASSERTM(_levels.size() < _virt_size, "cannot push new level, already have maximum number of levels");
        _levels.push_back(LevelInfo(level));
      }
 
      void push_level(
        const SystemMatrixType& system_matrix,
        const SystemFilterType& system_filter,
        std::shared_ptr<SolverType> coarse_solver)
      {
        push_level(std::make_shared<StdLevelType>(system_matrix, system_filter, coarse_solver));
      }
 
      void push_level(
        const SystemMatrixType& system_matrix,
        const SystemFilterType& system_filter,
        const TransferOperatorType& transfer_operator,
        std::shared_ptr<SolverType> smoother_pre,
        std::shared_ptr<SolverType> smoother_post,
        std::shared_ptr<SolverType> smoother_peak,
        std::shared_ptr<SolverType> coarse_solver = nullptr
        )
      {
        push_level(std::make_shared<StdLevelType>(system_matrix, system_filter, transfer_operator,
          smoother_pre, smoother_post, smoother_peak, coarse_solver));
      }
 
      Index size_physical() const
      {
        return Index(_levels.size());
      }
 
      Index size_virtual() const
      {
        return Index(_virt_size);
      }
 
      LevelType& get_level(Index lvl)
      {
        XASSERTM(lvl < size_physical(), "invalid level index");
        return *(_levels.at(std::size_t(lvl)).level);
      }
 
      const LevelType& get_level(Index lvl) const
      {
        XASSERTM(lvl < size_physical(), "invalid level index");
        return *(_levels.at(std::size_t(lvl)).level);
      }
 
      void init_symbolic()
      {
        XASSERTM(!_have_init_numeric, "cannot call init_symbolic(): init_numeric() has already been called");
        XASSERTM(!_have_init_symbolic, "cannot call init_symbolic(): init_symbolic() has already been called");
 
        // initialize all levels in forward order
        for(auto it = _levels.begin(); it != _levels.end(); ++it)
          it->init_symbolic();
 
        _have_init_symbolic = true;
      }
 
      void done_symbolic()
      {
        XASSERTM(!_have_init_numeric, "cannot call done_symbolic(): done_numeric() has not been called yet");
        XASSERTM(_have_init_symbolic, "cannot call done_symbolic(): init_symbolic() has not been called yet");
 
        // release all levels in reverse order
        for(auto it = _levels.rbegin(); it != _levels.rend(); ++it)
          it->done_symbolic();
 
        _have_init_symbolic = false;
      }
 
      void init_numeric()
      {
        XASSERTM(!_have_init_numeric, "cannot call init_numeric(): init_numeric() has already been called");
        XASSERTM(_have_init_symbolic, "cannot call init_numeric(): init_symbolic() has not been called yet");
 
        // initialize all levels in forward order
        for(auto it = _levels.begin(); it != _levels.end(); ++it)
          it->init_numeric();
 
        _have_init_numeric = true;
      }
 
      void done_numeric()
      {
        XASSERTM(_have_init_numeric, "cannot call done_numeric(): init_numeric() has not been called yet");
 
        // release all levels in reverse order
        for(auto it = _levels.rbegin(); it != _levels.rend(); ++it)
          it->done_numeric();
 
        _have_init_numeric = false;
      }
 
      void init()
      {
        init_symbolic();
        init_numeric();
      }
 
      void done()
      {
        done_numeric();
        done_symbolic();
      }
 
      bool have_symbolic() const
      {
        return _have_init_symbolic;
      }
 
      bool have_numeric() const
      {
        return _have_init_numeric;
      }
 
      void reset_timings()
      {
        for(auto& li : _levels)
          li.reset_timings();
      }
 
      double get_time_smooth(int lvl = -1) const
      {
        if(lvl < 0)
        {
          double t(0.0);
          for(const auto& li : _levels)
            t += li.time_smooth;
          return t;
        }
 
        XASSERT(lvl < int(_levels.size()));
        return _levels.at(std::size_t(lvl)).time_smooth;
      }
 
      double get_time_coarse(int lvl = -1) const
      {
        if(lvl < 0)
        {
          double t(0.0);
          for(const auto& li : _levels)
            t += li.time_coarse;
          return t;
        }
 
        XASSERT(lvl < int(_levels.size()));
        return _levels.at(std::size_t(lvl)).time_coarse;
      }
 
      double get_time_defect(int lvl = -1) const
      {
        if(lvl < 0)
        {
          double t(0.0);
          for(const auto& li : _levels)
            t += li.time_defect;
          return t;
        }
 
        XASSERT(lvl < int(_levels.size()));
        return _levels.at(std::size_t(lvl)).time_defect;
      }
 
      double get_time_transfer(int lvl = -1) const
      {
        if(lvl < 0)
        {
          double t(0.0);
          for(const auto& li : _levels)
            t += li.time_transfer;
          return t;
        }
 
        XASSERT(lvl < int(_levels.size()));
        return _levels.at(std::size_t(lvl)).time_transfer;
      }
    }; // class MultiGridHierarchy<...>
 
    template<
      typename SystemMatrix_,
      typename SystemFilter_,
      typename TransferOperator_>
    class MultiGrid :
      public SolverBase<typename SystemMatrix_::VectorTypeR>
    {
    public:
      typedef SolverBase<typename SystemMatrix_::VectorTypeR> BaseClass;
 
      typedef MultiGridHierarchy<SystemMatrix_, SystemFilter_, TransferOperator_> HierarchyType;
 
      typedef typename HierarchyType::LevelType LevelType;
      typedef typename HierarchyType::LevelInfo LevelInfo;
 
      typedef typename LevelType::SystemMatrixType MatrixType;
      typedef typename LevelType::SystemFilterType FilterType;
      typedef typename LevelType::SystemVectorType VectorType;
      typedef typename LevelType::TransferOperatorType TransferOperatorType;
      typedef typename LevelType::SolverType SolverType;
 
      typedef typename MatrixType::DataType DataType;
 
    protected:
      std::shared_ptr<HierarchyType> _hierarchy;
      MultiGridCycle _cycle;
      MultiGridAdaptCGC _adapt_cgc;
      Index _top_level;
      Index _crs_level;
 
      std::vector<int> _counters;
 
      std::vector<double> _toes;
      std::vector<double> _mpi_execs_reduction;
      std::vector<double> _mpi_execs_blas2;
      std::vector<double> _mpi_execs_blas3;
      std::vector<double> _mpi_execs_collective;
      std::vector<double> _mpi_waits_reduction;
      std::vector<double> _mpi_waits_blas2;
      std::vector<double> _mpi_waits_blas3;
      std::vector<double> _mpi_waits_collective;
 
    public:
      explicit MultiGrid(std::shared_ptr<HierarchyType> hierarchy, MultiGridCycle cycle, int top_level = 0, int crs_level = -1) :
        BaseClass(),
        _hierarchy(hierarchy),
        _cycle(cycle),
        _adapt_cgc(MultiGridAdaptCGC::Fixed),
        _top_level(Index(top_level)),
        _crs_level(crs_level >= 0 ? Index(crs_level) : Index(int(_hierarchy->size_virtual()) + crs_level))
      {
        XASSERTM(_crs_level < _hierarchy->size_virtual(), "invalid coarse level");
        XASSERTM(_top_level <= _crs_level, "invalid top/coarse level combination");
      }
 
      virtual ~MultiGrid()
      {
      }
 
      void set_adapt_cgc(MultiGridAdaptCGC adapt_cgc)
      {
        _adapt_cgc = adapt_cgc;
      }
 
      void set_cycle(MultiGridCycle cycle)
      {
        _cycle = cycle;
      }
 
      MultiGridCycle get_cycle() const
      {
        return _cycle;
      }
 
      void set_levels(int top_level, int crs_level)
      {
        _top_level = Index(top_level);
        _crs_level = (crs_level >= 0 ? Index(crs_level) : Index(int(_hierarchy->size_virtual()) + crs_level));
 
        XASSERTM(_crs_level < _hierarchy->size_virtual(), "invalid coarse level");
        XASSERTM(_top_level <= _crs_level, "invalid top/coarse level combination");
      }
 
      Index get_top_level() const
      {
        return _top_level;
      }
 
      Index get_crs_level() const
      {
        return _crs_level;
      }
 
      virtual String name() const override
      {
        return "MultiGrid-" + stringify(_cycle);
      }
 
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        // ensure that the hierarchy was initialized
        XASSERTM(_hierarchy->have_symbolic(), "init_symbolic() of multigrid hierarchy was not called yet");
 
        // get the virtual hierarchy size
        const std::size_t num_lvls = _hierarchy->size_virtual();
 
        // allocate counter vector
        _counters.resize(num_lvls, 0);
 
        // allocate statistics vectors
        _toes.resize(num_lvls, 0.0);
        _mpi_execs_reduction.resize(num_lvls, 0.0);
        _mpi_execs_blas2.resize(num_lvls, 0.0);
        _mpi_execs_blas3.resize(num_lvls, 0.0);
        _mpi_execs_collective.resize(num_lvls, 0.0);
        _mpi_waits_reduction.resize(num_lvls, 0.0);
        _mpi_waits_blas2.resize(num_lvls, 0.0);
        _mpi_waits_blas3.resize(num_lvls, 0.0);
        _mpi_waits_collective.resize(num_lvls, 0.0);
      }
 
      virtual void done_symbolic() override
      {
        _toes.clear();
        _mpi_execs_reduction.clear();
        _mpi_execs_blas2.clear();
        _mpi_execs_blas3.clear();
        _mpi_execs_collective.clear();
        _mpi_waits_reduction.clear();
        _mpi_waits_blas2.clear();
        _mpi_waits_blas3.clear();
        _mpi_waits_collective.clear();
        _counters.clear();
 
        BaseClass::done_symbolic();
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        // ensure that the hierarchy was initialized
        XASSERTM(_hierarchy->have_numeric(), "init_numeric() of multigrid hierarchy was not called yet");
      }
 
      virtual void done_numeric() override
      {
        BaseClass::done_numeric();
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
 
        // reset statistics counters
        for(std::size_t i(0); i <  std::size_t(_hierarchy->size_virtual()); ++i)
        {
          _toes.at(i) = 0.0;
          _mpi_execs_reduction.at(i) = 0.0;
          _mpi_execs_blas2.at(i) = 0.0;
          _mpi_execs_blas3.at(i) = 0.0;
          _mpi_execs_collective.at(i) = 0.0;
          _mpi_waits_reduction.at(i) = 0.0;
          _mpi_waits_blas2.at(i) = 0.0;
          _mpi_waits_blas3.at(i) = 0.0;
          _mpi_waits_collective.at(i) = 0.0;
        }
 
        // get a reference to the finest level
        LevelInfo& lvl_top = _hierarchy->_get_level_info(_top_level);
 
        // copy defect to RHS vector
        lvl_top.vec_rhs.copy(vec_def);
 
        // apply the corresponding cycle
        Status status(Status::undefined);
        switch(_cycle)
        {
        case MultiGridCycle::V:
          status = this->_apply_cycle_v();
          break;
 
        case MultiGridCycle::F:
          status = this->_apply_cycle_f();
          break;
 
        case MultiGridCycle::W:
          status = this->_apply_cycle_w();
          break;
 
        default:
          // whoops...
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::aborted, 1));
          status = Status::aborted;
          break;
        }
 
        // copy solution to correction vector
        vec_cor.copy(lvl_top.vec_sol);
 
        // propagate solver statistics
        for(std::size_t i(0); i <  std::size_t(_hierarchy->size_virtual()); ++i)
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionLevelTimings>(this->name(), Index(i),
            _toes.at(i), _mpi_execs_reduction.at(i), _mpi_execs_blas2.at(i), _mpi_execs_blas3.at(i), _mpi_execs_collective.at(i), _mpi_waits_reduction.at(i), _mpi_waits_blas2.at(i),
            _mpi_waits_blas3.at(i), _mpi_waits_collective.at(i)));
        }
 
        // okay
        Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::success, 1));
        return status;
      }
 
    protected:
      Status _apply_cycle_v()
      {
        // The V-cycle is pretty simple:
        // 1) Restrict from the top level to the coarse level
        // 2) Solve coarse level system
        // 3) Prolongate from the coarse level to the top level
 
        // restrict (apply top-level pre-smoother if it exists)
        this->_apply_rest(_top_level, true);
 
        // solve coarse level
        this->_apply_coarse();
 
        // prolongate (apply top-level post-smoother if it exists)
        this->_apply_prol(_top_level, true);
 
        // okay
        return Status::success;
      }
 
      Status _apply_cycle_f()
      {
        // The F-cycle is slightly more complex than the V-cycle:
        // 1) Restrict from the top level to the coarse level
        //
        // 2) For each intermediate level (that's everything between the
        //    coarse and top level) in ascending order:
        //    2.1) solve the coarse system
        //    2.2) prolongate to the intermediate level without applying
        //         the post-smoother on the intermediate level
        //    2.3) apply peak-smoother
        //    2.4) restrict to the coarse level without applying
        //         the pre-smoother on the intermediate level
        //
        // 3) Solve the coarse level system
        //
        // 4) Prolongate from the coarse level to the top level
 
        // restrict (apply top-level pre-smoother if it exists)
        this->_apply_rest(_top_level, true);
 
        // get the coarsest level on this process; that's either the
        // coarse level or the last level for this process:
        const Index last_level = Math::min(_crs_level, _hierarchy->size_physical());
 
        // ensure that we do not apply this if there is just one level
        if(last_level > Index(0))
        {
          // F-cycle intermediate peak-level loop
          for(Index peak_lvl(last_level-1); peak_lvl > _top_level; --peak_lvl)
          {
            // solve coarse level
            this->_apply_coarse();
 
            // prolongate to current peak-level without post-smoothing
            this->_apply_prol(peak_lvl, false);
 
            // apply peak-smoother
            this->_apply_smooth_peak(peak_lvl);
 
            // restrict from current peak-level without pre-smoothing
            this->_apply_rest(peak_lvl, false);
          }
        }
 
        // solve coarse level
        this->_apply_coarse();
 
        // prolongate (apply top-level post-smoother if it exists)
        this->_apply_prol(_top_level, true);
 
        // okay
        return Status::success;
      }
 
      Status _apply_cycle_w()
      {
        // The W-cycle is a tough nut to implement/understand:
        // As always, we start by restricting from the top level to the
        // coarse level and solving the coarse level system afterwards.
        //
        // The next part is the tricky bit:
        // We need to implement the "inner W-cycle part". In theory,
        // this is usually defined as a "recursive 2x V-cycle application",
        // but we do not want to use recursion here explicitly.
        // We know that the W-cycle performs 2^(top_level - crs_level)
        // coarse grid solves, so we use that for our for-loop.
        // Moreover, we manage a counter for each level above the coarse
        // level and count how often we have applied the peak-smoother.
        // Each iteration we choose the lowest level that has a peak-smooth
        // count equal to 0 as the next intermediate level, increment
        // its counter by 1 and reset all level counters below this level
        // back to 0. This may sound perfectly confusing, but it actually
        // does the job :) See the description of the multigrid cycles
        // in the documentation for further explanations and a fancy image
        // of the W-cycle level counters in action.
        //
        // Oh yes, and finally we also need to prolongate from the
        // coarse level to the top level, of course.
 
        // restrict (apply top-level pre-smoother if it exists)
        this->_apply_rest(_top_level, true);
 
        // solve coarse level
        this->_apply_coarse();
 
        // get the coarsest level on this process; that's either the
        // coarse level or the last level for this process:
        const Index last_level = Math::min(_crs_level, _hierarchy->size_physical());
 
        // reset all level counters to 0
        for(std::size_t i(_top_level); i <= std::size_t(last_level); ++i)
        {
          _counters[i] = 0;
        }
 
        // compute the total number of coarse-grid solves to perform:
        // that's = 2^(top_level - crs_level)
        const int num_cgs = 1 << (last_level - _top_level);
 
        // W-cycle loop: count the number of coarse-grid solves
        for(int cgs(1); cgs < num_cgs; ++cgs)
        {
          // find the next peak level; that's the lowest level
          // above the coarse level whose counter is 0
          std::size_t peak_lvl = std::size_t(last_level);
          while(peak_lvl > std::size_t(_top_level))
          {
            // if the counter is 0, we have our next peak level
            if(_counters[--peak_lvl] == 0)
              break;
          }
 
          // >>>>> DEBUG >>>>>
          //std::cout << stringify(cgs).pad_front(3) << " ";
          //std::cout << stringify(peak_lvl).pad_front(2) << ":";
          //for(std::size_t i(0); i <= std::size_t(top_lvl); ++i)
          //  std::cout << " " << _counters[i];
          // <<<<< DEBUG <<<<<
 
          // sanity check: do not go beyond our top level
          XASSERTM(peak_lvl >= _top_level, "W-cycle sanity check failed: invalid peak level");
 
          // reset all counters below our peak level to 0
          for(std::size_t i = std::size_t(last_level-1); i > peak_lvl; --i)
            _counters[i] = 0;
 
          // increment the counter of our peak level by 1,
          // indicating that we are going to perform an
          // peak-smoother step on this level
          ++_counters[peak_lvl];
 
          // >>>>> DEBUG >>>>>
          //std::cout << " >";
          //for(std::size_t i(0); i <= std::size_t(top_lvl); ++i)
          //  std::cout << " " << _counters[i];
          //std::cout << "\n";
          // <<<<< DEBUG <<<<<
 
          // prolongate to current peak-level without post-smoothing
          this->_apply_prol(Index(peak_lvl), false);
 
          // apply peak-smoother
          this->_apply_smooth_peak(Index(peak_lvl));
 
          // restrict from current peak-level without pre-smoothing
          this->_apply_rest(Index(peak_lvl), false);
 
          // solve coarse level
          this->_apply_coarse();
        }
 
        // sanity check: all level counters above the coarse level must be 1
        for(std::size_t i(last_level); i > std::size_t(_top_level); )
        {
          --i;
          XASSERTM(this->_counters[i] == 1, "W-cycle sanity check failed: invalid counters");
        }
 
        // prolongate (apply top-level post-smoother if it exists)
        this->_apply_prol(_top_level, true);
 
        // okay
        return Status::success;
      }
 
      Status _apply_coarse()
      {
        // no coarse level on this process?
        if(_crs_level >= _hierarchy->size_physical())
          return Status::success;
 
        // process the coarse grid level
        TimeStamp at;
        double mpi_exec_reduction_start(Statistics::get_time_mpi_execute_reduction());
        double mpi_exec_blas2_start(Statistics::get_time_mpi_execute_blas2());
        double mpi_exec_blas3_start(Statistics::get_time_mpi_execute_blas3());
        double mpi_exec_collective_start(Statistics::get_time_mpi_execute_collective());
        double mpi_wait_start_reduction(Statistics::get_time_mpi_wait_reduction());
        double mpi_wait_start_blas2(Statistics::get_time_mpi_wait_blas2());
        double mpi_wait_start_blas3(Statistics::get_time_mpi_wait_blas3());
        double mpi_wait_start_collective(Statistics::get_time_mpi_wait_collective());
 
        // get the coarse level info
        LevelInfo& lvl_crs = _hierarchy->_get_level_info(_crs_level);
 
        // get the system filter
        const FilterType& system_filter = lvl_crs.level->get_system_filter();
 
        // get the coarse grid solver
        std::shared_ptr<SolverType> coarse_solver = lvl_crs.level->get_coarse_solver();
 
        // if the have a coarse grid solver, apply it
        TimeStamp stamp_coarse;
        if(coarse_solver)
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionCallCoarseSolver>(this->name(), coarse_solver->name()));
          if(!status_success(coarse_solver->apply(lvl_crs.vec_sol, lvl_crs.vec_rhs)))
            return Status::aborted;
        }
        else
        {
          // simply copy the RHS thus emulating an identity solver
          lvl_crs.vec_sol.copy(lvl_crs.vec_rhs);
 
          // apply the correction filter
          system_filter.filter_cor(lvl_crs.vec_sol);
        }
        lvl_crs.time_coarse += stamp_coarse.elapsed_now();
 
        _toes.at(std::size_t(_crs_level)) += at.elapsed_now();
        double mpi_exec_reduction_stop(Statistics::get_time_mpi_execute_reduction());
        double mpi_exec_blas2_stop(Statistics::get_time_mpi_execute_blas2());
        double mpi_exec_blas3_stop(Statistics::get_time_mpi_execute_blas3());
        double mpi_exec_collective_stop(Statistics::get_time_mpi_execute_collective());
        double mpi_wait_stop_reduction(Statistics::get_time_mpi_wait_reduction());
        double mpi_wait_stop_blas2(Statistics::get_time_mpi_wait_blas2());
        double mpi_wait_stop_blas3(Statistics::get_time_mpi_wait_blas3());
        double mpi_wait_stop_collective(Statistics::get_time_mpi_wait_collective());
        _mpi_execs_reduction.at(std::size_t(_crs_level)) += mpi_exec_reduction_stop - mpi_exec_reduction_start;
        _mpi_execs_blas2.at(std::size_t(_crs_level)) += mpi_exec_blas2_stop - mpi_exec_blas2_start;
        _mpi_execs_blas3.at(std::size_t(_crs_level)) += mpi_exec_blas3_stop - mpi_exec_blas3_start;
        _mpi_execs_collective.at(std::size_t(_crs_level)) += mpi_exec_collective_stop - mpi_exec_collective_start;
        _mpi_waits_reduction.at(std::size_t(_crs_level)) += mpi_wait_stop_reduction - mpi_wait_start_reduction;
        _mpi_waits_blas2.at(std::size_t(_crs_level)) += mpi_wait_stop_blas2 - mpi_wait_start_blas2;
        _mpi_waits_blas3.at(std::size_t(_crs_level)) += mpi_wait_stop_blas3 - mpi_wait_start_blas3;
        _mpi_waits_collective.at(std::size_t(_crs_level)) += mpi_wait_stop_collective - mpi_wait_start_collective;
 
        return Status::success;
      }
 
      bool _apply_smooth_def(Index cur_lvl, SolverType& smoother)
      {
        // get the level info
        LevelInfo& lvl = _hierarchy->_get_level_info(cur_lvl);
 
        // get the system matrix and filter
        const MatrixType& system_matrix = lvl.level->get_system_matrix();
        const FilterType& system_filter = lvl.level->get_system_filter();
 
        // apply peak-smoother
        TimeStamp stamp_smooth;
        Statistics::add_solver_expression(std::make_shared<ExpressionCallSmoother>(this->name(), smoother.name()));
        smoother.apply(lvl.vec_cor, lvl.vec_def);
        //if(!status_success(smoother.apply(lvl.vec_cor, lvl.vec_def)))
          //return false;
        lvl.time_smooth += stamp_smooth.elapsed_now();
 
        // apply correction filter
        system_filter.filter_cor(lvl.vec_cor);
 
        // update solution vector
        lvl.vec_sol.axpy(lvl.vec_cor);
 
        // re-compute defect
        TimeStamp stamp_defect;
        system_matrix.apply(lvl.vec_def, lvl.vec_sol, lvl.vec_rhs, -DataType(1));
        lvl.time_defect += stamp_defect.elapsed_now();
 
        // apply defect filter
        system_filter.filter_def(lvl.vec_def);
 
        // okay
        return true;
      }
 
      Status _apply_smooth_peak(Index cur_lvl)
      {
        TimeStamp at;
        double mpi_exec_reduction_start(Statistics::get_time_mpi_execute_reduction());
        double mpi_exec_blas2_start(Statistics::get_time_mpi_execute_blas2());
        double mpi_exec_blas3_start(Statistics::get_time_mpi_execute_blas3());
        double mpi_exec_collective_start(Statistics::get_time_mpi_execute_collective());
        double mpi_wait_start_reduction(Statistics::get_time_mpi_wait_reduction());
        double mpi_wait_start_blas2(Statistics::get_time_mpi_wait_blas2());
        double mpi_wait_start_blas3(Statistics::get_time_mpi_wait_blas3());
        double mpi_wait_start_collective(Statistics::get_time_mpi_wait_collective());
 
        // get the level info
        LevelInfo& lvl = _hierarchy->_get_level_info(cur_lvl);
 
        // get the system matrix and filter
        const MatrixType& system_matrix = lvl.level->get_system_matrix();
        const FilterType& system_filter = lvl.level->get_system_filter();
 
        // compute defect
        TimeStamp stamp_defect;
        system_matrix.apply(lvl.vec_def, lvl.vec_sol, lvl.vec_rhs, -DataType(1));
        lvl.time_defect += stamp_defect.elapsed_now();
 
        // apply defect filter
        system_filter.filter_def(lvl.vec_def);
 
        // get the peak-smoother
        std::shared_ptr<SolverType> smoother_peak = lvl.level->get_smoother_peak();
 
        // do we have a peak-smoother?
        if(smoother_peak)
        {
          // apply peak-smoother
          if(!this->_apply_smooth_def(cur_lvl, *smoother_peak))
            return Status::aborted;
        }
        else
        {
          // no peak-smoother given
          // apply the pre- and post-smoothers instead
          std::shared_ptr<SolverType> smoother_pre = lvl.level->get_smoother_pre();
          std::shared_ptr<SolverType> smoother_post = lvl.level->get_smoother_post();
 
          if((smoother_pre)  && (!this->_apply_smooth_def(cur_lvl, *smoother_pre)))
            return Status::aborted;
          if((smoother_post) && (!this->_apply_smooth_def(cur_lvl, *smoother_post)))
            return Status::aborted;
        }
 
        TimeStamp bt;
        _toes.at(std::size_t(cur_lvl)) += bt.elapsed(at);
        double mpi_exec_reduction_stop(Statistics::get_time_mpi_execute_reduction());
        double mpi_exec_blas2_stop(Statistics::get_time_mpi_execute_blas2());
        double mpi_exec_blas3_stop(Statistics::get_time_mpi_execute_blas3());
        double mpi_exec_collective_stop(Statistics::get_time_mpi_execute_collective());
        double mpi_wait_stop_reduction(Statistics::get_time_mpi_wait_reduction());
        double mpi_wait_stop_blas2(Statistics::get_time_mpi_wait_blas2());
        double mpi_wait_stop_blas3(Statistics::get_time_mpi_wait_blas3());
        double mpi_wait_stop_collective(Statistics::get_time_mpi_wait_collective());
        _mpi_execs_reduction.at(std::size_t(cur_lvl)) += mpi_exec_reduction_stop - mpi_exec_reduction_start;
        _mpi_execs_blas2.at(std::size_t(cur_lvl)) += mpi_exec_blas2_stop - mpi_exec_blas2_start;
        _mpi_execs_blas3.at(std::size_t(cur_lvl)) += mpi_exec_blas3_stop - mpi_exec_blas3_start;
        _mpi_execs_collective.at(std::size_t(cur_lvl)) += mpi_exec_collective_stop - mpi_exec_collective_start;
        _mpi_waits_reduction.at(std::size_t(cur_lvl)) += mpi_wait_stop_reduction - mpi_wait_start_reduction;
        _mpi_waits_blas2.at(std::size_t(cur_lvl)) += mpi_wait_stop_blas2 - mpi_wait_start_blas2;
        _mpi_waits_blas3.at(std::size_t(cur_lvl)) += mpi_wait_stop_blas3 - mpi_wait_start_blas3;
        _mpi_waits_collective.at(std::size_t(cur_lvl)) += mpi_wait_stop_collective - mpi_wait_start_collective;
 
        return Status::success;
      }
 
      Status _apply_rest(const Index cur_lvl, bool cur_smooth)
      {
        const Index last_level = Math::min(_crs_level, _hierarchy->size_physical());
 
        // restriction loop: from current level to last level
        for(Index i(cur_lvl); i < last_level; ++i)
        {
          TimeStamp at;
          double mpi_exec_reduction_start(Statistics::get_time_mpi_execute_reduction());
          double mpi_exec_blas2_start(Statistics::get_time_mpi_execute_blas2());
          double mpi_exec_blas3_start(Statistics::get_time_mpi_execute_blas3());
          double mpi_exec_collective_start(Statistics::get_time_mpi_execute_collective());
          double mpi_wait_start_reduction(Statistics::get_time_mpi_wait_reduction());
          double mpi_wait_start_blas2(Statistics::get_time_mpi_wait_blas2());
          double mpi_wait_start_blas3(Statistics::get_time_mpi_wait_blas3());
          double mpi_wait_start_collective(Statistics::get_time_mpi_wait_collective());
 
          // get our fine and coarse levels
          LevelInfo& lvl_f = _hierarchy->_get_level_info(i);
 
          // get system matrix and filters
          const MatrixType& system_matrix   = lvl_f.level->get_system_matrix();
          const FilterType& system_filter_f = lvl_f.level->get_system_filter();
 
          // get our pre-smoother
          std::shared_ptr<SolverType> smoother = lvl_f.level->get_smoother_pre();
 
          // The following if-block ensures that we never apply the pre-smoother
          // when restricting from an intermediate peak-level in the F- or W-cycle,
          // as otherwise the following statements would discard the solution vector
          // approximation that has already been computed by prior restriction and
          // peak-smoothing operations.
          if(cur_smooth || (i > cur_lvl))
          {
            // apply pre-smoother if we have one
            if(smoother)
            {
              // apply pre-smoother
              TimeStamp stamp_smooth;
              Statistics::add_solver_expression(std::make_shared<ExpressionCallSmoother>(this->name(), smoother->name()));
              smoother->apply(lvl_f.vec_sol, lvl_f.vec_rhs);
              //if(!status_success(smoother->apply(lvl_f.vec_sol, lvl_f.vec_rhs)))
                //return Status::aborted;
              lvl_f.time_smooth += stamp_smooth.elapsed_now();
 
              // compute defect
              TimeStamp stamp_defect;
              system_matrix.apply(lvl_f.vec_def, lvl_f.vec_sol, lvl_f.vec_rhs, -DataType(1));
              lvl_f.time_defect += stamp_defect.elapsed_now();
            }
            else
            {
              // format solution
              lvl_f.vec_sol.format();
 
              // set our rhs as defect
              lvl_f.vec_def.copy(lvl_f.vec_rhs);
            }
          }
 
          // filter our defect
          system_filter_f.filter_def(lvl_f.vec_def);
 
          // get our transfer operator
          const TransferOperatorType* transfer_operator = lvl_f.level->get_transfer_operator();
          XASSERTM(transfer_operator != nullptr, "transfer operator is missing");
 
          // break loop?
          bool break_loop = false;
 
          // ghost transfer?
          if(transfer_operator->is_ghost())
          {
            // send restriction to parent processes and return
            TimeStamp stamp_rest;
            transfer_operator->rest_send(lvl_f.vec_def);
            lvl_f.time_transfer += stamp_rest.elapsed_now();
            break_loop = true;
          }
          else
          {
            // okay, get coarse level
            LevelInfo& lvl_c = _hierarchy->_get_level_info(i+1);
            const FilterType& system_filter_c = lvl_c.level->get_system_filter();
 
            // restrict onto coarse level
            //Statistics::add_solver_expression(std::make_shared<ExpressionRestriction>(this->name(), i));
            TimeStamp stamp_rest;
            transfer_operator->rest(lvl_f.vec_def, lvl_c.vec_rhs);
            lvl_f.time_transfer += stamp_rest.elapsed_now();
 
            // filter coarse defect
            system_filter_c.filter_def(lvl_c.vec_rhs);
          }
 
          // collect timings
          _toes.at((size_t)i) += at.elapsed_now();
          double mpi_exec_reduction_stop(Statistics::get_time_mpi_execute_reduction());
          double mpi_exec_blas2_stop(Statistics::get_time_mpi_execute_blas2());
          double mpi_exec_blas3_stop(Statistics::get_time_mpi_execute_blas3());
          double mpi_exec_collective_stop(Statistics::get_time_mpi_execute_collective());
          double mpi_wait_stop_reduction(Statistics::get_time_mpi_wait_reduction());
          double mpi_wait_stop_blas2(Statistics::get_time_mpi_wait_blas2());
          double mpi_wait_stop_blas3(Statistics::get_time_mpi_wait_blas3());
          double mpi_wait_stop_collective(Statistics::get_time_mpi_wait_collective());
          _mpi_execs_reduction.at((size_t)i) += mpi_exec_reduction_stop - mpi_exec_reduction_start;
          _mpi_execs_blas2.at((size_t)i) += mpi_exec_blas2_stop - mpi_exec_blas2_start;
          _mpi_execs_blas3.at((size_t)i) += mpi_exec_blas3_stop - mpi_exec_blas3_start;
          _mpi_execs_collective.at((size_t)i) += mpi_exec_collective_stop - mpi_exec_collective_start;
          _mpi_waits_reduction.at((size_t)i) += mpi_wait_stop_reduction - mpi_wait_start_reduction;
          _mpi_waits_blas2.at((size_t)i) += mpi_wait_stop_blas2 - mpi_wait_start_blas2;
          _mpi_waits_blas3.at((size_t)i) += mpi_wait_stop_blas3 - mpi_wait_start_blas3;
          _mpi_waits_collective.at((size_t)i) += mpi_wait_stop_collective - mpi_wait_start_collective;
 
          if(break_loop)
            break;
 
          // descent to prior level
        }
 
        return Status::success;
      }
 
      virtual Status _apply_prol(const Index cur_lvl, bool cur_smooth)
      {
        const Index last_level = Math::min(_crs_level, _hierarchy->size_physical());
 
        // prolongation loop: from last level to current level
        for(Index i(last_level); i > cur_lvl;)
        {
          --i;
          TimeStamp at;
          double mpi_exec_reduction_start(Statistics::get_time_mpi_execute_reduction());
          double mpi_exec_blas2_start(Statistics::get_time_mpi_execute_blas2());
          double mpi_exec_blas3_start(Statistics::get_time_mpi_execute_blas3());
          double mpi_exec_collective_start(Statistics::get_time_mpi_execute_collective());
          double mpi_wait_start_reduction(Statistics::get_time_mpi_wait_reduction());
          double mpi_wait_start_blas2(Statistics::get_time_mpi_wait_blas2());
          double mpi_wait_start_blas3(Statistics::get_time_mpi_wait_blas3());
          double mpi_wait_start_collective(Statistics::get_time_mpi_wait_collective());
 
          // get our fine level
          LevelInfo& lvl_f = _hierarchy->_get_level_info(i);
 
          // get our transfer operator
          const TransferOperatorType* transfer_operator = lvl_f.level->get_transfer_operator();
          XASSERTM(transfer_operator != nullptr, "transfer operator is missing");
 
          // ghost transfer?
          if(transfer_operator->is_ghost())
          {
            // receive prolongation
            TimeStamp stamp_prol;
            transfer_operator->prol_recv(lvl_f.vec_cor);
            lvl_f.time_transfer += stamp_prol.elapsed_now();
          }
          else
          {
            // get coarse level
            LevelInfo& lvl_c = _hierarchy->_get_level_info(i+1);
 
            // prolongate
            TimeStamp stamp_prol;
            transfer_operator->prol(lvl_f.vec_cor, lvl_c.vec_sol);
            lvl_f.time_transfer += stamp_prol.elapsed_now();
          }
 
          // get system matrix and filters
          const MatrixType& system_matrix   = lvl_f.level->get_system_matrix();
          const FilterType& system_filter_f = lvl_f.level->get_system_filter();
 
          // apply correction filter
          system_filter_f.filter_cor(lvl_f.vec_cor);
 
          // initialize omega for (adaptive) coarse grid correction
          DataType omega_cgc = DataType(1);
 
          // use some sort of adaptive cgc?
          if(_adapt_cgc != MultiGridAdaptCGC::Fixed)
          {
            // multiply coarse grid correction by system matrix
            TimeStamp stamp_defect;
            system_matrix.apply(lvl_f.vec_tmp, lvl_f.vec_cor);
            lvl_f.time_defect += stamp_defect.elapsed_now();
 
            // apply filter
            system_filter_f.filter_def(lvl_f.vec_tmp);
 
            // compute adaptive omega
            switch(_adapt_cgc)
            {
            case MultiGridAdaptCGC::MinEnergy:
              omega_cgc = lvl_f.vec_def.dot(lvl_f.vec_cor)
                        / lvl_f.vec_tmp.dot(lvl_f.vec_cor);
              break;
 
            case MultiGridAdaptCGC::MinDefect:
              omega_cgc = lvl_f.vec_def.dot(lvl_f.vec_tmp)
                        / lvl_f.vec_tmp.dot(lvl_f.vec_tmp);
              break;
 
            default:
              // This default block is required to make the GCC happy...
              break;
            }
          }
 
          // update our solution vector
          lvl_f.vec_sol.axpy(lvl_f.vec_cor, omega_cgc);
 
          // get our post-smoother
          std::shared_ptr<SolverType> smoother = lvl_f.level->get_smoother_post();
 
          // apply post-smoother if we have one
          if(smoother && (cur_smooth || (i > cur_lvl)))
          {
            // compute new defect
            if(_adapt_cgc != MultiGridAdaptCGC::Fixed)
            {
              // We have used some sort of adaptive coarse grid correction.
              // In this case, the current solution vector is
              //    sol_new = sol_old + omega * cor
              //
              // Moreover, we have the following vectors at hand:
              //    1) tmp = A*cor
              //    2) def_old = rhs - A*sol_old
              //
              // So we can apply the same "trick" as in the PCG method to update our
              // defect vector instead of recomputing it by a matrix-vector mult, i.e.:
              //   def_new = rhs - A * sol_new
              //           = rhs - A * (sol_old + omega * cor)
              //           = rhs - A *  sol_old - omega * A * cor
              //           = def_old            - omega * tmp
              lvl_f.vec_def.axpy(lvl_f.vec_tmp, -omega_cgc);
            }
            else
            {
              // we have to recompute the defect by (rhs - A*sol_new)
              TimeStamp stamp_defect;
              system_matrix.apply(lvl_f.vec_def, lvl_f.vec_sol, lvl_f.vec_rhs, -DataType(1));
              lvl_f.time_defect += stamp_defect.elapsed_now();
 
              // apply defect filter
              system_filter_f.filter_def(lvl_f.vec_def);
            }
 
            // apply post-smoother
            Statistics::add_solver_expression(std::make_shared<ExpressionCallSmoother>(this->name(), smoother->name()));
            TimeStamp stamp_smooth;
            smoother->apply(lvl_f.vec_cor, lvl_f.vec_def);
            //if(!status_success(smoother->apply(lvl_f.vec_cor, lvl_f.vec_def)))
              //return Status::aborted;
            lvl_f.time_smooth += stamp_smooth.elapsed_now();
 
            // update solution vector
            lvl_f.vec_sol.axpy(lvl_f.vec_cor);
          }
 
          _toes.at((size_t)i) += at.elapsed_now();
          double mpi_exec_reduction_stop(Statistics::get_time_mpi_execute_reduction());
          double mpi_exec_blas2_stop(Statistics::get_time_mpi_execute_blas2());
          double mpi_exec_blas3_stop(Statistics::get_time_mpi_execute_blas3());
          double mpi_exec_collective_stop(Statistics::get_time_mpi_execute_collective());
          double mpi_wait_stop_reduction(Statistics::get_time_mpi_wait_reduction());
          double mpi_wait_stop_blas2(Statistics::get_time_mpi_wait_blas2());
          double mpi_wait_stop_blas3(Statistics::get_time_mpi_wait_blas3());
          double mpi_wait_stop_collective(Statistics::get_time_mpi_wait_collective());
          _mpi_execs_reduction.at((size_t)i) += mpi_exec_reduction_stop - mpi_exec_reduction_start;
          _mpi_execs_blas2.at((size_t)i) += mpi_exec_blas2_stop - mpi_exec_blas2_start;
          _mpi_execs_blas3.at((size_t)i) += mpi_exec_blas3_stop - mpi_exec_blas3_start;
          _mpi_execs_collective.at((size_t)i) += mpi_exec_collective_stop - mpi_exec_collective_start;
          _mpi_waits_reduction.at((size_t)i) += mpi_wait_stop_reduction - mpi_wait_start_reduction;
          _mpi_waits_blas2.at((size_t)i) += mpi_wait_stop_blas2 - mpi_wait_start_blas2;
          _mpi_waits_blas3.at((size_t)i) += mpi_wait_stop_blas3 - mpi_wait_start_blas3;
          _mpi_waits_collective.at((size_t)i) += mpi_wait_stop_collective - mpi_wait_start_collective;
 
          // ascend to next level
        }
 
        return Status::success;
      }
    }; // class MultiGrid<...>
 
    template<
      typename SystemMatrix_,
      typename SystemFilter_,
      typename TransferOperator_>
    std::shared_ptr<MultiGrid<SystemMatrix_, SystemFilter_, TransferOperator_>> new_multigrid(
      std::shared_ptr<MultiGridHierarchy<SystemMatrix_, SystemFilter_, TransferOperator_>> hierarchy,
      MultiGridCycle cycle = MultiGridCycle::V,
      int top_level = 0,
      int crs_level = -1)
    {
      return std::make_shared<MultiGrid<SystemMatrix_, SystemFilter_, TransferOperator_>>
        (hierarchy, cycle, top_level, crs_level);
    }
  } // namespace Solver
} // namespace FEAT