Support for more than 3 indices

[TakeTwiceDailey]

Currently, it seems that @parallel and @parallel_indices are limited to launching kernels with 3 or less indices.

This can be worked around by instead launching over a linear index for simple kernels, like this simple copy operation:

@parallel_indices (linear_index) function copy!(Un,U)
    U[linear_index]  = Un[linear_index]
    return
end

which works with arbitrary rank arrays. My use case for this is a 3D simulation that is set up with many connected blocks. I would like one index to set the block number, and the other three indices to label the 3D point in the array of that block. I can set up simple kernels like the above, but I can't do more complicated things like a subset operation:

@parallel_indices (block,xi,yi,zi) function copy!(Un,U)
    U[block,xi,yi,zi]  = Un[block,xi,yi,zi]
    return
end

@parallel (1:6,5:n-4,5:n-4,5:n-4) copy!(Un,U)

that only affects the interior region of all 6 blocks for example. Instead, I have to write a loop that launches 6 separate @parallel kernels.

Another thing that would be nice is to launch a kernel over a given list of indices. Sometimes a simple Cartesian subset is not adequate, and I need to launch a kernel over a more complicated subset of indices in an array. I could just run a kernel over all indices and set up an if statement in the kernel to exit if it isn't in the given subset, but I am worried about branch divergence on GPU.

Thanks!

[omlins]

Hi @TakeTwiceDailey , sorry for the delay in the reply. Yes, there is a limit to three or less indices, as in the vendor libraries like CUDA. Adding support for n-dimensional indices does not seem to be possible without introducing overhead. In fact, I've seen that in KernelAbstractions, there are plans to support again three or less indices, exactly for this reason. Therefore we do not want to add support for n-dimensional indices.

There is normally an other solution than arrays with more than 3 indices. A great solution is to create, for example, a tuple or a named tuple of arrays. Another solution is using CellArrays, but when the tuple solution works, it's easier and better.

In your particular case, how many blocks do you have? And how big is the domain of the other three indices?

Regarding the second issue, we might be able to do some convenience solution for this at some point; however, in any case, working with non-contiguous indices will never be high-performant and should not be a key pattern in an application. Now that said, you can indeed just launch enough threads to reassure to capture all the indices, and for the indices you don't have to do anything return: if (index !in indices) return end. Or alternatively you could pass the list of n indices to your kernel and then launch the kernel as @parallel (1:n) kernel(indices,...) and map in the kernel the index: i = indices[ix], etc. This way you're launching not more threads than required. You might prefer the second solution.

[TakeTwiceDailey]

No worries, thanks for the reply! If performance dictates 3 or less indices, then I should conform to that, not the other way around for sure. Since I posted this, I was able to code in a more clever way so as to avoid the need for the subset operation I described. I do however rely on linear_index kernels to code the main Runge Kutta stages so as to allow rank 4 arrays, which are indeed CellArrays.

Since you asked, my setup includes 7 blocks, each of which can have arbitrarily fine grids. My issue when I posted this is that each of the blocks uses a particular kind of finite differencing stencil called SBP operators, whose stencils are the same in the interior, but change near the boundaries. I was trying to get the speed of the common code in the interior by launching only in that region, but I found that I can basically get the same speed by storing a small StaticArray of stencil coefficients and selecting one by proximity to the boundary.

As for the second point, I was worried about warp divergence on GPU, which should lead to a huge slowdown, which to my understanding happens when you have a lot of if logic in a kernel. As it turns out however, there is not a huge penalty to having at least a single if (in main region) elseif (in smaller region) end in the main kernel at least in my implementation, so I have settled on this.

It is good to learn where the computational limitations are, so I appreciate your reply.

[omlins]

Great, I'm glad you found a solution!