Parallel version of for loops and dfor

We offer two main flavours of parallel for macros: a parallel dfor and a dfor mapping onto SIMT, i.e.

bare of internal and external dependencies. The latter is subject to some restrictions:

• Don't use virtual function calls.
• Take into account that the loop body works with copies of local variables (you can still use pointers).
• Don't use synchronisation primitives.

All macros want to have a C++ identifier counter. This counter will be a well-defined integer vector within the loop body. The counters always start at 0, but you can shift them obviously internally to have another offset. The max argument is either an integer, or you can paste in a Vector. The usage for the latter looks similar to

parallelDforWithSchedulerInstructions(volume,
(tarch::la::Vector<2,int>({numberOfVolumesPerAxisInPatch+2*solversHaloSize,numberOfVolumesPerAxisInPatch})),
loopPlacement)

By using a vector, you can have rectangular iteration ranges and not only squares/cubes. The syntax with the additional brackets is a little bit annoying, but some LLVM front-ends otherwise get confused. Try without them, but don't be surprised if it doesn't work.

The mapping onto your variables of choice is very similar to what we see a lot in CUDA and SYCL kernels. In the simplest case, you write

parallelDforWithSchedulerInstructions(volume, ...) {
  int firstIndex = volume(0);
  int secondIndex = volume(1);
} endParallelDfor

In this context, it is also interesting to study tarch::la::slice() which allows you to extract whole subvectors from the index vector in one rush.

Generic realisation remarks

• I convert the max argument always into counter##Max and then work against the vector. This way, I can support both vectors and scalars. In the former case, the copy constructor is used. In the latter case the constructor accepting a scalar.
• I embed the whole loop into another scope (brackets). Otherwise, I run into issues if users write things like

  parallelDfor(i,3) {
    [...]
  } endParallelDfor
  parallelDfor(i,3) {
    [...]
  } endParallelDfor

  as the corresponding vector iMax then would be declared twice.

OpenMP

OpenMP relies on pragmas around the nested for loops, which is tricky as parallelDfor itself is a pragma. I hijack the construct

         _Pragma( "my embedded pragma" )

which is, to the best of my knowledge, a C99 instruction and not C++17. So I simply assume and hope that this works. Note that the _Pragma statement yields a pragma, i.e. you are not allowed to write

         _Pragma( "#pragma my embedded pragma" )

as this would unfold into

         #pragma #pragma my embedded pragma

Further to that, variables of the pragma are not replaced in the unfolded pragma. Therefore, I have to go this unorthodox way with the helper loop guard.

For the available OpenMP variants, it neither makes sense to have parallel fors nor does it work properly:

• If we run the kernel within a task (which we typically do), the task cannot spawn out onto multiple cores.
• I found no way to have multiple pragmas in a row where we toggle through the options via if.

If we use a taskloop in OpenMP, the default visibility of all variables will be taskfirst. Furthermore, the created copies will be subject to a copy constructor and their destructor. This means that we alter the visibility when we switch from parallel for or no parallelisation into a taskloop. Therefore, I manually change the default visibility to shared.

TBB

TBB is, in my opinion, the most straightforward implementation, as it supports d-dimensional ranges: It has loop collapsing built-in. I use the TBB ranges, embed a nested for loop into the range, and then puzzle the counter of interest together from the ranges indices again.

SYCL

Different to TBB, SYCL has this accelerator paradigm in mind, even though we submit to the host queue. Therefore, we have to capture the variables via =, we have to ensure all routines (such as dLinearised() are marked as GPU offloadable), we do not have any dynamic scheduling, grain sizes, and so forth. Most importantly, we cannot call virtual functions or functors from SYCL kernels.

For all of these reasons, we do not support parallel for for SYCL, and we do not support SIMT parallelism either on the host. In theory, the latter should work straightforwardly, but I got major issues with copy constructors et al trying to go down this route, so I eventually dismissed it as a stupid idea.

SIMT loops are always to be embedded into parallel loops. We do not embed the other way round. Codes like the batched kernels will use SIMT loops right from the start. So it would be important that the SIMT loops exploit parallelism. If we make compromises, we have to make compromises on the parallel loop level. Once SYCL abandons the need for call by value in the capture clause of the functors, we should be able to use the SYCLised simt macros for Peano's loops.