Faster Parallel Reductions on Kepler

Originally published at:

Parallel reduction is a common building block for many parallel algorithms. A presentation from 2007 by Mark Harris provided a detailed strategy for implementing parallel reductions on GPUs, but this 6-year old document bears updating. In this post I will show you some features of the Kepler GPU architecture which make reductions even faster: the shuffle (SHFL)…

Nice post!

A "bit twiddling" variant of multiple reductions is to encode multiple values in a single word and then implicitly perform multiple reductions in one 32-bit SHFL reduction (or scan).

For example, if you categorized the values of a warp into 5 types (or 4 + "no match") and wanted to know the total for each type you can encode this into 5 6-bit fields and run your warp reduction (or scan). Each 6-bit field will total from 0-32.

I haven't tested if this kind of reduction is much faster than 5 ballot()+popc() ops but it's close to the same instruction count (10) in isolation. But I'm pretty sure it's faster than a prefix sum scan because of the extra instructions.

I've found this warp hack useful for quickly categorizing and binning of values!

Great information just the right time.
I was just planning to experiment with shuffle for our code.
I used the warp reduce to do some pre-reduction in the kernel that generates the data for the global sum in addition to the pre-reduction in registers. And as the generating kernel now has less stores to global memory it is slightly faster and the following reduction now only is over an array that is 32x smaller.

Another point: I was wondering about the return in 'warpReduceSumTriple' as well as the reinterpret_cast(&var); without template arguments. I have not tried, but I believe the compiler might not like that code.

Great catches. The warpReduceSumTriple error was a copy/paste mistake, and the missing reinterpret_cast template argument was due to the editor stripping '<' / '>' to avoid malicious html. I have fixed both errors!

In view of the increased ILP I was wondering whether you tried to use the interleaved warpReduceSum by default also for single reductions by dividing the array in multiple sections. I.e. split it into 4 parts and then use a warpReduceSumMultiple which does shuffle for int4.
Of course this will result in a smaller number of threads but eventually the increased ILP will yield a better performance.

Yes, see the discussion of vector loads Justin made in the post (and the Pro Tip on the subject). No need to divide the array into sections -- just cast the array to an array of int2 or int4 and load one vector per thread.

I was playing with merging the blockReduceSum and using it together with atomics directly in a Kernel which calculates the dot_product of two vectors.

When my problem size is not a multiple of the blockSize - which happen when tuning blockSizes - the last block will have less 'active' warps then blockDim/warpSize. I do bounds checking in the Kernel with the typical idx (=blockDim.x*blockIdx.x+threadIdx.x) > number_of_elements.
Of course I can activate the missing warps and just set their values to zero. This however resulted in a somewhat increased register usage and was slower. Is there a way to make blockReduceSum safe even if not all threads in the block are actually executing it, i.e., have exited because they belong to the last block and their index is beyond the problem size?

Typically I will keep those threads around. I
never put an early return in a kernel. If i’m doing a reduction across the
entire block I would structure the main loop using the grid stride loops
discussed in an earlier blog post and then reduce across the block after
exiting the loop. The for loop naturally handles the bounds checking you have
mentioned. At the end of the loop all threads are still present for the

In addition, the grid-stride loops also divides the array into multiple sections so that each thread processes multiple elements.

A misprint in section "Atomic Power":
if (threadIdx.x & warpSize == 0)
should be
if (threadIdx.x % warpSize == 0)

Fixed, thanks! Actually it's more efficient to write:

if (threadIdx.x & (warpSize - 1) == 0)

Because the compiler doesn't know that warpSize is a power of 2.

In the blockReduceSum() function above, is the keyword 'static' in front of the shared memory declaration really needed? nvcc throws an error on this in my setup...

Also isn't a bit-shift operation faster than the '/'? If so, why in the shuffle reduce loop do you do a warpSize/2 instead of (warpSize>>1) (or for that matter, an offset /=2 instead of offset >>=1)?

I was just realizing that I needed a fast single-kernel reduction. This allows me to easily launch up to 32 reductions, each one using a different cuda stream, so it fills the GPU utilization, and do the global sync just once at the end when all kernels have finished.


Need to discuss the shared memory bank conflict in the previous document by Mark Harris "reduction.pdf". Why there exist no shared memory bank conflicts in Reduction#1 (Interleaved Addressing) example whereas Reduction#2 suffers from shared memory bank conflict?

Reduction without storing array in shared memory is really great for my app!
It looks like this example of shuffle reduction works faster with this little tweak:
__inline__ __device__ int warpReduce(int val) {
#pragma unroll // for "16", not "warpSize/2"
for (int offset = 16; offset > 0; offset /= 2)
{ ... }
return val;

Yes, if you hard-code the warpSize it will definitely be faster. Unfortunately warpSize is not a compile-time constant.

static is optional as all shared memory allocations are static. I included it to be explicit. Older compilers did not support this keyword despite using the semantics of it. You are probably using an older compiler.

A shift is indeed faster but compilers should be able to make that transformation. I didn't use a shift because it makes the example a little more complicated.

Bank conflicts only occur between threads in the same warp. __shfl is handling all the exchanges within the warp and does not suffer from shared memory bank conflicts. Thus there is no reason to discuss bank conflicts in this example.

I think Munesh is referring to my "Optimizing Parallel Reduction in CUDA" document. The reason, Munesh, is that in reduction #1, (active) threads index shared memory with a stride of 1, while in #2 they index shared memory with a stride of increasing multiples of 2, meaning increasing number of accesses map to the same bank.

Thanks Mark, I got it. I sketched the mapping on a paper. Sorry had to put the question here as I could not find any forum elsewhere.