Blelloch Exclusive Scan in GLSL Compute Shaders

The Blelloch exclusive prefix scan is a work-efficient parallel algorithm designed to compute prefix sums with \(O(N)\) total additions across \(O(\log N)\) parallel steps. This guide details how to implement the Blelloch algorithm inside a single GLSL compute shader workgroup using shared memory (shared) and workgroup execution barriers (barrier()). By splitting execution into an up-sweep (reduction) phase and a down-sweep phase, threads within a workgroup cooperatively calculate exclusive running sums entirely on the GPU.

Algorithm Phases

The Blelloch scan operates on an array of size \(N\) (where \(N\) is a power of two) mapped across \(N / 2\) invocation threads:

1. Up-Sweep (Reduce Phase): Traverses the array from leaves to root, building a binary reduction tree in-place. At each step, partial sums are stored at specific indices.
2. Clear the Root: The last element of the array (\(N - 1\)) is explicitly set to zero, serving as the identity element for exclusive prefix sums.
3. Down-Sweep Phase: Traverses back down the tree from the root to the leaves. At each step, each node passes its current value to its left child and updates its right child with the sum of its old value and the left child's value.

GLSL Implementation

Below is a complete GLSL compute shader function demonstrating an exclusive Blelloch scan on an array of size SCAN_SIZE = 1024 using a local workgroup size of 512 threads.

#version 450

#define WORKGROUP_SIZE 512
#define SCAN_SIZE (WORKGROUP_SIZE * 2)

layout(local_size_x = WORKGROUP_SIZE) in;

// Workgroup shared memory for local scanning
shared uint shared_data[SCAN_SIZE];

void blelloch_exclusive_scan(uint local_id) {
    // -------------------------------------------------------------
    // Phase 1: Up-Sweep (Reduction)
    // -------------------------------------------------------------
    uint stride = 1;
    for (uint d = SCAN_SIZE >> 1; d > 0; d >>= 1) {
        barrier();
        if (local_id < d) {
            uint ai = stride * (2 * local_id + 1) - 1;
            uint bi = stride * (2 * local_id + 2) - 1;
            shared_data[bi] += shared_data[ai];
        }
        stride <<= 1;
    }

    // -------------------------------------------------------------
    // Phase 2: Clear Root
    // -------------------------------------------------------------
    if (local_id == 0) {
        shared_data[SCAN_SIZE - 1] = 0;
    }

    // -------------------------------------------------------------
    // Phase 3: Down-Sweep
    // -------------------------------------------------------------
    for (uint d = 1; d < SCAN_SIZE; d <<= 1) {
        stride >>= 1;
        barrier();
        if (local_id < d) {
            uint ai = stride * (2 * local_id + 1) - 1;
            uint bi = stride * (2 * local_id + 2) - 1;

            uint temp = shared_data[ai];
            shared_data[ai] = shared_data[bi];
            shared_data[bi] += temp;
        }
    }
    barrier();
}

void main() {
    uint local_id = gl_LocalInvocationID.x;

    // Load elements into shared memory (each thread loads 2 elements)
    shared_data[2 * local_id]     = 1; // Replace with input buffer value
    shared_data[2 * local_id + 1] = 1; // Replace with input buffer value

    // Perform scan in-place
    blelloch_exclusive_scan(local_id);

    // shared_data now holds the exclusive prefix sum of the elements
}

Execution Details

• Index Calculation: The indices ai and bi target the left and right children at each level of the virtual binary tree.
• Barriers: Calls to barrier() ensure that all active threads read the previous iteration's results before any thread writes to shared memory for the current step.
• Exclusive vs. Inclusive: The algorithm produces an exclusive scan because the last element is set to 0 before the down-sweep. If an inclusive scan is required, the original array elements can be added back, or each thread can shift results left by one position and append the reduction sum at the end.