Loop.h

Go to the documentation of this file.

// This file is part of the Peano project. For conditions of distribution and
// use, please see the copyright notice at www.peano-framework.org
#pragma once
 
#include <bitset>
 
#include "config.h"
#include "peano4/utils/Globals.h"
#include "tarch/compiler/CompilerSpecificSettings.h"
#include "tarch/accelerator/accelerator.h"
#include "tarch/la/Vector.h"
 
 
#if defined(SharedTBB)
  #if defined(SharedTBBExtension) and !defined(TBB_PREVIEW_BLOCKED_RANGE_ND)
    #include "tbb/blocked_rangeNd.h"
  #else
    #include <oneapi/tbb/blocked_range.h>
  #endif
 
  #include <oneapi/tbb/blocked_range2d.h>
  #include <oneapi/tbb/blocked_range3d.h>
  #include <oneapi/tbb/parallel_for.h>
 
  #ifdef TBB_OLD_BLOCKED_RANGE_ND
  namespace tbb {
  namespace detail {
  namespace d1 {
    using blocked_nd_range = blocked_rangeNd;
  }}}
  #endif
#endif
 
#if defined(SharedSYCL)
#include "tarch/multicore/sycl/Core.h"
 
#define CPUGPUMethod SYCL_EXTERNAL
#else
#define CPUGPUMethod
#endif
 
namespace peano4 {
  namespace utils {
    enum class LoopPlacement { Serial, Nested, SpreadOut };
 
    typedef std::bitset<Dimensions> LoopDirection;
 
    CPUGPUMethod void dInc(tarch::la::Vector<Dimensions, int>& counter, int max);
 
    void dDec(tarch::la::Vector<Dimensions, int>& counter, int max);
 
    void dInc(tarch::la::Vector<Dimensions, int>& counter, const tarch::la::Vector<Dimensions, int>& max);
 
    void dIncByVector(tarch::la::Vector<Dimensions, int>& counter, int max, int increment);
 
    void dIncByScalar(tarch::la::Vector<Dimensions, int>& counter, int max, int increment);
 
    void dInc(tarch::la::Vector<Dimensions, int>& counter, int max, int doNotExamine);
 
    void dInc(tarch::la::Vector<Dimensions, int>& counter, int max, LoopDirection& direction);
 
    CPUGPUMethod int dCmp(const tarch::la::Vector<Dimensions, int>& counter, int max);
 
    int dCmp(const tarch::la::Vector<Dimensions, int>& counter, const tarch::la::Vector<Dimensions, int>& max);
 
    bool dCmpLinearOrder(
      const tarch::la::Vector<Dimensions, int>& counter, const tarch::la::Vector<Dimensions, int>& max
    );
 
    CPUGPUMethod int dLinearised(const tarch::la::Vector<Dimensions, int>& counter, int max);
    CPUGPUMethod int dLinearised(const std::bitset<Dimensions>& counter);
 
    int d2Linearised(const tarch::la::Vector<2, int>& counter, int max);
 
    int d3Linearised(const tarch::la::Vector<3, int>& counter, int max);
 
    int dLinearisedWithoutLookup(const tarch::la::Vector<Dimensions, int>& counter, int max);
 
    tarch::la::Vector<Dimensions, int> dDelinearised(int value, int max);
 
    tarch::la::Vector<Dimensions, int> dDelinearisedWithoutLookup(int value, int max);
 
 
    void setupLookupTableForDLinearised();
    void setupLookupTableForDDelinearised();
 
    CPUGPUMethod tarch::la::Vector<Dimensions, int> dStartVector();
 
    tarch::la::Vector<Dimensions, int> dStartVector(int dim, int value);
 
    tarch::la::Vector<Dimensions, int> dStartVector(int max, const LoopDirection& direction);
  }
}
 
 
std::string toString(peano4::utils::LoopPlacement);
 
 
#define dfor(counter, max) \
  for (tarch::la::Vector<Dimensions, int> counter = peano4::utils::dStartVector(); peano4::utils::dCmp(counter, max); \
       peano4::utils::dInc(counter, max))
 
 
#define parallelFor(counter, max) \
  parallelForWithSchedulerInstructions(counter, max, ::peano4::utils::LoopPlacement::SpreadOut)
 
#define simtFor(counter, max) \
  simtForWithSchedulerInstructions(counter, max, ::peano4::utils::LoopPlacement::SpreadOut)
 
