code/day77-clusters/cluster_reduction.cuThis is the source used by the lesson and its recorded evidence. Compile commands and expected output live in the directory README.
// SPDX-License-Identifier: MIT
//
// Day 77: thread block clusters and distributed shared memory.
//
// One reduction, five schedules. Every mode reduces the same 64 MiB of
// ones with day 24's two-loads-plus-sequential-tree block kernel, then
// combines the per-block partials into one float five different ways:
//
// atomic every block atomicAdds its partial to one global address.
// 32,767 blocks, 32,767 atomics on the same word.
// dsm-1 the cluster path with cluster size 1. Every block is rank 0
// of its own cluster, so the atomic count does not move; this
// row prices the cluster launch machinery itself.
// dsm-2/4/8 blocks in a cluster write their partial into their own
// shared memory, cluster.sync(), and rank 0 reads the others'
// partials through cluster.map_shared_rank() before issuing
// one atomic per cluster. At 8, the global atomic count drops
// from 32,767 to 4,096.
//
// The input is all 1.0f and the count is under 2^24, so every partial and
// every running total is a whole number float represents exactly, whatever
// order the atomics land in. That is day 24's trick and it is why the
// check below is an exact comparison with no tolerance.
//
// What this program does not check: that any tail is faster than any
// other. Both tails hang off the same 64 MiB of DRAM reads, and whether
// 28,671 fewer single-address atomics are visible over that is exactly
// the claim the run exists to test.
//
// Needs compute capability 9.0: clusters and distributed shared memory
// exist on Hopper and every Blackwell, including RTX 50 cards (CC 12.0),
// and on nothing older. On the course's Tesla T4 this program compiles
// with -arch=sm_90 and refuses to run, with a real branch, before it
// allocates anything.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_90 -o cluster_reduction \
// cluster_reduction.cu
// Run: ./cluster_reduction
//
// -arch=sm_90, not sm_90a: nothing here is architecture-specific, and
// sm_90 embeds compute_90 PTX that JIT-compiles forward onto CC 10.x,
// 11.0 and 12.x cards. An sm_90a binary loads on nothing but Hopper.
//
// UNVERIFIED: not yet compiled or run on real hardware. Do not publish any
// output as this program's output until it has run. See
// research/REVIEW-PROCESS.md.
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <cooperative_groups.h>
#include <cuda_runtime.h>
// The one error macro. This file is standalone, the way a Compiler Explorer
// embed is, so it carries its own verbatim copy. `err_` carries a trailing
// underscore so it cannot collide with a variable at the call site, and the
// do/while makes the macro one statement so it survives a braceless `if`.
#define CUDA_CHECK(call) \
do { \
cudaError_t err_ = (call); \
if (err_ != cudaSuccess) { \
std::fprintf(stderr, "CUDA error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, cudaGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
namespace cg = cooperative_groups;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// 8 is the largest portable cluster size: "a maximum of 8 thread blocks in
// a cluster is supported as a portable cluster size in CUDA" (programming
// guide 2.2.1, CUDA 12.6.2 archive, checked 2026-09-01). Larger sizes are
// an architecture-specific opt-in this program does not take; main() prints
// what cudaOccupancyMaxPotentialClusterSize says this card would allow.
constexpr int kMaxClusterBlocks = 8;
// Day 24's element count. 2^24 minus 611: the sum is 16,776,605, under
// 2^24, so every partial sum lands on a float exactly, and 611 is 13 x 47,
// so no block size divides the count and every bounds check runs.
constexpr size_t kElems = 16777216ull - 611ull;
// Each block consumes two elements per thread, day 24's version 4.
constexpr int kElemsPerBlock = 2 * kThreadsPerBlock;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert((kThreadsPerBlock & (kThreadsPerBlock - 1)) == 0,
"the tree halves the stride, so the block size must be a "
"power of two");
// Day 24's two-loads block reduction, as a device function both kernels
// share so the tails are the only thing that differs between them.
//
// One thread: loads elements i and i + blockDim.x, adds them, then runs
// the sequential-addressing tree over the block's shared tile.
//
// One warp: each global load is 32 consecutive floats, 128 contiguous
// bytes, four 32-byte sectors. The tree's active lanes read consecutive
// shared words, one bank each, so nothing serializes.
//
// Launch assumption: exactly kThreadsPerBlock threads per block. Every
// thread reaches every __syncthreads(); the guards cover loads, never
// barriers. After the call, tile[0] holds the block's partial.
__device__ void reduceBlockToTile(const float* __restrict__ in, size_t n,
float* tile) {
const unsigned int tid = threadIdx.x;
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) * 2 + tid;
float sum = (i < n) ? in[i] : 0.0f;
if (i + blockDim.x < n) {
sum += in[i + blockDim.x];
}
tile[tid] = sum;
__syncthreads();
for (unsigned int s = kThreadsPerBlock / 2; s > 0; s >>= 1) {
if (tid < s) {
tile[tid] += tile[tid + s];
}
__syncthreads();
}
}
// The baseline tail. One global atomic per block, all 32,767 of them
// landing on the same word of global memory.
//
// Launch assumption: plain <<<>>>, no cluster.
__global__ void reduceBlockAtomic(const float* __restrict__ in,
float* __restrict__ sum, size_t n) {
__shared__ float tile[kThreadsPerBlock];
reduceBlockToTile(in, n, tile);
// snippet: atomic-tail
if (threadIdx.x == 0) {
atomicAdd(sum, tile[0]);
}
// end snippet
}
// The cluster tail. Blocks in a cluster combine their partials through
// distributed shared memory, and only rank 0 goes to global memory.
//
// Launch assumption: launched through cudaLaunchKernelEx with the
// cudaLaunchAttributeClusterDimension attribute set and a grid that is a
// multiple of the cluster size. cg::this_cluster() is only meaningful
// under that launch (a non-cluster launch on CC 9.0+ reads as a 1x1x1
// cluster). The two cluster.sync() calls are both load bearing: the first
// makes every block's tile[0] visible before any remote read, the second
// keeps every block alive until rank 0 has finished reading, because a
// block's shared memory ceases to exist when the block exits.
__global__ void reduceClusterDsm(const float* __restrict__ in,
float* __restrict__ sum, size_t n) {
__shared__ float tile[kThreadsPerBlock];
reduceBlockToTile(in, n, tile);
// snippet: cluster-combine
cg::cluster_group cluster = cg::this_cluster();
// Every block's partial is in its own tile[0]. The barrier makes all
// of them visible to every block in the cluster, and guarantees no
// block has raced ahead into the combine.
cluster.sync();
if (cluster.block_rank() == 0 && threadIdx.x == 0) {
float total = 0.0f;
for (int r = 0; r < static_cast<int>(cluster.num_blocks()); ++r) {
// The same shared array, seen through block r's mapping. Rank
// 0's own tile comes back unchanged when r is its own rank.
const float* remote = cluster.map_shared_rank(tile, r);
total += remote[0];
}
atomicAdd(sum, total);
}
// Without this barrier a non-zero rank may exit while rank 0 is still
// reading its shared memory, and the read lands in a dead address
// space. The failure is a race, not an error message.
cluster.sync();
// end snippet
}
// Launches the cluster kernel with a runtime cluster size, day 58's
// cudaLaunchKernelEx machinery with the cluster attribute in place of the
// PDL one. The grid is padded up to a multiple of the cluster size, which
// the launch requires; padded blocks read nothing (their whole range is
// past n) and contribute an exact 0.0f.
static void launchClusterDsm(int clusterSize, int blocks, const float* d_in,
float* d_sum, size_t n) {
const int padded = ((blocks + clusterSize - 1) / clusterSize) * clusterSize;
// snippet: cluster-launch
cudaLaunchAttribute attrs[1];
attrs[0].id = cudaLaunchAttributeClusterDimension;
attrs[0].val.clusterDim.x = static_cast<unsigned int>(clusterSize);
attrs[0].val.clusterDim.y = 1;
attrs[0].val.clusterDim.z = 1;
cudaLaunchConfig_t cfg = {};
cfg.gridDim = dim3(static_cast<unsigned int>(padded));
cfg.blockDim = dim3(kThreadsPerBlock);
cfg.dynamicSmemBytes = 0;
cfg.stream = nullptr; // the legacy default stream
cfg.attrs = attrs;
cfg.numAttrs = 1;
CUDA_CHECK(cudaLaunchKernelEx(&cfg, reduceClusterDsm, d_in, d_sum, n));
// end snippet
}
// Times a launch with CUDA events and returns the mean milliseconds per
// run. A host clock around a launch measures the launch, because launches
// are asynchronous. Day 9 takes that apart. Copy this helper verbatim.
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
// Warm up this kernel, not just the first kernel in the program. Lazy
// module loading has been the default since CUDA 12.2 on Linux, so the
// first launch of each kernel pays its own load.
for (int i = 0; i < kWarmupRuns; ++i) {
launch();
}
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start));
for (int i = 0; i < kTimedRuns; ++i) {
launch();
}
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms / kTimedRuns;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
std::printf("GPU: %s (compute capability %d.%d), %d SMs\n", prop.name,
prop.major, prop.minor, prop.multiProcessorCount);
// The capability gate. A real branch, before any allocation: clusters
// need CC 9.0, and on an older card the cluster launch would fail at
// run time anyway. Saying so plainly beats a driver error.
if (prop.major < 9) {
std::fprintf(stderr,
"this GPU is compute capability %d.%d; thread block "
"clusters need 9.0 or newer (Hopper, or any "
"Blackwell including RTX 50)\n",
prop.major, prop.minor);
return EXIT_FAILURE;
}
int clusterLaunch = 0;
CUDA_CHECK(cudaDeviceGetAttribute(&clusterLaunch, cudaDevAttrClusterLaunch,
device));
if (clusterLaunch == 0) {
std::fprintf(stderr,
"cudaDevAttrClusterLaunch is 0 on this device; "
"cluster launch is not supported here\n");
return EXIT_FAILURE;
}
// What the biggest cluster this kernel could launch at is, on this
// card, according to the occupancy API. Eight is the portable upper
// bound; smaller GPUs and MIG configurations may report less.
int maxClusterSize = 0;
{
cudaLaunchConfig_t cfg = {};
cfg.gridDim = dim3(1);
cfg.blockDim = dim3(kThreadsPerBlock);
cfg.dynamicSmemBytes = 0;
CUDA_CHECK(cudaOccupancyMaxPotentialClusterSize(
&maxClusterSize, reduceClusterDsm, &cfg));
}
std::printf("max potential cluster size for this kernel: %d\n",
maxClusterSize);
// Every failure below records itself and falls through to the one
// cleanup block at the bottom, so no path can return with device
// memory still allocated.
int failures = 0;
const size_t inputBytes = kElems * sizeof(float);
const int blocks =
static_cast<int>((kElems + kElemsPerBlock - 1) / kElemsPerBlock);
const float want = static_cast<float>(kElems);
float* d_in = nullptr;
float* d_sum = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, inputBytes));
CUDA_CHECK(cudaMalloc(&d_sum, sizeof(float)));
std::vector<float> h_in(kElems, 1.0f);
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), inputBytes, cudaMemcpyHostToDevice));
std::printf("\n%zu elements, %.0f MiB, every element 1.0f, %d blocks\n",
kElems, static_cast<double>(inputBytes) / (1024.0 * 1024.0),
blocks);
// Correctness before any timing: one untimed run per mode against the
// exact expected sum. The accumulator is zeroed before each, and a
// wrong total names the mode and the difference in elements.
struct Mode {
const char* name;
int clusterSize; // 0 means the plain atomic baseline
};
const Mode modes[] = {{"atomic baseline", 0},
{"dsm, cluster of 1", 1},
{"dsm, cluster of 2", 2},
{"dsm, cluster of 4", 4},
{"dsm, cluster of 8", kMaxClusterBlocks}};
const int kModes = static_cast<int>(sizeof(modes) / sizeof(modes[0]));
for (int m = 0; m < kModes; ++m) {
if (modes[m].clusterSize > maxClusterSize) {
std::printf("%s: SKIP (device maximum is %d)\n", modes[m].name,
maxClusterSize);
continue;
}
CUDA_CHECK(cudaMemset(d_sum, 0, sizeof(float)));
if (modes[m].clusterSize == 0) {
reduceBlockAtomic<<<blocks, kThreadsPerBlock>>>(d_in, d_sum,
kElems);
CUDA_CHECK(cudaGetLastError());
} else {
launchClusterDsm(modes[m].clusterSize, blocks, d_in, d_sum, kElems);
}
CUDA_CHECK(cudaDeviceSynchronize());
float got = 0.0f;
CUDA_CHECK(
cudaMemcpy(&got, d_sum, sizeof(float), cudaMemcpyDeviceToHost));
if (got != want) {
std::fprintf(stderr,
"%s summed to %.9g, want %.9g; that is a "
"difference of %.9g elements\n",
modes[m].name, got, want, want - got);
++failures;
}
}
// The timing table. The timed loops accumulate into d_sum without
// resetting it, because a reset inside the loop would sit inside the
// measurement; the total is garbage during timing and nothing reads
// it. Correctness was settled above on the untimed runs.
std::printf("\n%-22s %14s %12s %10s %9s\n", "mode", "global atomics",
"time (ms)", "GB/s", "vs atomic");
std::printf("%-22s %14s %12s %10s %9s\n", "----------------------",
"--------------", "----------", "--------", "---------");
float atomicMs = 0.0f;
for (int m = 0; m < kModes; ++m) {
const int c = modes[m].clusterSize;
if (c > maxClusterSize) {
std::printf("%-22s %14s %12s %10s %9s\n", modes[m].name, "skipped",
"-", "-", "-");
continue;
}
float ms = 0.0f;
long atomics = 0;
if (c == 0) {
ms = timeKernel([&] {
reduceBlockAtomic<<<blocks, kThreadsPerBlock>>>(d_in, d_sum,
kElems);
});
atomics = blocks;
atomicMs = ms;
} else {
ms = timeKernel(
[&] { launchClusterDsm(c, blocks, d_in, d_sum, kElems); });
const int padded = ((blocks + c - 1) / c) * c;
atomics = padded / c;
}
const double gbps = static_cast<double>(inputBytes) /
(static_cast<double>(ms) * 1.0e-3) / 1.0e9;
std::printf("%-22s %14ld %12.3f %10.1f %9.2f\n", modes[m].name, atomics,
ms, gbps, atomicMs / ms);
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_sum));
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}