PyOpenCL - Multi-dimensional reduction kernel

Question

I'm a total newbie to OpenCL.

I'm trying to code a reduction kernel that sums along one axis for a multi-dimensional array. I have stumbled upon that code which comes from here: https://tmramalho.github.io/blog/2014/06/16/parallel-programming-with-opencl-and-python-parallel-reduce/

__kernel void reduce(__global float *a, __global float *r, __local float *b) {
        uint gid = get_global_id(0);
        uint wid = get_group_id(0);
        uint lid = get_local_id(0);
        uint gs = get_local_size(0);

        b[lid] = a[gid];

        barrier(CLK_LOCAL_MEM_FENCE);
        
        for(uint s = gs/2; s > 0; s >>= 1) {
          if(lid < s) {
            b[lid] += b[lid+s];
          }
          barrier(CLK_LOCAL_MEM_FENCE);
        }
        if(lid == 0) r[wid] = b[lid];
}

I don't understand the for loop part. I get that uint s = gs/2 means that we split the array in half, but then it is a complete mystery. Without understanding it, I can't really implement another version for taking the maximum of an array for instance, even less for multi-dimensional arrays.

Furthermore, as far as I understand, the reduce kernel needs to be rerun another time if "N is bigger than the number of cores in a single unit".

Could you give me further explanations on that whole piece of code? Or even guidance on how to implement it for taking the max of an array?

Complete code can be found here: https://github.com/tmramalho/easy-pyopencl/blob/master/008_localreduce.py

take a look at the last diagram here github.com/mateuszbuda/GPUExample — Elad Maimoni

chapelo chapelo · Accepted Answer · 2021-01-16T00:53:29

Your first question about the meaning of the for loop:

for(uint s = gs/2; s > 0; s >>= 1)

It means that you divide the local size gs by 2, and keep dividing by 2 (the shift part s >>= 1 is equivalent to s = s/2) while s > 0, in other words, until s = 1. This algorithm depends on your array's size being a power of 2, otherwise you'd have to deal with the excess of a power of 2 until you have reduced the whole array, or you'd have to fill your array with neutral values for the reduction until completing a power of 2 size.

Your second concern when N is bigger than the capacity of your GPU, you are right: you have to run your reduction in portions that fit and then merge the results.

Finally, when you ask for guidance on how to implement a reduction to get the max of an array, I would suggest the following:

For a simple reduction like max or sum, try using numpy, especially if you are dealing with programming the reduction by axis.
If you think that the GPU would give you an advantage, try first using pyopencl's Multidimensional Array functionality, e.g. max.
If the reduction is more math intensive, try using pyopencl's Parallel Algorithms, e.g. reduction

I think that the whole point of using pyopencl is to avoid dealing with the underlying GPU's architecture. Otherwise, it is easier to deal with CUDA or HIP directly instead of OpenCL.

PyOpenCL - Multi-dimensional reduction kernel

1 Answers