CUDA stream compaction with a prefix sum

Someone on the NVIDIA forums, describing how he packs particles into bins in a molecular dynamics code:

"In my application (MD) fast atomics on Fermi (even faster on Kepler) I find that an atomic op version of this is significantly faster. I just run one thread per element, determine the bin, atomicAdd to that bin's size counter, and then write the data to the appropriate location in the bin's data (location returned by the atomicAdd). Due to indeterministic thread execution, items end up in a different order each time the kernel is run."

DrAnderson42, 8 February 2013, https://forums.developer.nvidia.com/t/opposite-of-stream-compaction/28201 (checked 2026-08-30)

He is right about both halves and he says so plainly. Most people who write that kernel only notice the first half. The counter is fast, the answer is a different arrangement every run, and nothing in the code says so.

This page is about the version that does not have that property, and about why it is the most reusable thing in this module.

The counter you cannot have, and the sum that replaces it

Start from the loop everybody already knows:

size_t m = 0;
for (size_t i = 0; i < n; ++i) {
    if (in[i] != 0) {
        out[m++] = in[i];
    }
}

m++ is the whole problem. Every iteration reads what the previous one wrote, which is exactly the shape you cannot spread across eight million threads.

So look at what m holds. When the loop reaches element i, m is the number of survivors strictly before i. That is not a counter, it is a function of i, and you can compute it for every i at once.

Write down a 1 for every element that survives and a 0 for every element that does not, take the exclusive prefix sum of that array, and entry i is the slot element i belongs in.

Three steps, and each of them is a kernel you have already met:

1. Predicate. One comparison per element, no thread looking at any other.
2. Scan. The prefix sum from day 31, which does not know or care what it is scanning.
3. Scatter. Every surviving thread writes once, to the slot it was handed.

A forum regular states it in one sentence: "The general approach is to create a flag array indicating which elements to copy. Then create a prefix sum of the flag array to find the output positions of elements to copy. Last, perform the copy."

(striker159, 27 October 2022, https://forums.developer.nvidia.com/t/how-to-put-specific-elements-from-one-array-to-another-array-use-cuda/232010 , checked 2026-08-30.) An NVIDIA moderator on the same thread calls it "the canonical method".

The scan also provides the total with no extra pass. The top of the last level is the number of survivors, already sitting in global memory on the device, so nothing downstream has to wait for a copy to the host before it can run.

Why the fast version answers a different question

The intuition that sends you to atomicAdd is a good one. m is a counter that every iteration increments, the GPU has an instruction for exactly that, and day 26 showed you how to use it. The translation is two lines and it is correct: every survivor gets a distinct slot and the final counter is the right total.

What it silently drops is the order. Blocks reach the counter in whatever sequence the scheduler runs them, so element eight million can land ahead of element three.

Many common tests still pass: the count, sum, and set membership all match. The failure surfaces later, in whatever consumes the output and assumed the input's order was still there.

That is not an argument against atomics. DrAnderson42's bins do not care, and he says so. It is an argument for knowing which of the two you built, because they are different functions and only one preserves order.

The second reason to know the scan version is that a partition, a split by any predicate, and day 35's radix sort are all the same three steps with a different line in the middle.

Three phases, and a floor to price them against

Full program in code/day33-compaction/compaction.cu. The pattern is three of its kernels. The rest is the scan those three use and the discipline around the measurement.

One predicate kernel per problem, one scan, one scatter. markNonZero and markEven differ only in the line that computes the predicate, and everything downstream of them is shared. Those kernels and both scatters walk the array with the grid-stride loop from day 8, so the grid is a tuning knob rather than a correctness requirement.

Only the tile scan needs a fixed launch shape. Reusing the same scan and scatter matters more here than raw speed.

A copy is the baseline. copyInts reads every int and writes every int with no predicate, no scan and no scatter, from the same grid over the same buffers in the same process. Borrowing a copy bandwidth from day 11 would fold that day's block shape and buffer size into this page's answer.

Every row is checked before it is timed, and the two compactions are checked differently: the scanned one against the host element for element, the atomic one only for holding the same values. A test that demanded the atomic version's order would fail on scheduling, which is worse than no test.

Events, and a warm-up per kernel. Day 9 covers why a host clock around a launch measures the launch.

The scatter completes the pattern:

__global__ void scatterCompact(const int* __restrict__ in,
                               const int* __restrict__ rank,
                               int* __restrict__ out, size_t n) {
    const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
    for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
         i < n; i += step) {
        const int value = in[i];
        if (value != 0) {
            out[rank[i]] = value;
        }
    }
}

