GPU Parallel Reduction with Shared Memory

1. The Parallel Reduction Strategy

Imagine you need to sum a very large array of numbers. In a sequential program, you would use a simple loop. In parallel, we can do better. The standard approach is a tree-based reduction. Threads work in pairs to sum elements, then half of those threads sum the results, and so on, until a single sum remains.

A conceptual view of parallel reduction. An array of N elements is reduced to N/2 partial sums, then N/4, and so on, until a single result is produced.

To implement this efficiently on a GPU, we'll use shared memory to hold the intermediate sums for each block. Let's walk through the process step-by-step, starting with a basic implementation and progressively optimizing it.

The First Bottlenecks: Divergence and Bank Conflicts

As the article points out, the reduce0 kernel has two major flaws.

1. Warp Divergence: The condition if (tid % (2*s) == 0) causes threads within the same warp to take different execution paths, leading to serialization and poor performance.
2. Shared Memory Bank Conflicts: Even after fixing the divergence (as shown in the article's REDUCE-1), the memory access pattern itself is problematic. This is called interleaved addressing, where threads access memory locations that are far apart. This often leads to multiple threads in a warp trying to access different addresses within the same shared memory bank, causing the accesses to be serialized.

2. A Better Way: Sequential Addressing

To solve the bank conflict issue, we need to change our memory access pattern. Instead of threads accessing distant elements, we'll have them access adjacent elements. This is called sequential addressing.

The key change in REDUCE-2 is the loop structure:

// Previous, interleaved approach
for(unsigned int s = 1; s < blockDim.x; s *= 2) {
    int index = 2 * s * tid;
    if (index < blockDim.x) {
        sdata[index] += sdata[index + s];
    }
    __syncthreads();
}

// New, sequential approach
for(unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
    if (tid < s) {
        sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
}

In the sequential approach, the stride s starts large and halves each time. Active threads (tid < s) add the value from the second half of the active array sdata[tid + s] to their own value sdata[tid]. This ensures that all memory accesses within a warp are to contiguous memory locations, avoiding bank conflicts.

3. Advanced Optimizations: Putting All Threads to Work

The sequential addressing approach is much better, but it's still not perfect. The condition if (tid < s) means that in the very first iteration, half of the threads in the block are idle. We can do better.

Making Use of Idle Threads

Instead of each thread loading one element, we can have each thread load two elements and perform the first step of the reduction right away. This simple change halves the amount of shared memory needed and immediately engages all threads in useful work.

This technique is described in the REDUCE-3 section of the article you are reading. The core change is in the data loading part of the kernel:

// Before: one load per thread
sdata[tid] = g_in_data[i];
__syncthreads();

// After: two loads and one add per thread
unsigned int i = blockIdx.x*(blockDim.x*2) + threadIdx.x;
sdata[tid] = g_in_data[i] + g_in_data[i + blockDim.x];
__syncthreads();

Now, all threads contribute to the first level of the reduction tree, making better use of the GPU's resources.

Unrolling the Final Warp

As the reduction proceeds, the number of active threads decreases. When the stride s becomes less than or equal to 32, all remaining active threads are within a single warp.

Recall from our previous lesson that threads within a warp execute in lockstep (are SIMD synchronous). This has a powerful implication: for operations entirely within a single warp, you do not need __syncthreads(). The if (tid < s) check also becomes redundant because all threads in the warp execute the instruction anyway. We can exploit this by "unrolling" the last 6 iterations of the loop (for strides 32, 16, 8, 4, 2, and 1) into a series of explicit additions. This eliminates loop overhead and unnecessary synchronization calls, providing a significant performance boost.

This is a more advanced technique, so let's watch a video that explains it clearly.

4. The Optimized Reduction Kernel

By combining these optimizations, we arrive at a highly efficient parallel reduction kernel. Let's look at the structure that brings these ideas together. This corresponds to a combination of REDUCE-3 and REDUCE-4 from the article.

// Helper function for the final warp-level reduction
// 'volatile' is crucial to ensure writes to shared memory are not optimized away by the compiler
__device__ void warpReduce(volatile int* sdata, int tid) {
    sdata[tid] += sdata[tid + 32];
    sdata[tid] += sdata[tid + 16];
    sdata[tid] += sdata[tid + 8];
    sdata[tid] += sdata[tid + 4];
    sdata[tid] += sdata[tid + 2];
    sdata[tid] += sdata[tid + 1];
}

// Main reduction kernel
__global__ void reduce_kernel(int *g_in_data, int *g_out_data) {
    extern __shared__ int sdata[];

    unsigned int tid = threadIdx.x;
    // Each thread performs the first reduction step during the load from global memory
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    sdata[tid] = g_in_data[i] + g_in_data[i + blockDim.x];
    __syncthreads();

    // Do parallel reduction in shared memory, stopping before the last warp
    for (unsigned int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }

    // Unrolled reduction for the final warp
    if (tid < 32) {
        warpReduce(sdata, tid);
    }

    // Thread 0 writes the block's final result back to global memory
    if (tid == 0) {
        g_out_data[blockIdx.x] = sdata[0];
    }
}

This kernel represents a sophisticated implementation that is far more performant than our initial naive attempt. It effectively manages memory access patterns, maximizes thread utilization, and reduces instruction overhead by understanding the fundamental properties of the CUDA execution model.

Conclusion

This lesson concludes our module on CUDA programming. By walking through the implementation and optimization of a parallel reduction, you have seen how all the core concepts—kernels, thread hierarchy, memory spaces, and synchronization—come together to solve a real problem.

Key Takeaways:

• Parallel Reduction is a fundamental tree-based algorithm for aggregating data in parallel.
• Shared Memory is essential for holding intermediate results within a block, avoiding costly global memory traffic.
• Optimization is a Process: Performance is achieved by methodically identifying and eliminating bottlenecks, including:
  • Bank Conflicts: Solved by using sequential addressing.
  • Idle Threads: Mitigated by performing more work per thread, such as the first add during load.
  • Instruction Overhead: Reduced by unrolling loops, especially for the final warp where __syncthreads() is not needed.
• Understanding the hardware's behavior (e.g., warp-synchronous execution) is key to unlocking advanced optimizations.

You now have a solid foundation in the principles of GPU programming. In our next module, "Infrastructure for Scalable Systems," we will broaden our perspective, moving from the architecture of a single node to the principles behind building large-scale distributed and cloud computing systems. We will begin by exploring the components of a commodity cluster.