CUDA Gather and Scatter

Defining Gather and Scatter

At their core, gather and scatter are about indirect addressing. They use an auxiliary array of indices to determine where to read data from (gather) or write data to (scatter).

Let's start with a clear, code-based definition.

As you can see from the reading:

• Gather: An output element at a contiguous index i is assigned a value from a potentially non-contiguous location idx[i] in an input array. The indirection happens on the right-hand side (the read operation).
  output[i] = input[map[i]]
• Scatter: An input element from a contiguous index i is written to a potentially non-contiguous location map[i] in an output array. The indirection happens on the left-hand side (the write operation).
  output[map[i]] = input[i]

The following diagram provides a helpful visual for these operations. The Index array corresponds to the idx or map array in our pseudocode.

Another useful way to think about this, as mentioned in the forum post, is in terms of memory access patterns:

• Gather transforms data from non-contiguous storage into a contiguous block.
• Scatter transforms data from a contiguous block into non-contiguous storage.

Performance and Implementation in CUDA

While conceptually simple, these patterns have significant performance implications, especially on a throughput-oriented architecture like a GPU. Since you are focused on low-level performance, understanding these implications is crucial.

The primary issue is memory coalescing, which we covered in Module 5. When threads in a warp access memory, performance is highest when their accesses can be "coalesced" into a single, wide transaction. This happens when they access consecutive memory locations.

Gather Kernel

Let's write a simple CUDA kernel for the gather operation. Each thread will be responsible for computing one element of the output array.

__global__ void gather_kernel(float* output, const float* input, const int* map, int N) {
    // Calculate the global thread ID
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid < N) {
        // Each thread writes to its own contiguous location in 'output'.
        // It reads from a location in 'input' dictated by the map.
        output[tid] = input[map[tid]];
    }
}

Let's analyze the memory accesses in this kernel from the perspective of a single warp:

• Writes to output: The write output[tid] = ... is perfectly coalesced, as consecutive threads write to consecutive memory locations.
• Reads from map: The read map[tid] is also perfectly coalesced.
• Reads from input: The read input[map[tid]] is the critical part. Since the map array can contain arbitrary indices, the access pattern into input is unpredictable. It is highly likely to be uncoalesced, leading to multiple memory transactions and reduced memory bandwidth.

This is why gather is sometimes referred to as a potentially "bad pattern" if performance is paramount and the map is random. However, it is often unavoidable.

Scatter Kernel and Race Conditions

Now for the scatter kernel.

__global__ void scatter_kernel(float* output, const float* input, const int* map, int N) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid < N) {
        // Each thread reads from its own contiguous location in 'input'.
        // It writes to a location in 'output' dictated by the map.
        output[map[tid]] = input[tid];
    }
}

The memory access analysis is the inverse of gather:

• Reads from input and map: These are coalesced.
• Writes to output: The write output[map[tid]] = ... is likely uncoalesced, just as the read was in the gather kernel.

However, scatter has an additional, more serious correctness issue: race conditions. What if the map is not a permutation? For example, what if map[5] and map[10] both equal 42? In this case, thread 5 and thread 10 would attempt to write to output[42] simultaneously. The final value at that location would be indeterminate.

For a simple scatter like the one above, we must assume the map is a permutation, meaning every index appears exactly once. If duplicates are possible, we must use atomic operations (e.g., atomicAdd, atomicExch) to ensure correctness. Using atomics turns the scatter into a different primitive, resembling the histogram we implemented earlier.

Application: Stream Compaction (Parallel Filter)

Gather and scatter are not just standalone operations; they are powerful building blocks. A classic example is stream compaction, also known as filtering. The goal is to take a large array and create a new, smaller array containing only the elements that satisfy a certain condition (a predicate).

Imagine you have an array [1, 8, 3, 6, 5, 4] and you want to keep only the even numbers. The desired output is [8, 6, 4]. This can be done efficiently in parallel by combining map, scan, and scatter.

The process involves three steps:

1. Map (Predicate Evaluation): Create a temporary flags array of 0s and 1s. For each element in the input, a thread writes a 1 to the flags array if the element should be kept, and a 0 otherwise.

   • Input: [1, 8, 3, 6, 5, 4]
   • Flags: [0, 1, 0, 1, 0, 1]

2. Exclusive Scan: Perform an exclusive scan on the flags array. The resulting offsets array tells us the destination index for each element that we are keeping. The total number of elements to keep is the last value of the inclusive scan, or the sum of all flags.

   • Flags: [0, 1, 0, 1, 0, 1]
   • Offsets (Exclusive Scan Result): [0, 0, 1, 1, 2, 2]

3. Scatter (Conditional): A kernel iterates from i = 0 to N-1. If flags[i] is 1, it scatters input[i] to output[offsets[i]].

   • input[1] (which is 8) goes to output[offsets[1]] (which is output[0]).
   • input[3] (which is 6) goes to output[offsets[3]] (which is output[1]).
   • input[5] (which is 4) goes to output[offsets[5]] (which is output[2]).
   • Final Output: [8, 6, 4]

This elegant composition of primitives allows a complex, data-dependent reorganization to be performed in a structured, parallel way.

The final scatter step of the compaction algorithm can be implemented with the following kernel. Notice the conditional logic before the write operation.

__global__ void compact_scatter_kernel(
    float* output, 
    const float* input, 
    const int* flags, 
    const int* offsets, 
    int N) 
{
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    if (tid < N) {
        // Only scatter the elements we decided to keep.
        if (flags[tid] == 1) {
            // The scan result gives us the correct destination index.
            int destination_idx = offsets[tid];
            output[destination_idx] = input[tid];
        }
    }
}

Conclusion

In this lesson, we have added the final fundamental primitives, gather and scatter, to our parallel programming toolkit. You now understand how to move data according to arbitrary index maps and how these operations form the basis of more complex algorithms.

Key Takeaways:

• Gather is an indirect read (output[i] = input[map[i]]). Its performance is often limited by uncoalesced reads from the input array.
• Scatter is an indirect write (output[map[i]] = input[i]). Its performance can be limited by uncoalesced writes, and more importantly, it introduces a potential for race conditions if the destination indices are not unique.
• Correctly handling non-unique indices in a scatter requires atomic operations, which changes its behavior to be more like a histogram.
• Gather and scatter are powerful building blocks. We saw how they can be combined with scan to implement a highly efficient parallel stream compaction algorithm.

With a solid grasp of reduction, scan, histogramming, gather, and scatter, you are now equipped with the core components used in a vast range of high-performance algorithms.

This lesson concludes our module on parallel primitives. We have focused on what to compute and how to structure the computation. In our next module, we will shift our focus to performance analysis. We will begin by learning how to measure a kernel's efficiency, classify it as memory-bound or compute-bound, and use professional profiling tools to scientifically diagnose and fix bottlenecks.