The word scatter oversells it. Destinations rise with i, so one warp's survivors land in a run of at most 32 consecutive slots, and a warp that keeps everything writes 128 contiguous bytes exactly like a copy. This is not the random write pattern that ruins coalescing.

The partition is the same kernel with one piece of arithmetic added:

__global__ void scatterPartition(const int* __restrict__ in,
                                 const int* __restrict__ rank,
                                 int* __restrict__ out, size_t n,
                                 const int* __restrict__ totalEven) {
    const size_t base = static_cast<size_t>(*totalEven);
    const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
    for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
         i < n; i += step) {
        const int value = in[i];
        const size_t evensBefore = static_cast<size_t>(rank[i]);
        if ((value & 1) == 0) {
            out[evensBefore] = value;
        } else {
            out[base + (i - evensBefore)] = value;
        }
    }
}

rank[i] is the number of even elements before i, so i - rank[i] is the number of odd ones. The second counter this pattern appears to need never exists, and totalEven arrives as a pointer into device memory rather than as an int, so the count never has to reach the host.

The comparison kernel is two lines shorter than either:

__global__ void compactAtomic(const int* __restrict__ in, int* __restrict__ out,
                              unsigned int* __restrict__ count, size_t n) {
    const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
    for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
         i < n; i += step) {
        const int value = in[i];
        if (value != 0) {
            const unsigned int slot = atomicAdd(count, 1u);
            out[slot] = value;
        }
    }
}

Note. The tile scan is the one kernel in the file whose pointers are not __restrict__, and the first draft got that wrong by copying the habit from every other kernel in the file. The scan runs in place, so the caller passes one array as both the input and the output. __restrict__ promises the memory is reached through that pointer only, which is a promise the caller has already broken, and the compiler is entitled to believe it. In place is still safe, for a different reason: a thread reads in[i] into a register before the first barrier and writes out[i] after the last one, and no other thread touches index i.

The honest caveat is that three phases means three passes. The program writes an array of ranks to global memory and reads it back, where a production implementation keeps the per-tile ranks in shared memory and passes a carry between tiles. That costs traffic on a memory-bound kernel, and it is paid on purpose here, because the three steps being separate is the thing this page is trying to show.

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)
40 SMs, 49152 bytes of shared memory per block
Input: 8387997 ints (32.00 MiB per buffer)
Grid: 1024 blocks x 256 threads for the stride kernels
Scan: 256 elements per tile, 32766 then 128 then 1 tiles
Mean of 10 runs after 3 warm-ups. Host copies are outside every timing.

Where the time goes, on the 50.0 percent input
row                                   launches   time (ms)   x copy
-----------------------------         --------  ----------  -------
copy in to out, no compaction                1       0.300     1.00
mark the predicate                           1       0.298     0.99
mark, then scan                              6       1.400     4.67
scatter alone, ranks ready                   1       0.443     1.48
compaction, all three phases                 7       1.843     6.14
partition even then odd                      7       1.892     6.31
compaction, one atomic counter               2       0.498     1.66

The copy row moved 33551988 bytes each way at 223.7 GB/s.
The partition put 6275981 even elements first.

Survival rate sweep
  kept    survivors   scan (ms)   atomic (ms) scan/atom   out of order
 -----  -----------  ----------   -----------  --------  -------------
 10.2%       852058       1.724         0.457      3.77         848797
 50.0%      4194345       1.843         0.498      3.70        4177744
 89.9%      7537218       1.772         0.521      3.40        7507495

Every row above was checked before it was timed. The scanned output
had to match the host element for element; the atomic output only had
to hold the same values, and the last column counts how many of its
slots hold something else.

Compaction costs about six times a straight copy, and the scan takes most of that time. Marking the predicate is free at 0.99x copy, inside the noise of the copy row itself. Add the scan and it becomes 4.67x.

The scatter on its own is 1.48x. All three phases together are 6.14x, across seven kernel launches.

That decomposition is the point of the table. If compaction is too slow, the scan is what to attack, and day 39 shows what a tuned library scan does to that row.

Partition costs almost the same as compaction, 6.31x against 6.14x, because it is the same machinery run once with the predicate and once with its negation, sharing the scan. Two outputs for three percent more time.

Run it yourself

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

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

No embed on this page, and the blocker is time, not memory. The tile scan uses 2 KiB of shared memory, but the program allocates about 128 MiB on each side, runs seven configurations for ten timed runs after three warm-ups each, and checks 8,387,997 elements against a host reference on every input.

