100 DAYS OF GPU Course preview
Day 31Module 4
in-technical-review

CUDA prefix sum, the Hillis-Steele scan

The parallel prefix sum most people meet first is GPU Gems 3, chapter 39, and its opening algorithm is a loop you can write straight from the description:

for (unsigned int offset = 1; offset < blockDim.x; offset *= 2) {
    if (tid >= offset) {
        tile[tid] += tile[tid - offset];
    }
    __syncthreads();
}

The chapter says outright that this does not work: "Because not all threads run simultaneously for arrays larger than the warp size, Algorithm 1 will not work, because it performs the scan in place on the array."

Right about the outcome, dated about the cause, and quiet about the failure. Day 14 measured a kernel with this shape of race and got zero wrong answers at 32 threads a block and 17,312 at 64. Your small test passes and you ship the bug.

By the end of this page you can write a block-wide inclusive scan and an exclusive one, and say which of the two owes the next kernel a number.

Two scans, and which one hands you the total

A scan is a running total. Given 3 1 4 1, the inclusive scan is 3 4 8 9 and the exclusive scan is 0 3 4 8.

GPU Gems puts the difference in one sentence: in an exclusive scan "each element j of the result is the sum of all elements up to but not including j in the input array. In an inclusive scan, all elements including j are summed."

That is not bookkeeping. Stream compaction on day 33 asks, for every element it keeps, how many kept elements sit to its left; that number is the element's slot in the output and it must not count the element itself, so it is an exclusive scan. Day 35 builds a radix sort on the same primitive for the same reason.

The second difference costs more time. An inclusive scan ends holding the total; an exclusive scan does not. The last element of an inclusive scan over a tile is the sum of the tile.

The last element of an exclusive scan is the sum of everything but the last element, so getting the total back costs another read and another add. Day 32 wants one total per block, so which of the two you build decides whether that number is already available.

Hillis-Steele itself is one idea repeated. At the first step every element adds the element one place to its left, at the second the element two places to its left, then four, then eight. After the step at offset d an element holds the sum of the 2 * d elements ending at it, or of everything up to it where there are fewer, so log2(n) steps finish a tile of n and a block of 256 threads scanning its own tile takes eight.

Eight steps, and 1,793 additions against a sequential scan's 255. The general form is n log2(n) - (n - 1) against n - 1.

GPU Gems states the cost: "The algorithm performs O(n log2 n) addition operations. Remember that a sequential scan performs O(n) adds. Therefore, this naive implementation is not work-efficient."

Hillis-Steele is short, it is the easiest correct parallel scan to write, and day 32 replaces it with one that does linear work.

Why one array is not enough

Look again at tile[tid] += tile[tid - offset]. Thread tid reads the word belonging to thread tid - offset in the same step in which that thread is writing it, and the only __syncthreads() sits below both.

Whether you get the old value or the new one depends on which thread got there first. That is a data race, and day 14 is the lesson about the barrier that removes it.

GPU Gems blames warp scheduling, and in 2007 that explained the failure: inside one warp the 32 lanes issued together, so they all read before any of them wrote. Volta ended that rule.

With independent thread scheduling every lane carries its own program counter and nothing obliges two lanes to sit on the same instruction. A kernel that is right only because one warp happens to stay together is not right, it is untested, and day 14's zero at 32 threads is what untested looks like.

The fix is to stop writing where you are reading. Give the block two halves of a shared array, read the half holding the current state, write the other, and swap. One barrier a step is enough, because it separates every read of a half from the next write to that half.

Double buffering then creates a bug the in-place version cannot have. In place, a thread with nothing to add is correct to do nothing: its value is already in the array.

Across two halves, doing nothing leaves the destination holding whatever was there two steps ago, and that stale value is what the next step reads. Every thread has to write every step, including the ones with no addition to perform.

What the program checks, and what it refuses to

Full program in code/day31-scan-1/scan.cu. It runs four kernels at four block sizes, 32, 64, 256 and 1024, over the same input.

Three rules hold it together.

Every element is (i % 8) + 1, so the check is exact. The largest prefix a tile can reach is 4,608, a whole number well under 2^24 that a float holds exactly, and the element count leaves a partial tile at every block size in the sweep. The comparison is != with no tolerance at all: a mismatch means the wrong set of elements was added, never a rounding difference.

The two broken kernels are reported, not asserted. The program prints how many elements each got wrong and never gates on the answer being wrong. One of them fails through a race, and nothing promises a race the same outcome twice.

