Multi-Pass Prefix Sum on Arbitrary Buffers in GLSL

This article explains how to compute a parallel prefix sum (scan) over arbitrarily sized buffers using multi-pass GLSL compute shader dispatches. While a single compute workgroup can efficiently scan data within its local memory limits, larger buffers require a hierarchical reduce-then-scan strategy. By breaking data into blocks, scanning intermediate block sums across consecutive dispatches, and propagating offsets back to the source elements, GPU compute pipelines can process inputs of any size while maintaining \(O(N)\) work efficiency and accommodating non-power-of-two lengths.

The Workgroup Limit Challenge

In modern graphics APIs like OpenGL or Vulkan, a single compute shader workgroup is bounded by hardware limits—typically 1024 invocations—and local shared memory capacity. Algorithms such as the Blelloch or Hillis-Steele scan can achieve high throughput inside a workgroup using shared memory and barrier() synchronization. However, hardware cannot synchronize across distinct workgroups within a single dispatch.

To process an arbitrary buffer size \(N\) that exceeds the capacity of a single workgroup, the computation must be split across multiple compute shader dispatches. The CPU host orchestrates these dispatches, using execution boundaries and global memory barriers to ensure correct data dependencies.

The Three-Pass Hierarchical Pipeline

A standard multi-pass scan relies on a block-wise divide-and-conquer strategy, generally requiring three distinct phases:

Pass 1: Local Scan & Reduction -> Pass 2: Block Sum Scan -> Pass 3: Uniform Add

Pass 1: Local Scan and Block Reduction

The primary buffer of size \(N\) is divided into uniform blocks of size \(B\) (where \(B\) is typically \(2 \times \text{workgroup size}\) to maximize per-thread throughput).

  • Each workgroup loads its assigned block from a Shader Storage Buffer Object (SSBO) into shared memory.
  • The workgroup performs an inclusive or exclusive scan locally.
  • The scanned values are written back to the primary buffer.
  • The total sum of each block is written to an auxiliary "block sums" buffer of size \(M = \lceil N / B \rceil\).

Pass 2: Recursive Scan of Block Sums

The auxiliary buffer containing \(M\) elements represents the aggregate sums of each block.

  • If \(M \le B\), a single workgroup dispatches to perform a local scan directly on this auxiliary buffer.
  • If \(M > B\), Pass 1 and Pass 2 are applied recursively: the auxiliary buffer is partitioned into blocks, its intermediate block sums are stored in a second auxiliary buffer, and so on, until the reduced size fits within a single workgroup.
  • Once the top of this hierarchy is reached and scanned, the scanned offsets bubble back down.

Pass 3: Uniform Add (Distribution)

Once the auxiliary buffer holds the prefix sums of all preceding blocks, a final dispatch adds each block’s base offset to its corresponding elements in the primary buffer.

  • Workgroup \(i\) loads the scalar prefix sum from index \(i\) of the auxiliary buffer.
  • Every thread in workgroup \(i\) adds this scalar value to its corresponding elements in the primary buffer.
  • Element \(0\) of the first block requires no offset (or \(0\)).

Handling Arbitrary Buffer Sizes

Arbitrary buffer sizes introduce elements that do not cleanly divide by the block size \(B\). To manage non-power-of-two sizes and arbitrary element counts without allocating padded memory buffers, the shader implementations utilize runtime boundary clamping:

  1. Grid Dispatch Calculation:
    The number of workgroups dispatched for a buffer of size \(N\) with a block capacity of \(B\) elements per workgroup is: \[\text{numWorkgroups} = \left\lfloor \frac{N + B - 1}{B} \right\rfloor\]

  2. Inactive Element Neutralization:
    When elements are loaded from the SSBO into shared memory, invocations check whether their global index exceeds \(N\):

    uint globalIdx = workgroupID * blockSize + localThreadID;
if (globalIdx < u_TotalElements) {
    sharedData[localThreadID] = inputData[globalIdx];
} else {
    sharedData[localThreadID] = 0; // Identity element for addition
}

    Padding inactive elements with the identity value (\(0\) for addition) ensures that out-of-bounds lanes do not corrupt the scan or the resulting block sums.

  3. Guarded Writes:
    When writing the locally scanned results back to the primary SSBO in Pass 1 and applying the base offsets in Pass 3, threads must check bounds again:

    if (globalIdx < u_TotalElements) {
    outputData[globalIdx] = sharedData[localThreadID] + blockOffset;
}

CPU-Side Synchronization

Because workgroups write to global SSBO memory in one phase that is required as input in the next, memory coherence must be enforced between dispatches.

In OpenGL, this requires invoking memory barriers on the host:

// Pass 1: Local Scan
glDispatchCompute(numWorkgroups, 1, 1);
glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT);

// Pass 2: Scan Block Sums
glDispatchCompute(auxWorkgroups, 1, 1);
glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT);

// Pass 3: Distribute Offsets
glDispatchCompute(numWorkgroups, 1, 1);
glMemoryBarrier(GL_SHADER_STORAGE_BARRIER_BIT);

In Vulkan, pipeline barriers with srcStageMask = VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT and dstStageMask = VK_PIPELINE_STAGE_COMPUTE_SHADER_BIT paired with matching VK_ACCESS_SHADER_WRITE_BIT to VK_ACCESS_SHADER_READ_BIT access masks are placed between each dispatch command within the command buffer.

By combining conditional data-loading with the three-pass hierarchical structure and host-side synchronization, GLSL compute pipelines can execute parallel prefix sums across multi-gigabyte buffers of arbitrary size while maintaining numerical correctness.