code/day92-nccl/allreduce_grad.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 92: NCCL and collective operations.
//
// One gradient buffer per visible GPU, summed across all of them with
// ncclAllReduce, then checked against a double-precision host sum of the
// buffers the devices actually produced. Every rank must finish holding
// the same numbers. That is what makes it an all-reduce and not a
// reduce, and the program checks it as a separate claim.
//
// Single process, one thread, one communicator per device built in one
// ncclCommInitAll call. Every collective here sits inside
// ncclGroupStart / ncclGroupEnd, because a single thread driving several
// devices has no choice: "every NCCL call may have to block, waiting for
// other threads/ranks to arrive, before effectively posting the NCCL
// operation on the given stream" (NCCL user guide, Group Calls, checked
// 2026-09-01). Drop the group and rank 0's call waits for a rank this
// same thread has not reached yet, which is a hang, not an error.
//
// Two or more GPUs are needed to measure communication. With one visible
// device the program builds a one-rank clique and runs anyway: the
// all-reduce then reduces one contribution, so the output equals the
// input, which grades the API wiring and the buffer contract and says
// nothing about a link. The banner it prints says exactly that.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o allreduce_grad \
// allreduce_grad.cu -lnccl
// Run: ./allreduce_grad
//
// PARTIALLY VERIFIED 2026-09-02: the one-rank fallback passed on a Tesla T4
// with NCCL 2.30.7.1 and correctly reported both bandwidth fields as n/a.
// The required two-GPU Kaggle run remains pending. See evidence/.
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <vector>
#include <cuda_runtime.h>
#include <nccl.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)
// Same shape for the library. A collective that returns ncclInvalidUsage
// and is not read reappears later as a hang on a different line, so every
// NCCL call in this file is wrapped exactly the way every CUDA call is.
// snippet: nccl-check
#define NCCL_CHECK(call) \
do { \
ncclResult_t err_ = (call); \
if (err_ != ncclSuccess) { \
std::fprintf(stderr, "NCCL error %s:%d: %s: %s\n", __FILE__, \
__LINE__, #call, ncclGetErrorString(err_)); \
std::exit(EXIT_FAILURE); \
} \
} while (0)
// end snippet
// 4 Mi floats plus 611 is 16 MiB per rank, small enough for a 16 GB T4 to
// hold four buffers and large enough that a PCIe transfer is not pure
// latency. The 611 is there so the kernel's bounds guard runs on every
// launch instead of never.
constexpr size_t kElems = 4ull * 1024ull * 1024ull + 611ull;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kMaxRanks = 8;
constexpr float kRelTolerance = 1e-5f;
constexpr int kMinCapability = 75;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
// grad[i] = (rank + 1) / ((i % 97) + 3). One thread owns one element.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes, one coalesced store, and nothing
// is read from global memory at all.
//
// Launch assumption: gridDim.x * blockDim.x >= n. The values are
// reciprocals of small integers, so hardly any of them is exact in binary
// floating point. That is deliberate: a buffer of whole numbers would
// survive any summation order and would not exercise the tolerance a real
// gradient reduction needs.
// snippet: fill-kernel
__global__ void fillLocalGradient(float* __restrict__ grad, size_t n,
int rank) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
grad[i] =
static_cast<float>(rank + 1) / static_cast<float>((i % 97) + 3);
}
}
// end snippet
// Sums the per-rank buffers element by element in double.
//
// It sums the copies the devices produced, not the kernel's formula. A
// reference that recomputed `(rank + 1) / ((i % 97) + 3)` would be
// checking the fill against itself; this one checks the collective, which
// is the only thing on trial here. It never allocates: the caller owns
// every buffer.
static void sumRanksCpu(const std::vector<std::vector<float>>& locals,
double* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
double total = 0.0;
for (size_t r = 0; r < locals.size(); ++r) {
total += static_cast<double>(locals[r][i]);
}
out[i] = total;
}
}
// The relative tolerance, derived rather than picked.
//
// Every output element is a sum of `ranks` floats, so the reduction depth
// is the rank count. EXERCISE-DESIGN.md gives rtol_K = max(1e-5, 4 * eps *
// sqrt(K)) with eps = 2^-23; at K = 8 the second term is about 1.3e-6, so
// the table value still wins and the sqrt term only starts to bite in the
// hundreds of ranks. Writing it out means the number on the page comes
// from the depth instead of from habit.
static float toleranceFor(int ranks) {
const float eps = 1.1920929e-7f; // 2^-23, the f32 spacing at 1.0
const float bound = 4.0f * eps * std::sqrt(static_cast<float>(ranks));
return (bound > kRelTolerance) ? bound : kRelTolerance;
}
// Returns the first index where got and want differ by more than the
// relative tolerance, or n if they agree everywhere. Returning the index
// rather than a bool is the point: "wrong at 512" names the chunk of the
// ring that lost data, "wrong" does not.
static size_t firstMismatch(const float* got, const double* want, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const double scale = (want[i] == 0.0) ? 1.0 : std::fabs(want[i]);
if (std::fabs(static_cast<double>(got[i]) - want[i]) >
static_cast<double>(relTolerance) * scale) {
return i;
}
}
return n;
}
int main() {
int visible = 0;
CUDA_CHECK(cudaGetDeviceCount(&visible));
if (visible < 1) {
std::fprintf(stderr, "no CUDA device visible\n");
return EXIT_FAILURE;
}
const int ranks = (visible > kMaxRanks) ? kMaxRanks : visible;
int ncclVersion = 0;
NCCL_CHECK(ncclGetVersion(&ncclVersion));
std::printf("NCCL headers %d.%d.%d, library version code %d\n", NCCL_MAJOR,
NCCL_MINOR, NCCL_PATCH, ncclVersion);
std::vector<int> devs(static_cast<size_t>(ranks));
for (int r = 0; r < ranks; ++r) {
devs[static_cast<size_t>(r)] = r;
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, r));
std::printf("rank %d: %s (compute capability %d.%d)\n", r, prop.name,
prop.major, prop.minor);
if (prop.major * 10 + prop.minor < kMinCapability) {
std::fprintf(stderr,
"device %d is compute capability %d.%d; this course "
"targets 7.5 and newer\n",
r, prop.major, prop.minor);
return EXIT_FAILURE;
}
}
if (ranks == 1) {
std::printf(
"\nONE VISIBLE GPU. Building a one-rank communicator, which is\n"
"the fallback path: the all-reduce sums a single contribution,\n"
"so the output equals the input. It grades the API wiring and\n"
"the buffer contract. No byte crosses a link, the ring has no\n"
"neighbour, and the bandwidth columns below are not a\n"
"measurement of anything. Kaggle's T4 x2 is the free two-GPU\n"
"tier this program is written for.\n");
}
// Everything below records failures and falls through to one cleanup
// block, so no path returns with device memory, streams, events or
// communicators still alive.
int failures = 0;
const size_t bytes = kElems * sizeof(float);
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
const size_t nr = static_cast<size_t>(ranks);
std::vector<ncclComm_t> comms(nr, nullptr);
std::vector<cudaStream_t> streams(nr, nullptr);
std::vector<cudaEvent_t> starts(nr, nullptr);
std::vector<cudaEvent_t> stops(nr, nullptr);
std::vector<float*> d_grad(nr, nullptr);
std::vector<float*> d_out(nr, nullptr);
std::vector<std::vector<float>> h_local(nr, std::vector<float>(kElems));
std::vector<std::vector<float>> h_out(nr, std::vector<float>(kElems));
std::vector<double> h_want(kElems);
// One call builds the whole clique, which is what makes this the
// single-process shape. The multi-process shape is ncclGetUniqueId on
// one rank, broadcast by MPI or a socket, then ncclCommInitRank on
// each; the collective calls below are identical either way.
// snippet: init
NCCL_CHECK(ncclCommInitAll(comms.data(), ranks, devs.data()));
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaStreamCreate(&streams[k]));
CUDA_CHECK(cudaMalloc(&d_grad[k], bytes));
CUDA_CHECK(cudaMalloc(&d_out[k], bytes));
fillLocalGradient<<<blocks, kThreadsPerBlock, 0, streams[k]>>>(
d_grad[k], kElems, r);
CUDA_CHECK(cudaGetLastError());
}
// end snippet
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaEventCreate(&starts[k]));
CUDA_CHECK(cudaEventCreate(&stops[k]));
CUDA_CHECK(cudaMemcpyAsync(h_local[k].data(), d_grad[k], bytes,
cudaMemcpyDeviceToHost, streams[k]));
CUDA_CHECK(cudaStreamSynchronize(streams[k]));
}
sumRanksCpu(h_local, h_want.data(), kElems);
// The collective, issued once and checked before anything is timed.
// Inside a group these calls may return without having enqueued
// anything: the user guide says stream operations "can return without
// having enqueued the operation on the stream" and that
// cudaStreamSynchronize "can therefore be called only after
// ncclGroupEnd returns". So the synchronise loop sits below the group,
// never inside it.
// snippet: allreduce
NCCL_CHECK(ncclGroupStart());
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
NCCL_CHECK(ncclAllReduce(d_grad[k], d_out[k], kElems, ncclFloat,
ncclSum, comms[k], streams[k]));
}
NCCL_CHECK(ncclGroupEnd());
for (int r = 0; r < ranks; ++r) {
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaStreamSynchronize(streams[r]));
}
// end snippet
const float rtol = toleranceFor(ranks);
std::printf(
"\n%zu elements, %zu MiB per rank, %d rank(s), relative "
"tolerance %.3g\n",
kElems, bytes / (1024 * 1024), ranks, rtol);
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaMemcpy(h_out[k].data(), d_out[k], bytes,
cudaMemcpyDeviceToHost));
const size_t bad =
firstMismatch(h_out[k].data(), h_want.data(), kElems, rtol);
if (bad != kElems) {
std::fprintf(stderr, "rank %d wrong at %zu: got %.9g, want %.9g\n",
r, bad, static_cast<double>(h_out[k][bad]),
h_want[bad]);
++failures;
}
}
// The second claim, and the one the word all is doing: every rank ends
// up with the same bytes. A reduce would leave rank 0 right and the
// others holding their own input, which passes the check above on rank
// 0 alone.
for (int r = 1; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
if (std::memcmp(h_out[k].data(), h_out[0].data(), bytes) != 0) {
std::fprintf(stderr,
"rank %d's result is not byte-identical to rank 0's; "
"that is a reduce, not an all-reduce\n",
r);
++failures;
}
}
// Warm up the collective before timing it. NCCL sets up its channels
// and loads its kernels on the first call of a given shape, and that
// cost belongs to neither the link nor the algorithm.
for (int w = 0; w < kWarmupRuns; ++w) {
NCCL_CHECK(ncclGroupStart());
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
NCCL_CHECK(ncclAllReduce(d_grad[k], d_out[k], kElems, ncclFloat,
ncclSum, comms[k], streams[k]));
}
NCCL_CHECK(ncclGroupEnd());
}
for (int r = 0; r < ranks; ++r) {
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaStreamSynchronize(streams[r]));
}
// Events, not a host clock, and one pair per device, because the work
// is spread across devices and there is no single stream that sees all
// of it. The reported time is the slowest rank's, since a collective
// is not finished until every rank is.
for (int r = 0; r < ranks; ++r) {
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaEventRecord(starts[static_cast<size_t>(r)],
streams[static_cast<size_t>(r)]));
}
for (int t = 0; t < kTimedRuns; ++t) {
NCCL_CHECK(ncclGroupStart());
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
NCCL_CHECK(ncclAllReduce(d_grad[k], d_out[k], kElems, ncclFloat,
ncclSum, comms[k], streams[k]));
}
NCCL_CHECK(ncclGroupEnd());
}
for (int r = 0; r < ranks; ++r) {
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaEventRecord(stops[static_cast<size_t>(r)],
streams[static_cast<size_t>(r)]));
}
float slowestMs = 0.0f;
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaEventSynchronize(stops[k]));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, starts[k], stops[k]));
ms /= kTimedRuns;
std::printf("rank %d: %.3f ms per all-reduce\n", r, ms);
if (ms > slowestMs) {
slowestMs = ms;
}
}
// Algorithm bandwidth is buffer over time. Bus bandwidth multiplies it
// by 2(N-1)/N, the bytes a ring actually moves per rank: (N-1)/N in
// the reduce-scatter and (N-1)/N again in the all-gather. nccl-tests
// defines it as "B = S/t * (2*(n-1)/n) = algbw * (2*(n-1)/n)", and the
// point of the correction is that bus bandwidth is comparable across
// rank counts while algorithm bandwidth is not.
// snippet: bandwidth
if (ranks == 1) {
std::printf("\none-rank API path: %.3f ms\n", slowestMs);
std::printf("algorithm bandwidth: n/a, no collective traffic\n");
std::printf("bus bandwidth: n/a, 2(N-1)/N is 0 at one rank\n");
} else {
const double seconds = static_cast<double>(slowestMs) * 1.0e-3;
const double algBw = static_cast<double>(bytes) / seconds / 1.0e9;
const double factor =
2.0 * static_cast<double>(ranks - 1) / static_cast<double>(ranks);
std::printf("\nslowest rank: %.3f ms, algorithm bandwidth %.1f GB/s\n",
slowestMs, algBw);
std::printf("bus bandwidth: %.1f GB/s (2(N-1)/N = %.3f)\n",
algBw * factor, factor);
}
// end snippet
for (int r = 0; r < ranks; ++r) {
const size_t k = static_cast<size_t>(r);
CUDA_CHECK(cudaSetDevice(r));
CUDA_CHECK(cudaFree(d_grad[k]));
CUDA_CHECK(cudaFree(d_out[k]));
CUDA_CHECK(cudaEventDestroy(starts[k]));
CUDA_CHECK(cudaEventDestroy(stops[k]));
CUDA_CHECK(cudaStreamDestroy(streams[k]));
NCCL_CHECK(ncclCommDestroy(comms[k]));
}
if (failures != 0) {
std::fprintf(stderr, "%d check(s) failed\n", failures);
return EXIT_FAILURE;
}
std::printf("all %d rank(s) agree with the double-precision sum\n", ranks);
return EXIT_SUCCESS;
}