100 DAYS OF GPU Course preview
Day 25Module 3
in-technical-review

Optimizing a CUDA reduction to the memory ceiling

NVIDIA's reduction walkthrough, the deck this whole ladder comes from, ends its fifth version with a callout in capitals:

IMPORTANT: For this to be correct, we must use the "volatile" keyword!

Mark Harris, Optimizing Parallel Reduction in CUDA, slide 22, https://developer.download.nvidia.com/assets/cuda/files/reduction.pdf (checked 2026-08-30).

That was true in 2007. On every GPU this course targets, it is undefined behaviour. It usually still returns the right answer, which makes the race hard to spot.

Day 24 built v4, a shared-memory tree over two elements per thread. This page adds many elements per thread, finishes the last warp in registers instead of with volatile, moves the block size into a template argument, loads four floats per instruction, and adds a second pass so the total stays on the device.

It scores all seven rows against a copy kernel timed in the same process. That copy rate is the right limit for this reduction.

What a reduction actually spends its time on

Summing n floats reads n floats and writes one. It uses no second array and does one add per element, so the kernel is memory-bound at every size worth timing. The useful measure is the fraction of the card's memory bandwidth it reaches.

This program first times a plain copy over its own buffer in the same process, then scores every row against it. Day 11 measured 232.9 GB/s for a copy on the same node at a different size. That number is a useful check, but it is not the denominator here because both parts of the ratio must come from the same run.

The kernel has a load stage and a fold stage. Rebuilt here over 2^26 floats, day 24's v4 gives every thread two elements, so it launches 131,072 blocks. Each block then runs an eight-level fold through shared memory for 512 elements.

A grid-stride loop over a fixed grid of 1024 blocks lets each thread first add 256 elements into a register. The fold still has eight levels, but now covers 65,536 elements per block instead of 512. Its cost per element falls by a factor of 128.

Every version after it changes the fold. Version five stops using the block-wide __syncthreads() barrier once 32 values remain in one warp. Versions six and seven reduce the instruction count.

If the grid-stride loop makes the fold cost too small to measure, those later changes cannot save much time.

For the last 64 values, lane 0's shared-memory path uses 18 accesses and 6 barriers, while the shuffle tail uses 2 shared reads and 5 register exchanges.

The rule that stopped being true

The deck states the old reasoning plainly:

"Instructions are SIMD synchronous within a warp. That means when s <= 32: We don't need to __syncthreads()"

So the tail becomes six unguarded lines over a volatile pointer:

__device__ void warpReduce(volatile int* sdata, int tid) {
    sdata[tid] += sdata[tid + 32];
    sdata[tid] += sdata[tid + 16];
    sdata[tid] += sdata[tid + 8];
    sdata[tid] += sdata[tid + 4];
    sdata[tid] += sdata[tid + 2];
    sdata[tid] += sdata[tid + 1];
}

