Optimizing reduction in HIP#
Reduction is a fundamental operation that uses a parallel pattern to combine a range of input values into a single scalar using a binary operation such as sum, maximum, or product. It appears throughout GPU workloads and is a practical vehicle for studying shared memory, thread divergence, bank conflicts, and memory bandwidth optimization.
This tutorial walks through a series of HIP kernels for summing a large array of floating-point values, with each step addressing a specific performance bottleneck. Starting from a naive interleaved-addressing kernel, each version eliminates a specific bottleneck: thread divergence, shared-memory bank conflicts, redundant synchronization barriers, and finally memory bandwidth.
The complete source file for all kernels is available at:
Note
All kernels target a large array of single-precision floating-point values
and are validated against a sequential reference sum. They compile with
amdclang++ -O3 -std=c++17 and run on any ROCm-supported GPU architecture.
Prerequisites#
Before starting this tutorial, ensure the following are in place.
ROCm installed and
amdclang++available onPATH.Familiarity with the HIP execution model (grids, blocks, wavefronts) and its mapping to AMD GPU hardware (dispatches, workgroups, wavefronts).
Application tracing and profiling using rocprofv3 installed for performance analysis.
Reduction fundamentals#
In functional programming, reduction is known as
fold. In the
C++ standard library, it appears as std::accumulate and, since C++17, as
std::reduce. A reduction takes a range of inputs, an identity element, and
a binary operation, and repeatedly applies the operation until a single value
remains. For addition, the identity is 0.
On a GPU, a parallel reduction takes a tree-shaped form. In each round, threads pair up and evaluate two values, reducing them to one, halving the number of live values until one remains.
Each block independently reduces its portion of the input to a single partial result. Repeating this process on successive output arrays until only one element remains completes the device-wide reduction without requiring global synchronization within a single kernel launch.
Compile and run:
amdclang++ -O3 -std=c++17 reduction.hip -o reduction
./reduction
Profile wall-clock time with rocprofv3:
rocprofv3 --kernel-trace --output-format csv -- ./reduction
The --kernel-trace CSV reports End_Timestamp - Start_Timestamp (both in
nanoseconds) for each kernel dispatch. Compare the duration column across
kernel variants to measure the impact of each optimization step.
Naive kernel#
The naive kernel loads one element per thread into shared memory and applies a
tree-shaped reduction using an interleaved addressing pattern. At each step,
only threads whose index is a multiple of 2 * stride are active and
participate in the reduction.
The following kernel implements this pattern:
This causes thread divergence. Within a wavefront, all lanes must execute the
same instruction at the same time. When some lanes take a branch, and others do
not, the hardware must execute both paths serially with the inactive lanes
masked off. In the interleaved pattern, at least one lane in each wavefront
hits the if statement at every level of the tree, so wavefronts remain
active long after the majority of their lanes have stopped doing useful work.
1__global__ void reduction_naive(const float *input, float *output, int n)
2{
3 extern __shared__ float sdata[];
4
5 int tid = threadIdx.x;
6 int gid = blockIdx.x * blockDim.x + tid;
7
8 sdata[tid] = (gid < n) ? input[gid] : 0.0f;
9 __syncthreads();
10
11 for (int stride = 1; stride < blockDim.x; stride *= 2)
12 {
13 if (tid % (2 * stride) == 0)
14 {
15 sdata[tid] += sdata[tid + stride];
16 }
17 __syncthreads();
18 }
19
20 if (tid == 0)
21 {
22 output[blockIdx.x] = sdata[0];
23 }
24}
What to observe#
Profile the naive kernel and record the following counter.
Counter |
What to look for |
|---|---|
Kernel duration (kernel-trace CSV) |
|
Reducing thread divergence#
Reduce divergence by reassigning which threads are active so that
inactive threads accumulate uniformly toward the upper end of the thread index
range. Once an entire wavefront is inactive, it can skip directly to
__syncthreads(), rather than executing the branch with all lanes masked.
This pattern, however, introduces a new problem: bank conflicts.
Resolving bank conflicts#
With AMD GPUs, shared memory (Local Data Share, or LDS) is organized into banks of 4 bytes each. On CDNA GPUs each Compute Unit has 32 banks. On RDNA GPUs the Work Group Processor has 64 banks, sub-divided into two sets of 32 banks each affiliated with a pair of SIMD32 units; a wavefront executes on one SIMD32 and maps its accesses to the affiliated 32-bank set. A bank conflict occurs when two or more threads in the same wavefront access different addresses that map to the same bank, causing those accesses to be serialized.
The reduced-divergence pattern still causes conflicts because the stride between active threads’ memory accesses doesn’t align with the bank layout. The fix is to set the stride to half the block size and halve each step, so that active threads always occupy the lower half of the remaining live range and access consecutive banks.
Note
To avoid bank conflicts, read and write shared memory in a coalesced manner, where each lane in a wavefront accesses a consecutive location. For more details, see the data share operations chapter of the CDNA3 ISA or RDNA3 ISA.
1__global__ void reduction_sequential(const float *input, float *output, int n)
2{
3 extern __shared__ float sdata[];
4
5 int tid = threadIdx.x;
6 int gid = blockIdx.x * blockDim.x + tid;
7
8 sdata[tid] = (gid < n) ? input[gid] : 0.0f;
9 __syncthreads();
10
11 for (int stride = blockDim.x / 2; stride > 0; stride >>= 1)
12 {
13 if (tid < stride)
14 {
15 sdata[tid] += sdata[tid + stride];
16 }
17 __syncthreads();
18 }
19
20 if (tid == 0)
21 {
22 output[blockIdx.x] = sdata[0];
23 }
24}
What to observe#
Compare the following counters against the naive kernel baseline.
Counter |
What to look for |
|---|---|
Kernel duration (kernel-trace CSV) |
|
|
Should be at or near zero, confirming that consecutive active threads access consecutive LDS banks. |
Wavefront reduction#
Every __syncthreads() operation in the reduction loop is necessary while threads from
different wavefronts are cooperating. Within a single wavefront, however,
threads execute in lockstep: they all advance through instructions together, so
a write by one lane is immediately visible to all other lanes in the same
wavefront without a barrier. Once the active thread count drops to one
wavefront, the remaining barriers are unnecessary overhead.
The wavefront reduction kernel exploits this by restructuring the algorithm into two phases. First, each wavefront independently reduces its own slice of shared memory without any barriers.
Lane 0 of each wavefront then writes its partial result into a compact staging area,
and a single __syncthreads() makes all partial results visible. The first
wavefront then reduces the staging area, again without barriers.
1template<int WarpSize>
2__device__ __forceinline__ void warp_reduce(float *sdata, int tid, int lane)
3{
4 if (WarpSize == 64) { if (lane < 32) { sdata[tid] += sdata[tid + 32]; } }
5 if (lane < 16) { sdata[tid] += sdata[tid + 16]; }
6 if (lane < 8) { sdata[tid] += sdata[tid + 8]; }
7 if (lane < 4) { sdata[tid] += sdata[tid + 4]; }
8 if (lane < 2) { sdata[tid] += sdata[tid + 2]; }
9 if (lane < 1) { sdata[tid] += sdata[tid + 1]; }
10}
1template<int WarpSize>
2__global__ __launch_bounds__(block_size)
3void reduction_warp(const float *input, float *output, int n)
4{
5 extern __shared__ float sdata[];
6
7 int tid = threadIdx.x;
8 int gid = blockIdx.x * blockDim.x + tid;
9 int lane = tid % WarpSize;
10 int warpid = tid / WarpSize;
11
12 sdata[tid] = (gid < n) ? input[gid] : 0.0f;
13
14 warp_reduce<WarpSize>(sdata, tid, lane);
15
16 int num_warps = blockDim.x / WarpSize;
17 if (lane == 0)
18 {
19 sdata[blockDim.x + warpid] = sdata[tid];
20 }
21 __syncthreads();
22
23 if (warpid == 0)
24 {
25 sdata[tid] = (lane < num_warps) ? sdata[blockDim.x + lane] : 0.0f;
26 warp_reduce<WarpSize>(sdata, tid, lane);
27 }
28
29 if (tid == 0)
30 {
31 output[blockIdx.x] = sdata[0];
32 }
33}
Note
The wavefront-level reduction shown here uses shared memory. On AMD GPUs, the same result can be achieved without shared memory traffic using shuffle instructions or Data Parallel Primitives (DPP), which exchange values between lanes entirely in registers. These techniques are covered in the Introduction to Compiler builtins chapter.
WGP mode and CU mode on RDNA GPUs#
For a full list of supported GPUs, see the ROCm system requirements.
On CDNA GPUs, the array is organized as a set of Compute Unit (CU) pipelines. Each CU contains four SIMD64 units and its own Local Data Share (LDS), which threads from wavefronts running on that CU can access. CDNA does not offer Work Group Processor mode as RDNA does, so the following information does not apply.
On RDNA GPUs, the array is organized as a set of Work Group Processor (WGP) pipelines. Each WGP contains two CUs, each with two SIMD32 units. The LDS is attached to the WGP, so threads from different wavefronts can access the same LDS if they run on CUs within the same WGP.
Wavefronts are dispatched in one of two modes. These control whether wavefronts are distributed across two SIMD32s within a single CU (CU mode) or across all four SIMD32s within a WGP (WGP mode).
CU mode executes two wavefronts per block on a CU and provides only half the LDS to each wavefront. Independence between CUs can improve performance for workloads that avoid inter-wavefront communication.
WGP mode executes four wavefronts per block on a WGP with a shared LDS. It can increase occupancy and improve performance for workloads without heavy inter-wavefront communication, but it can degrade performance for programs that rely on atomics or extensive inter-wavefront communication through shared memory.
The wavefront reduction kernel communicates partial results between wavefronts through
shared memory, making it sensitive to this distinction. The inter-wavefront staging
step benefits from CU mode because all wavefronts in the block are guaranteed to
share the same LDS instance. Compile with -mcumode to enable CU mode on
RDNA GPUs. Memory-bandwidth-bound kernels such as the vectorized loads
kernel are unaffected by this setting.
What to observe#
Compare the following counters against the sequential-addressing kernel.
Counter |
What to look for |
|---|---|
Kernel duration (kernel-trace CSV) |
|
|
Reduction in LDS stall cycles (fewer barriers mean less time waiting for shared memory to become consistent). |
Vectorized loads#
All kernels so far issue one 32-bit load per thread. The GPU memory system can
issue 128-bit loads at the same cost, so replacing four scalar loads with a
single float4 read quadruples the data moved per instruction. Each thread
loads four consecutive elements, reduces them to a scalar in registers, and
then enters the same wavefront reduction as before.
1template<int WarpSize>
2__global__ __launch_bounds__(block_size)
3void reduction_float4(const float *input, float *output, int n)
4{
5 extern __shared__ float sdata[];
6
7 int tid = threadIdx.x;
8 int gid = blockIdx.x * (blockDim.x * 4) + tid;
9 int lane = tid % WarpSize;
10 int warpid = tid / WarpSize;
11
12 float val = 0.0f;
13 if (gid + 3 * blockDim.x < n)
14 {
15 float4 v = reinterpret_cast<const float4 *>(input)[blockIdx.x * blockDim.x + tid];
16 val = v.x + v.y + v.z + v.w;
17 }
18 else
19 {
20 if (gid < n) { val += input[gid]; }
21 if (gid + blockDim.x < n) { val += input[gid + blockDim.x]; }
22 if (gid + 2 * blockDim.x < n) { val += input[gid + 2 * blockDim.x]; }
23 if (gid + 3 * blockDim.x < n) { val += input[gid + 3 * blockDim.x]; }
24 }
25
26 int num_warps = blockDim.x / WarpSize;
27 sdata[tid] = val;
28 warp_reduce<WarpSize>(sdata, tid, lane);
29
30 if (lane == 0)
31 {
32 sdata[blockDim.x + warpid] = sdata[tid];
33 }
34 __syncthreads();
35
36 if (warpid == 0)
37 {
38 sdata[tid] = (lane < num_warps) ? sdata[blockDim.x + lane] : 0.0f;
39 warp_reduce<WarpSize>(sdata, tid, lane);
40 }
41
42 if (tid == 0)
43 {
44 output[blockIdx.x] = sdata[0];
45 }
46}
Note
The float4 reinterpret cast requires the input pointer to be 16-byte
aligned. Allocations from hipMalloc satisfy this requirement.
What to observe#
Compare the following counters against the wavefront reduction kernel.
Counter |
What to look for |
|---|---|
Kernel duration (kernel-trace CSV) |
|
|
Reduction in VMEM instruction cycles (fewer instructions for the same total data moved). |
Device-level reduction#
The block kernels above each produce one partial result per block. To reduce an entire array to a scalar, apply the block kernel repeatedly on successively smaller arrays, alternating between two scratch buffers, until only one value remains. For an input of N elements processed B elements per block, this completes in logB(N) passes.
1 {
2 float *d_partial_a = nullptr;
3 float *d_partial_b = nullptr;
4 HIP_CHECK(hipMalloc(&d_partial_a, sizeof(float) * nb1));
5 HIP_CHECK(hipMalloc(&d_partial_b, sizeof(float) * nb1));
6
7 int current = input_size;
8 const float *src = d_input;
9 float *dst = d_partial_a;
10
11 while (current > 1)
12 {
13 int nb = (current + block_size * 4 - 1) / (block_size * 4);
14 if (warpsizehost == 64)
15 {
16 reduction_float4<64><<<nb, block_size, smem_warp>>>(src, dst, current);
17 }
18 else
19 {
20 reduction_float4<32><<<nb, block_size, smem_warp>>>(src, dst, current);
21 }
22 HIP_CHECK(hipGetLastError());
23 current = nb;
24 if (current > 1)
25 {
26 src = dst;
27 dst = (dst == d_partial_a) ? d_partial_b : d_partial_a;
28 }
29 }
30
31 HIP_CHECK(hipDeviceSynchronize());
32 float result = 0.0f;
33 HIP_CHECK(hipMemcpy(&result, dst, sizeof(float), hipMemcpyDeviceToHost));
34
35 // Compute the scalar reference by reducing the per-block CPU partials.
36 float ref_scalar = 0.0f;
37 for (float v : ref_bs4)
38 {
39 ref_scalar += v;
40 }
41 float rel_err = std::fabs(result - ref_scalar) / (std::fabs(ref_scalar) + 1.0f);
42
43 std::cout << " Device reduction: ";
44 std::cout << (rel_err < 1e-3f ? "PASSED" : "FAILED") << std::endl;
45
46 HIP_CHECK(hipFree(d_partial_a));
47 HIP_CHECK(hipFree(d_partial_b));
48 }
Further reading#
The following resources provide deeper coverage of the tools and hardware referenced in this tutorial.
rocPRIM - production-quality reduction primitives that handle edge cases and automatically apply architecture-specific tuning.
rocprofv3 documentation - detailed guide to timeline and counter profiling.
AMD GPU architecture guides (ISA references) - VGPR budgets, LDS bank geometry, and wavefront scheduling details for each architecture family.