Table of Contents
Reduction (CUDA)
Summing an array on a GPU has no direct equivalent of OpenMP's reduction clause. Reduction in CUDA is written by hand as a tree: each block reduces its own slice in shared memory, then the partial results are combined.
__global__ void reduce(const float *in, float *out, int n) { __shared__ float s[256]; int t = threadIdx.x; int i = blockIdx.x * blockDim.x + t; s[t] = (i < n) ? in[i] : 0.0f; __syncthreads(); for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) { if (t < stride) s[t] += s[t + stride]; __syncthreads(); } if (t == 0) out[blockIdx.x] = s[0]; }
Each iteration halves the number of active threads, so a block of 256 finishes in 8 steps instead of 256. The __syncthreads() inside the loop is required, and it sits outside the if so that every thread reaches it (Synchronization (CUDA)).
The stride starts large and halves. Doing it the other way round, with neighbouring threads accumulating into adjacent slots, produces bank conflicts and leaves the active threads scattered across warps instead of packed into the first few.
The kernel emits one partial per block, so it does not finish the job. Two options: launch it again on the partials until one value remains, or have each block finish with a single atomicAdd into a global accumulator. The second is simpler and usually fine, since the contention is one atomic per block rather than one per element.
Once stride drops below 32, all the surviving threads are in one warp and the shared memory round trip can be replaced by warp shuffles:
for (int off = 16; off > 0; off >>= 1) v += __shfl_down_sync(0xffffffff, v, off);
In production, use CUB's BlockReduce or Thrust's reduce rather than hand-rolling any of this. The hand-written version is worth understanding because the same tree pattern shows up in scan, sort, and histogram kernels.