Volta removed the premise. With independent thread scheduling, the 32 lanes of a warp each have a program counter and need not run the same instruction at the same time. NVIDIA says: "The correctness of such a program depends on implicit warp-synchronous behavior, which may change from one hardware architecture to another", and "thread convergence is guaranteed only within explicitly synchronous warp-level primitives" (https://developer.nvidia.com/blog/using-cuda-warp-level-primitives/ , checked 2026-08-30).

volatile stops the compiler from caching the value. It does not make the lanes wait for one another.

The replacement is the warp shuffle from day 23. A lane reads another lane's register directly, and the mask names the lanes that take part. NVIDIA rewrote its own reduction sample this way, but the old deck remains online.

The deck's speedups do not apply to every GPU. It ran on a G80 with a quoted peak of 86.4 GB/s, while its v4 kernel reached 17.377 GB/s. The next slide says: "At 17 GB/s, we're far from bandwidth bound."

That kernel used about one fifth of the G80's quoted bandwidth, so cutting instructions helped. The table below tests the same changes on a T4.

Five rungs, and a ceiling measured in the same process

Full program in code/day25-reduction-2/reduction_2.cu.

Three rules, and the second is why this page quotes no number from another one.

Every reduction row reads the same bytes, and every GB/s comes from one constant. No row can look fast by reading less, which is the failure day 11 exists to warn about.

The ceiling is measured here, not quoted. A copy kernel over the same buffer runs first, in the same process on the same card, so the fraction column compares two numbers from one run.

Every row is checked before it is timed, and the check is exact. Every input element is 1.0f, so the answer is the element count and the sum counts how many accumulations happened. A kernel that drops one element in 2^26 fails; the same kernel on random data would be off by a relative 1.5e-8 and sail through any tolerance worth writing.

The first two rows are the baseline and the cascade. Day 24's kernel is rebuilt here rather than quoted from that page, because its problem size is different and a ratio across two pages is not a ratio. Then the same tree gets a grid-stride loop, which is the deck's own last rung arriving first.

Rung five replaces the volatile tail with this:

__device__ float warpReduceSum(float val) {
    for (unsigned int offset = kWarpSize / 2; offset > 0; offset >>= 1) {
        val += __shfl_down_sync(kFullMask, val, offset);
    }
    return val;
}

and stops the shared tree while 64 values are still live, so warp 0 folds them in pairs and then finishes in registers:

    // Stop at 64 live values. Below that the survivors are one warp, and a
    // warp does not need a block-wide barrier to talk to itself.
    for (unsigned int half = blockDim.x / 2; half > kWarpSize; half >>= 1) {
        if (tid < half) {
            tile[tid] += tile[tid + half];
        }
        __syncthreads();
    }

    if (tid < kWarpSize) {
        const float val = warpReduceSum(tile[tid] + tile[tid + kWarpSize]);
        if (tid == 0) {
            out[blockIdx.x] = val;
        }
    }

Version six replaces the three uses of blockDim.x with a template argument and sizes the tile from that argument. Only four lines change.

blockDim.x is a run-time value, so ptxas cannot know the loop's trip count. A template argument gives ptxas that value at compile time, so it emits the levels as straight-line code. This is loop unrolling.

// The second exception to the no-templates-before-module-8 rule in
// CUDA-CODE-STYLE.md, after timeKernel. It is not a style choice: v6 of the
// canonical ladder IS template unrolling, because the block size has to be a
// compile-time constant for the compiler to erase the loop and the bounds
// tests. A runtime blockDim.x cannot do it. The comparison against a library
// reduction belongs to day 39, where CCCL is introduced.
template <unsigned int kBlockSize>
__global__ void reduceUnrolled(const float* __restrict__ in,
                               float* __restrict__ out, size_t n) {
    __shared__ float tile[kBlockSize];

    const unsigned int tid = threadIdx.x;
    const size_t step = gridDim.x * static_cast<size_t>(kBlockSize);

Rung seven is the only rung above the cascade that touches the part of the kernel that dominates it. A vectorized load reads a float4, so one instruction fetches 16 bytes and a warp covers 512 contiguous bytes instead of 128:

    const size_t n4 = n / 4;
    const float4* in4 = reinterpret_cast<const float4*>(in);

    float acc = 0.0f;
    for (size_t i = t; i < n4; i += step) {
        const float4 v = in4[i];
        acc += v.x + v.y + v.z + v.w;
    }

    // The 0 to 3 elements that do not fill a float4, one each to the first
    // few threads of the grid.
    const size_t tailStart = n4 * 4;
    if (tailStart + t < n) {
        acc += in[tailStart + t];
    }

Rung eight is not a kernel. Every row above it leaves partial sums for somebody else to add up, and that somebody has been the host. A second launch over those partials makes the answer a single device float:

    reduceVectorized<kThreadsPerBlock>
        <<<kBlocks, kThreadsPerBlock>>>(d_in, d_partials, kElems);
    CUDA_CHECK(cudaGetLastError());
    reduceUnrolled<kThreadsPerBlock><<<1, kThreadsPerBlock>>>(
        d_partials, d_sum, static_cast<size_t>(kBlocks));
    CUDA_CHECK(cudaGetLastError());

Note. The rows are not all doing the same job, and the table says so in a column. Five of them stop at partial sums, so they are timing an unfinished reduction. v8 is the only row that leaves one number on the device.

Hardware. There is a one-instruction warp sum, __reduce_add_sync, and it is "Supported by devices of compute capability 8.x or higher" (https://docs.nvidia.com/cuda/cuda-programming-guide/05-appendices/cpp-language-extensions.html , checked 2026-08-30). The verification node is a Tesla T4 at 7.5, so this page uses __shfl_down_sync, which runs everywhere from Kepler onward and costs five instructions instead of one.

Results

Measured. Tesla T4, driver 595.84, CUDA 12.6 (V12.6.85), built with nvcc -O3 -arch=sm_75. Captured 2026-08-30 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

n = 67108864 floats, 256 MiB read per reduction
v4 runs 131072 blocks of 256 threads, two elements each. Every row below it
runs 1024 blocks of 256 threads, 256 elements each.
kernel                     time (ms)       GB/s   % of copy  produces
------------------------   ---------   --------   ---------  --------
copy (ceiling)                 2.669      201.2       100.0  n floats
v4 two loads/thread            1.321      203.3       101.0  partials
v4c cascaded loads             1.125      238.5       118.6  partials
v5 shuffle tail                1.132      237.1       117.9  partials
v6 template unrolled           1.136      236.4       117.5  partials
v7 float4 loads                0.956      280.8       139.6  partials
v8 two passes                  0.960      279.6       139.0  one float

The copy row's GB/s counts the bytes it reads and the bytes it
writes. Every reduction row counts only the 256 MiB it reads,
because that is all a reduction moves. A reduction above 100 is not
a broken measurement, it is a kernel that does not have to write.
The last two rows are the only ones that leave a single number on
the device. v4 leaves 131072 partials and the four rows above v8 leave
1024, for somebody else to finish. What v8 costs is that somebody.

One line of that transcript is wrong and stays as it was printed. The program's footer says "the last two rows" leave a single number on the device. Only v8 does; v7 leaves 1,024 partials.

The wording was left over from a CUB row cut before the run (DECISIONS.md, 2026-08-30), it is corrected in reduction_2.cu, and the transcript is left as the run produced it rather than edited. The evidence file says the same thing at the bottom.

The corrected footer appears in the CUDA 13.0 run. Version 4 moved from 1.321 to 2.889 ms and from 101.0% to 45.5% of the copy row. Versions 4c through 6 measured about 212 GB/s, while versions 7 and 8 stayed close to their CUDA 12.6 results at 278.4 and 277.0 GB/s.

All sums still passed. Only v8 produced one device result, and the order of the results from 4c onward did not change.

Read the percentage column carefully, because it goes above 100 and that is not a mistake. The copy row moves bytes in both directions and its GB/s counts both. A reduction only reads, so its row counts the 256 MiB it reads and nothing else.

A reduction at 139 percent of the copy row is not faster than the memory system. It moves half as many bytes per element.

Cascading the loads was worth more than every instruction-level trick after it. Version 4c reaches 238.5 GB/s against version 4's 203.3, simply by giving each thread 256 elements instead of 2. Versions 5 and 6, the shuffle tail and the template unroll, measure 237.1 and 236.4: within noise of 4c and of each other.

These results show what each change did on this GPU. The classic deck's large gains came from cutting instruction cost on a G80 that was far below its memory bandwidth. This T4 is memory-bound from version 4c onward, so removing more instructions does not change the measured time.

Version 7 moves the needle again because it changes the memory access, not the arithmetic. float4 loads reach 280.8 GB/s. Wider loads mean fewer, larger memory transactions, which is the same lever day 11 pulled.

Run it yourself

Colab's free T4 or any card you own. The build line is in the repo's README:

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

There is no Compiler Explorer embed. The program allocates two 256 MiB device buffers and a 256 MiB host vector, then runs eighty timed launches over seven rows. That work exceeds Compiler Explorer's 20 second run limit.

The program writes out each version without a new library. Day 39 adds the library comparison and introduces CCCL.

Exercise

The grid is fixed at 1024 blocks. Make that a runtime knob, sweep it over 64, 256, 1024, 4096 and 16384 blocks with v8, and report each one's fraction of the copy row. Then say in one sentence why the curve has a flat middle rather than a peak.

Keep the count a power of two. The file's exactness argument needs every thread to take the same number of elements, and a static_assert will stop the build the moment that stops being true.

Time: 30 to 45 minutes. Submit: the five fractions, the block count you would ship, and the sentence.

Check: the harness re-runs the exact sum check at every grid size, so a configuration that drops elements fails. It prints the result and the expected total.

The harness reports each grid size's fraction of the copy row whether the sum passes or not, because the curve is the result to study. It checks the sum, not the time. A time limit would decide the question that the sweep asks you to test.

Hint 1

Adding blocks changes two costs in opposite ways. More blocks let more of the card work, but they also run the eight-level fold more times for the same input. Which cost matters at 64 blocks, and which matters at 16384?

Hint 2

At 64 blocks the grid is 16,384 threads on a T4 that can keep 40,960 resident. At 16384 blocks each thread reads 16 elements and then folds eight levels. Work out the folds per element at both ends and you have the curve.

Solution

The shape, not a number, because absolute values move by an order of magnitude between a T4 and an H100.

Too few blocks and the card is idle. At 64 blocks the grid is 16,384 threads against 40,960 resident slots, so most of the machine has nothing to do and the fraction of the copy row is poor no matter how good the kernel is.

Too many blocks and the fold stops being free. Every block pays eight tree levels and a barrier per level for whatever it loaded, so at 16384 blocks that cost is spread over 16 elements per thread instead of 256, and the thing rungs five to seven were shaving is suddenly back in the measurement.

Between those values, a range of block counts gives the same result. Pick a value in that range. There is no reason to tune it further.

Pitfalls

You copy the volatile tail out of a tutorial and it works. It may return the right answer in many runs. Since Volta, warp lanes need not run the same instruction at the same time, so the tail has a race.

Use __shfl_down_sync, or __syncwarp() if you must stay in shared memory. Do not rely on lockstep execution.

You pass a mask that names a lane which is not there. __shfl_down_sync takes a mask of participating lanes, and the tail on this page passes 0xffffffff only because it runs under if (tid < 32), where all 32 lanes of warp 0 are active. Put a shuffle inside a divergent branch with a full mask and the result is undefined, not merely wrong.

Your float4 cast is not 16-byte aligned. cudaMalloc returns a pointer aligned to at least 256 bytes, so casting the base of an allocation is safe. Casting d_in + 1, or a pointer into the middle of a struct, is not, and it fails as cudaErrorMisalignedAddress rather than as a wrong answer.

Your vectorised kernel silently drops up to three elements. Dividing the element count by four does not handle the remainder. If the count is not a multiple of four, process the remaining elements with scalar loads.

A result that misses three elements out of millions may still pass a loose tolerance. Check the exact sum with a known input.

You compare a reduction that stops at partial sums with one that does not. Five of the rows on this page leave a buffer of partials behind. A row that finishes the job has done strictly more work, so the comparison is only fair once you say which is which, which is what the produces column is for.

You profile as a normal user and get nothing. ncu reports ERR_NVGPUCTRPERM, in full: The user running <tool_name/application_name> does not have permission to access NVIDIA GPU Performance Counters or the Hardware Event System on the target device (https://developer.nvidia.com/nvidia-development-tools-solutions-err_nvgpuctrperm-permission-issue-performance-counters , checked 2026-08-30).

sudo ncu gets past it where you have root, because RmProfilingAdminOnly defaults to 1 on a stock driver. Neither Colab nor Kaggle gives you root, so on those tiers the timing table is the whole measurement and the Nsight Compute commands in the repo's README are for a machine you own.

Go deeper

Next

Day 26 sums the same array with a single global atomicAdd and finds where that crosses the tree you just built. Day 27 explains why the two-pass version above is the safe way to finish a multi-block reduction and a grid-wide flag is not, and day 28 removes the second launch with a cooperative grid sync. Day 39 is where CUB arrives and DeviceReduce::Sum goes head to head with the ladder you just built, on scan and sort as well.