#define parallelDfor(counter, max) \
  parallelDforWithSchedulerInstructions(counter, max, ::peano4::utils::LoopPlacement::SpreadOut)
 
#if Dimensions == 2
#define parallelDforWithSchedulerInstructions(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions2d(counter, max, loopPlacement)
 
#define parallelDPlusOneForWithSchedulerInstructions(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions3d(counter, max, loopPlacement)
 
#define parallelDPlusTwoForWithSchedulerInstructions(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions4d(counter, max, loopPlacement)
 
#define simtDforWithSchedulerInstructions(counter, max, loopPlacement) \
  simtDforWithSchedulerInstructions2d(counter, max, loopPlacement)
 
#define simtDPlusOneForWithSchedulerInstructions(counter, max, loopPlacement) \
  simtDforWithSchedulerInstructions3d(counter, max, loopPlacement)
 
#define simtDPlusTwoForWithSchedulerInstructions(counter, max, loopPlacement) \
  simtDforWithSchedulerInstructions4d(counter, max, loopPlacement)
#elif Dimensions == 3
#define parallelDforWithSchedulerInstructions(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions3d(counter, max, loopPlacement)
 
#define parallelDPlusOneForWithSchedulerInstructions(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions4d(counter, max, loopPlacement)
 
#define parallelDPlusTwoForWithSchedulerInstructions(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions5d(counter, max, loopPlacement)
 
#define simtDforWithSchedulerInstructions(counter, max, loopPlacement) \
  simtDforWithSchedulerInstructions3d(counter, max, loopPlacement)
 
#define simtDPlusOneForWithSchedulerInstructions(counter, max, loopPlacement) \
  simtDforWithSchedulerInstructions4d(counter, max, loopPlacement)
 
#define simtDPlusTwoForWithSchedulerInstructions(counter, max, loopPlacement) \
  simtDforWithSchedulerInstructions5d(counter, max, loopPlacement)