Nothing is timed. A block scan of 256 elements is smaller than the launch that carries it, so a clock here would measure launch overhead, which is day 9's subject.

The in-place version is in the file, so you can run it rather than imagine it:

__global__ void scanInclusiveInPlace(const float* __restrict__ in,
                                     float* __restrict__ out, size_t n) {
    __shared__ float tile[kMaxThreadsPerBlock];

    const unsigned int tid = threadIdx.x;
    const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;

    // The guard covers the load and the store, never the barrier. Zero is the
    // identity for addition, so a thread past the end of the array carries a
    // value that changes nobody's prefix.
    tile[tid] = (i < n) ? in[i] : 0.0f;
    __syncthreads();

    for (unsigned int offset = 1; offset < blockDim.x; offset *= 2) {
        if (tid >= offset) {
            tile[tid] += tile[tid - offset];
        }
        __syncthreads();
    }

    if (i < n) {
        out[i] = tile[tid];
    }
}

And the fix. src always names the half holding the complete state, so it names the answer when the loop ends and there is no final flip to get wrong:

__global__ void scanInclusiveDoubleBuffer(const float* __restrict__ in,
                                          float* __restrict__ out, size_t n) {
    __shared__ float tile[2 * kMaxThreadsPerBlock];

    const unsigned int tid = threadIdx.x;
    const unsigned int width = blockDim.x;
    const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + tid;

    tile[tid] = (i < n) ? in[i] : 0.0f;
    __syncthreads();

    unsigned int src = 0;
    for (unsigned int offset = 1; offset < width; offset *= 2) {
        const unsigned int dst = 1 - src;
        if (tid >= offset) {
            tile[dst * width + tid] =
                tile[src * width + tid] + tile[src * width + tid - offset];
        } else {
            tile[dst * width + tid] = tile[src * width + tid];
        }
        __syncthreads();
        src = dst;
    }

    if (i < n) {
        out[i] = tile[src * width + tid];
    }
}

That else is the whole of the second kernel in the file, which drops it and gets the answer wrong at every block size. Nothing else here is subtle: the load from global memory is a coalesced run of 32 consecutive floats, and the shared reads are words tid and tid - offset, consecutive across a warp and so spread over 32 distinct banks. There is no bank conflict in this kernel and day 15's arithmetic has nothing to say about it.

The exclusive scan is that same kernel with one changed store. After the loop, instead of writing its own inclusive prefix, each thread writes its left neighbour's, and thread 0 writes the zero the exclusive scan starts with:

    if (i < n) {
        out[i] = (tid > 0) ? tile[src * width + tid - 1] : 0.0f;
    }

Note. That last read needs no barrier of its own. The final iteration of the loop ended in a __syncthreads(), nothing writes the tile after it, and a shared read that races with no write is not a race.

One honest limit: every tile holds the same repeating pattern, so the check catches a dropped, doubled or misplaced element and would miss a kernel that is wrong only on some other arrangement of values.

Results

Measured. Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), built with nvcc -O3 -arch=sm_75. Captured 2026-09-01 on the project's verification node; full transcript in the page's evidence file.

GPU: Tesla T4 (compute capability 7.5), warp size 32
Max threads per block: 1024, shared memory per block: 49152 B

Work per block, computed at compile time
 threads    steps    scan adds   sequential
 -------    -----    ---------   ----------
      32        5          129           31
      64        6          321           63
     256        8         1793          255
    1024       10         9217         1023

1049187 elements, in[i] = (i % 8) + 1, one tile per block
threads  kernel                       mismatches   first bad  gated
-------  ---------------------------  ----------   ---------  -----
     32  inclusive, in place                   0           -     no
     32  inclusive, no copy               622954           2     no
     32  inclusive, double buffered            0           -    yes
     32  exclusive, double buffered            0           -    yes
     64  inclusive, in place               70320        5600     no
     64  inclusive, no copy               704924           1     no
     64  inclusive, double buffered            0           -    yes
     64  exclusive, double buffered            0           -    yes
    256  inclusive, in place              170530        3488     no
    256  inclusive, no copy               823776           1     no
    256  inclusive, double buffered            0           -    yes
    256  exclusive, double buffered            0           -    yes
   1024  inclusive, in place                8819         864     no
   1024  inclusive, no copy               901653           1     no
   1024  inclusive, double buffered            0           -    yes
   1024  exclusive, double buffered            0           -    yes

both double-buffered kernels matched the reference at every block size

The race is real, and it hides at one warp. The in-place kernel reads data[i - stride] and writes data[i] in the same step, so without two buffers some threads read a value their neighbour has already overwritten. Its row at 32 threads has zero mismatches.