That leaves no margin under Compiler Explorer's 20 second run cap. The three-level scan reaches 16,777,216 elements before it needs a fourth level.

Exercise

Turn the two-way partition into a three-way split, on an input that holds negatives, zeros and positives: every negative first, then every zero, then every positive, with the input's order preserved inside each group. Use one predicate and one scan per boundary rather than one per group, and say in one sentence why the third group needs no scan of its own.

Time: 30 to 45 minutes. Submit: your scatterThreeWay kernel and the sentence.

Check: the harness builds the three-way reference on the host, compares your whole output element for element with no tolerance, and prints the first index that disagrees with both values. It checks the two boundaries separately, so a split that is right inside each group but puts a boundary one slot early is reported as a boundary error rather than as a mismatch somewhere in the middle. It prints both boundary positions whether it passes or not.

Hint 1

You already know how to answer "how many elements before i satisfy P". The question is how many different predicates you have to ask that question about, and it is fewer than three.

Hint 2

The two-way version used one scan and got the second group from i - rank[i], because every element was in exactly one of two groups. With three groups, i - negRank[i] counts everything that is not negative. What second count do you need to split that remainder, and where does the base of the third group come from?

Solution

Two predicates, two scans, three destinations. Let a[i] be the number of negatives before i and b[i] the number of zeros before i, with totals A and B. Then a negative element goes to a[i], a zero goes to A + b[i], and a positive goes to A + B + (i - a[i] - b[i]), because everything before i that is neither negative nor zero is positive.

The third group needs no scan because its rank is arithmetic on the other two. That generalises: a k-way stable partition costs k-1 scans, never k, and the last group is always the leftover. A 4-bit radix pass is a 16-way partition, so day 35 pays for 15 counts and gets the sixteenth free.

The shipped program times the two-way partition at 1.892 ms, 6.31x the copy row. The three-way version adds one more scan and shares the rest, so measure it yourself and see how close to 6.31x it lands.

Pitfalls

Your compaction is correct and your output order is not. A global counter under atomic contention gives every survivor a distinct slot and the right total. The scheduler chooses the arrangement.

Check the order explicitly or use the scan. Day 26 covers what the counter costs.

You used an inclusive scan. The destination is the number of survivors strictly before i. An inclusive scan counts i as well, so every survivor lands one slot late and the last one writes past the end of the output.

Either scan exclusively or subtract the element's own flag.

You scattered in place. out[rank[i]] = in[i] with out and in the same array looks safe because rank[i] is never greater than i.

It is unsafe because thread i writes to slot rank[i], and the thread that owns that index may not have read it yet. Compaction needs a second buffer or a second kernel.

You marked an in-place kernel's pointers __restrict__. Passing one array as two restricted parameters tells the compiler two things that cannot both be true. It may hoist a load past a store and the answer changes with the optimisation level, which is the worst failure mode there is.

You copied the survivor count back to the host to size the next launch. The count is already on the device at the top of the scan. Reading it back turns seven asynchronous launches into a pipeline with a stall in the middle, and day 9 shows what that stall is worth.

You reported the atomic version as a speedup. It computes a different function. Publish the ratio next to the count of out-of-order slots or do not publish the ratio.

Go deeper

• CUDA Programming Guide, "Writing SIMT Kernels", for the atomic functions and their memory ordering: https://docs.nvidia.com/cuda/cuda-programming-guide/02-basics/writing-cuda-kernels.html (checked 2026-08-30)
• The forum thread where an NVIDIA moderator writes out the three-kernel structure and then says "In practice, writing a prefix sum is an operation I wouldn't want to write myself, although I probably could with enough effort": https://forums.developer.nvidia.com/t/how-to-put-specific-elements-from-one-array-to-another-array-use-cuda/232010 (checked 2026-08-30)
• Programming Massively Parallel Processors, 4th edition, chapter 11, "Prefix sum (scan): An introduction to work efficiency in parallel algorithms", and chapter 13 on sorting, where the same partition runs once per bit: https://shop.elsevier.com/books/programming-massively-parallel-processors/hwu/978-0-323-91231-0

Next

Day 34 takes the same array and asks a different question of every element, one that reads its neighbours, so the scan goes away and a halo comes back. Day 35 is this page run once per bit: a radix pass is a stable partition on one bit of the key, which is the two-way split you just wrote, and sorting is what you get for repeating it. Day 39 is where the course stops hand-writing these and measures what the libraries cost.