#endif
 
 
#if defined(SharedOMP) and defined(OpenMPTaskGroup)
#include "omp/Loop.h"
#elif defined(SharedOMP) and !defined(OpenMPTaskGroup)
#include "omp/LoopWithoutTaskLoop.h"
#elif defined(SharedTBB)
#include "tbb/Loop.h"
#elif defined(SharedCPP)
//and __cplusplus >= 202101L
#include "cpp/Loop.h"
#elif defined(SharedSYCL)
#include "sycl/Loop.h"
#else
#define parallelForWithSchedulerInstructions(counter, max, loopPlacement) \
  for (int counter = 0; counter < max; counter++) {
 
#define parallelDforWithSchedulerInstructions2d(counter, max, loopPlacement) \
  { \
    tarch::la::Vector<2, int> counter##Max(max); \
    for (int counter##0 = 0; counter##0 < counter##Max(0); counter##0 ++) \
      for (int counter##1 = 0; counter##1 < counter##Max(1); counter##1 ++) { \
        tarch::la::Vector<2, int> counter = {counter##0, counter##1};
 
#define parallelDforWithSchedulerInstructions3d(counter, max, loopPlacement) \
  { \
    tarch::la::Vector<3, int> counter##Max(max); \
    for (int counter##0 = 0; counter##0 < counter##Max(0); counter##0 ++) \
      for (int counter##1 = 0; counter##1 < counter##Max(1); counter##1 ++) \
        for (int counter##2 = 0; counter##2 < counter##Max(2); counter##2 ++) { \
          tarch::la::Vector<3, int> counter = {counter##0, counter##1, counter##2};
 
#define parallelDforWithSchedulerInstructions4d(counter, max, loopPlacement) \
  { \
    tarch::la::Vector<4, int> counter##Max(max); \
    for (int counter##0 = 0; counter##0 < counter##Max(0); counter##0 ++) \
      for (int counter##1 = 0; counter##1 < counter##Max(1); counter##1 ++) \
        for (int counter##2 = 0; counter##2 < counter##Max(2); counter##2 ++) \
          for (int counter##3 = 0; counter##3 < counter##Max(3); counter##3 ++) { \
            tarch::la::Vector<4, int> counter = {counter##0, counter##1, counter##2, counter##3};
 
#define parallelDforWithSchedulerInstructions5d(counter, max, loopPlacement) \
  { \
    tarch::la::Vector<5, int> counter##Max(max); \
    for (int counter##0 = 0; counter##0 < counter##Max(0); counter##0 ++) \
      for (int counter##1 = 0; counter##1 < counter##Max(1); counter##1 ++) \
        for (int counter##2 = 0; counter##2 < counter##Max(2); counter##2 ++) \
          for (int counter##3 = 0; counter##3 < counter##Max(3); counter##3 ++) \
            for (int counter##4 = 0; counter##4 < counter##Max(4); counter##4 ++) { \
              tarch::la::Vector<5, int> counter = {counter##0, counter##1, counter##2, counter##3, counter##4};
 
#define endParallelFor }
#define endParallelDfor \
  } \
  }
 
 
#define simtForWithSchedulerInstructions(counter, max, loopPlacement) \
  for (int counter = 0; counter < max; counter++) {
 
#define endSimtFor endParallelFor
 
#define simtDforWithSchedulerInstructions2d(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions2d(counter, max, loopPlacement)
 
#define simtDforWithSchedulerInstructions3d(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions3d(counter, max, loopPlacement)
 
#define simtDforWithSchedulerInstructions4d(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions4d(counter, max, loopPlacement)
 
#define simtDforWithSchedulerInstructions5d(counter, max, loopPlacement) \
  parallelDforWithSchedulerInstructions5d(counter, max, loopPlacement)
 
#define endSimtDfor endParallelDfor
 
#endif
 
 
#define dfor4(counter) \
  for (int counter##Scalar = 0; counter##Scalar < FOUR_POWER_D; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = Dimensions - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 4; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
 
#define dfor3(counter) \
  for (int counter##Scalar = 0; counter##Scalar < ThreePowerD; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = Dimensions - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 3; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
 
#define dfor5(counter) \
  for (int counter##Scalar = 0; counter##Scalar < FIVE_POWER_D; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = Dimensions - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 5; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
#define dfor7(counter) \
  for (int counter##Scalar = 0; counter##Scalar < SEVEN_POWER_D; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = Dimensions - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 7; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
 
#define dfor9(counter) \
  for (int counter##Scalar = 0; counter##Scalar < NINE_POWER_D; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = Dimensions - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 9; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
#define d2for(counter, max) \
  for (int counter##Scalar = 0; counter##Scalar < tarch::la::aPowI(2, max); counter##Scalar++) { \
    tarch::la::Vector<2, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = 2 - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= max; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
#define d2for2(counter) \
  for (int counter##Scalar = 0; counter##Scalar < tarch::la::aPowI(2, 2); counter##Scalar++) { \
    tarch::la::Vector<2, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = 2 - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 2; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
#define d3for2(counter) \
  for (int counter##Scalar = 0; counter##Scalar < tarch::la::aPowI(3, 2); counter##Scalar++) { \
    tarch::la::Vector<3, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = 3 - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 2; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
 
#define d3for(counter, max) \
  for (int counter##Scalar = 0; counter##Scalar < tarch::la::aPowI(3, max); counter##Scalar++) { \
    tarch::la::Vector<3, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = 3 - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= max; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
 
#define d2for3(counter) \
  for (int counter##Scalar = 0; counter##Scalar < tarch::la::aPowI(2, 3); counter##Scalar++) { \
    tarch::la::Vector<2, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = 2 - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 3; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
#define d3for3(counter) \
  for (int counter##Scalar = 0; counter##Scalar < tarch::la::aPowI(3, 3); counter##Scalar++) { \
    tarch::la::Vector<3, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = 3 - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 3; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
 
/*
 * bit flipping used for Dimensions = 2, and Dimensions = 3
 * for more information about the idea principle used refer to https://opt-patterns.wiki.tum.de/dfor
 */
#if Dimensions == 2
#define dfor2(counter) \
  for (int counter##Scalar = 0, AA##counter = 0, BB##counter = 0; counter##Scalar < TwoPowerD; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    counter(0)  = AA##counter; \
    counter(1)  = BB##counter; \
    AA##counter = !AA##counter; \
    BB##counter = !(AA##counter ^ BB##counter);
 
#elif Dimensions == 3
#define dfor2(counter) \
  for (int counter##Scalar = 0, AA##counter = 0, BB##counter = 0, CC##counter = 0; counter##Scalar < TwoPowerD; \
       counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    counter(0)  = AA##counter; \
    counter(1)  = BB##counter; \
    counter(2)  = CC##counter; \
    AA##counter = !AA##counter; \
    BB##counter = !(AA##counter ^ BB##counter); \
    CC##counter = CC##counter || (!AA##counter && !BB##counter && !CC##counter);
 
#else
#define dfor2(counter) \
  for (int counter##Scalar = 0; counter##Scalar < TwoPowerD; counter##Scalar++) { \
    tarch::la::Vector<Dimensions, int> counter; \
    { \
      int copy##counter##Scalar = counter##Scalar; \
      for (int counter##ddd = Dimensions - 1; counter##ddd >= 0; counter##ddd--) { \
        int counter##aPowI = 1; \
        for (int counter##jjj = 0; counter##jjj < counter##ddd; counter##jjj++) { \
          counter##aPowI *= 2; \
        } \
        counter(counter##ddd) = copy##counter##Scalar / counter##aPowI; \
        copy##counter##Scalar -= counter(counter##ddd) * counter##aPowI; \
      } \
    }
 
#endif
 
 
#define enddforx }
 
 
#define dfore(counter, max, dim, value) \
  for (tarch::la::Vector<Dimensions, int> counter = peano4::utils::dStartVector(dim, value); \
       peano4::utils::dCmp(counter, max); \
       peano4::utils::dInc(counter, max, dim))
 
 
#define zfor(counter, max, direction) \
  { \
    for (tarch::la::Vector<Dimensions, int> counter = peano4::utils::dStartVector(max, direction); \
         peano4::utils::dCmp(counter, max); \
         peano4::utils::dInc(counter, max, direction)) {
 
 
/*
 * zfor3 is an optimized version of zfor for max = 3
 * A lookup table is used for dim=2 and dim=3, for higher dimensions
 * the standard zfor is used instead
 */
#if Dimensions == 2
// first entry is direction, second entry is cell, third entry is tuples
static const int lookupzfor3[4][9][2] = {
  {{2, 2}, {1, 2}, {0, 2}, {0, 1}, {1, 1}, {2, 1}, {2, 0}, {1, 0}, {0, 0}},
  {{0, 2}, {1, 2}, {2, 2}, {2, 1}, {1, 1}, {0, 1}, {0, 0}, {1, 0}, {2, 0}},
  {{2, 0}, {1, 0}, {0, 0}, {0, 1}, {1, 1}, {2, 1}, {2, 2}, {1, 2}, {0, 2}},
  {{0, 0}, {1, 0}, {2, 0}, {2, 1}, {1, 1}, {0, 1}, {0, 2}, {1, 2}, {2, 2}}};
 
#define zfor3(counter, direction) \
  { \
    tarch::la::Vector<Dimensions, int> counter; \
    int                                counter##initDir = static_cast<int>(direction.to_ulong()); \
    for (int counter##i = 0; counter##i < 9; ++counter##i) { \
      counter(0) = lookupzfor3[counter##initDir][counter##i][0]; \
      counter(1) = lookupzfor3[counter##initDir][counter##i][1];
 
#elif Dimensions == 3
static const int lookupzfor3[8][27][3] = {
  {{2, 2, 2}, {1, 2, 2}, {0, 2, 2}, {0, 1, 2}, {1, 1, 2}, {2, 1, 2}, {2, 0, 2}, {1, 0, 2}, {0, 0, 2},
   {0, 0, 1}, {1, 0, 1}, {2, 0, 1}, {2, 1, 1}, {1, 1, 1}, {0, 1, 1}, {0, 2, 1}, {1, 2, 1}, {2, 2, 1},
   {2, 2, 0}, {1, 2, 0}, {0, 2, 0}, {0, 1, 0}, {1, 1, 0}, {2, 1, 0}, {2, 0, 0}, {1, 0, 0}, {0, 0, 0}},
  {{0, 2, 2}, {1, 2, 2}, {2, 2, 2}, {2, 1, 2}, {1, 1, 2}, {0, 1, 2}, {0, 0, 2}, {1, 0, 2}, {2, 0, 2},
   {2, 0, 1}, {1, 0, 1}, {0, 0, 1}, {0, 1, 1}, {1, 1, 1}, {2, 1, 1}, {2, 2, 1}, {1, 2, 1}, {0, 2, 1},
   {0, 2, 0}, {1, 2, 0}, {2, 2, 0}, {2, 1, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 0}, {1, 0, 0}, {2, 0, 0}},
  {{2, 0, 2}, {1, 0, 2}, {0, 0, 2}, {0, 1, 2}, {1, 1, 2}, {2, 1, 2}, {2, 2, 2}, {1, 2, 2}, {0, 2, 2},
   {0, 2, 1}, {1, 2, 1}, {2, 2, 1}, {2, 1, 1}, {1, 1, 1}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {2, 0, 1},
   {2, 0, 0}, {1, 0, 0}, {0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {2, 1, 0}, {2, 2, 0}, {1, 2, 0}, {0, 2, 0}},
  {{0, 0, 2}, {1, 0, 2}, {2, 0, 2}, {2, 1, 2}, {1, 1, 2}, {0, 1, 2}, {0, 2, 2}, {1, 2, 2}, {2, 2, 2},
   {2, 2, 1}, {1, 2, 1}, {0, 2, 1}, {0, 1, 1}, {1, 1, 1}, {2, 1, 1}, {2, 0, 1}, {1, 0, 1}, {0, 0, 1},
   {0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {2, 1, 0}, {1, 1, 0}, {0, 1, 0}, {0, 2, 0}, {1, 2, 0}, {2, 2, 0}},
  {{2, 2, 0}, {1, 2, 0}, {0, 2, 0}, {0, 1, 0}, {1, 1, 0}, {2, 1, 0}, {2, 0, 0}, {1, 0, 0}, {0, 0, 0},
   {0, 0, 1}, {1, 0, 1}, {2, 0, 1}, {2, 1, 1}, {1, 1, 1}, {0, 1, 1}, {0, 2, 1}, {1, 2, 1}, {2, 2, 1},
   {2, 2, 2}, {1, 2, 2}, {0, 2, 2}, {0, 1, 2}, {1, 1, 2}, {2, 1, 2}, {2, 0, 2}, {1, 0, 2}, {0, 0, 2}},
  {{0, 2, 0}, {1, 2, 0}, {2, 2, 0}, {2, 1, 0}, {1, 1, 0}, {0, 1, 0}, {0, 0, 0}, {1, 0, 0}, {2, 0, 0},
   {2, 0, 1}, {1, 0, 1}, {0, 0, 1}, {0, 1, 1}, {1, 1, 1}, {2, 1, 1}, {2, 2, 1}, {1, 2, 1}, {0, 2, 1},
   {0, 2, 2}, {1, 2, 2}, {2, 2, 2}, {2, 1, 2}, {1, 1, 2}, {0, 1, 2}, {0, 0, 2}, {1, 0, 2}, {2, 0, 2}},
  {{2, 0, 0}, {1, 0, 0}, {0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {2, 1, 0}, {2, 2, 0}, {1, 2, 0}, {0, 2, 0},
   {0, 2, 1}, {1, 2, 1}, {2, 2, 1}, {2, 1, 1}, {1, 1, 1}, {0, 1, 1}, {0, 0, 1}, {1, 0, 1}, {2, 0, 1},
   {2, 0, 2}, {1, 0, 2}, {0, 0, 2}, {0, 1, 2}, {1, 1, 2}, {2, 1, 2}, {2, 2, 2}, {1, 2, 2}, {0, 2, 2}},
  {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {2, 1, 0}, {1, 1, 0}, {0, 1, 0}, {0, 2, 0}, {1, 2, 0}, {2, 2, 0},
   {2, 2, 1}, {1, 2, 1}, {0, 2, 1}, {0, 1, 1}, {1, 1, 1}, {2, 1, 1}, {2, 0, 1}, {1, 0, 1}, {0, 0, 1},
   {0, 0, 2}, {1, 0, 2}, {2, 0, 2}, {2, 1, 2}, {1, 1, 2}, {0, 1, 2}, {0, 2, 2}, {1, 2, 2}, {2, 2, 2}}};
 
#define zfor3(counter, direction) \
  { \
    tarch::la::Vector<Dimensions, int> counter; \
    int                                counter##initDir = static_cast<int>(direction.to_ulong()); \
    for (int counter##i = 0; counter##i < 27; ++counter##i) { \
      counter(0) = lookupzfor3[counter##initDir][counter##i][0]; \
      counter(1) = lookupzfor3[counter##initDir][counter##i][1]; \
      counter(2) = lookupzfor3[counter##initDir][counter##i][2];
 
#else
#define zfor3(counter, direction) zfor(counter, 3, direction)
 
#endif
 
 
#define endzfor \
  } \
  }