At 64 it gets 70,320 elements wrong, at 256 it gets 170,530 wrong, 16 percent of the array, and at 1024 it is back down to 8,819. The zero at one warp is the risk because a kernel that passes every test at 32 threads can still be wrong.

This is the race from day 14 inside a scan.

The missing copy fails sooner and more clearly. The no-copy kernel has both buffers and skips the else branch, so idle threads never carry their value forward, and it gets 622,954 elements wrong at 32 threads, first failing at index 2.

That is 59 percent of the array at every block size, starting near the front. A bug that loud gets caught. The race above is easier to miss.

Hillis-Steele does more work than a sequential scan, by design. The compile-time table says a 256-thread tile costs 1793 adds against 255 done sequentially, which is seven times the arithmetic. It buys parallelism with work.

A GPU can run those adds across many lanes. Day 32 shows the version that stops paying this cost.

Run it yourself

Use a CUDA GPU that supports the lesson's minimum compute capability. The build line is in the repo's README:

nvcc -std=c++17 -O3 -arch=sm_75 -o scan scan.cu

No Compiler Explorer embed. Not because of shared memory, which is 8 KiB per block, but because the program scans 1,049,187 elements sixteen times on the device and once per block size on the host, which sits too close to Compiler Explorer's 20 second run cap to pin a lesson to. Cut the element count to a few thousand and one block size and it fits.

Exercise

Build the exclusive scan the other way round: shift the input right by one before the scan instead of shifting the result after it, so thread tid loads element tid - 1 and thread 0 loads zero. Then make both versions also write each block's total to a second output array, and say in one sentence which of the two hands you that total for free.

Time: 25 to 40 minutes. Submit: the new kernel, both totals arrays, and the sentence.

Check: the harness runs both exclusive kernels at every block size and compares each against the CPU reference exactly, printing the first index where either disagrees. It then checks both totals arrays against a per-tile sum, which is the half that catches the interesting mistake, and prints the two side by side whether they pass or fail.

Hint 1

Write down what is sitting in the tile when the loop ends, in each version. Not what you copy out of it, what is in it. The two are scans of different arrays.

Hint 2

A tile's total is the last element of an inclusive scan of that tile. Which of your two versions ever performs an inclusive scan of the real input, and which one scans the input with its last element pushed off the end?

Solution

Shift-then-scan never puts the last element into a prefix. Thread tid loaded element tid - 1, so the tile holds an inclusive scan of elements 0 to n - 2 and the last input element was added to nothing. The final word is the total minus that element, so recovering the total costs another read and another add.

Scan-then-shift scans the real tile first, so its final word is already the total. The shift changes what each thread copies out, not what the block computed. Both routes are correct exclusive scans; the difference is what else you hold when the kernel ends, which is what day 32 comes to collect.

The framing that outlives this exercise: when one primitive is built out of another, pick the order that leaves the by-product you will need, because recomputing it later costs a pass over memory you have already read.

Pitfalls

Your scan is right at 32 threads and wrong at 256. That is the in-place race, and the block size is the only thing that changed.

Write to the half you are not reading. Compute Sanitizer's racecheck reports it whether or not the answer came out wrong, needs no root and works on a hosted tier: compute-sanitizer --tool racecheck ./scan.

Day 62 reads that output on a scan.

You double buffer and the low elements go stale. The threads with no element to their left stopped writing, because in the in-place version that was correct. Every thread writes every step; the ones with nothing to add copy themselves across.

Your tile is full of numbers nobody wrote. The load filled one half and the missing copy reads the other, so the block scans whatever the SM happened to be holding. compute-sanitizer --tool initcheck names it.

Without the tool the symptom is an answer that changes when an unrelated kernel runs first.

You add a second __syncthreads() to be safe. With two halves it buys nothing, and every thread in the block waits for the slowest. Work out which read it separates from which write before you add one.

You keep the scan in place and reach for volatile. Same instinct as the volatile reduction tail day 24 takes apart, and the same failure.

volatile stops the compiler caching a value in a register. The compiler was never the problem; the other threads are.

You time a single-block scan and conclude scan is slow. A launch that scans 256 elements is mostly launch. Scan is worth timing when a grid covers millions of elements, which is day 32.

Go deeper

Next

Day 32 replaces this scan with Blelloch's, which does O(n) additions by walking a tree up and then back down, and spends the block totals you just met on scanning a million elements across three launches. Day 33 uses the exclusive scan as an output index, and day 35 uses it for a radix sort.