code/day91-multi-gpu/multi_gpu.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 91: using more than one GPU.
//
// One vector add over 16,777,827 floats (64 MiB per input array), run two
// ways and priced against each other:
//
// whole array everything on device 0. One stream, one allocation set,
// one kernel. This is every program in the course so far.
// split the array cut into one shard per visible device, each
// shard with its own stream, its own allocations and its
// own copies, all in flight at once.
//
// The element count is odd, so the two shards are not the same size and
// the split arithmetic has to carry a tail. That is deliberate: an
// element count that divides by two hides the one bookkeeping bug this
// day exists to teach.
//
// THE PATH THIS PROGRAM TAKES DEPENDS ON THE MACHINE. With two or more
// devices visible it runs the two-GPU split and the peer-to-peer table.
// With one device it runs the same code with a single shard, prints that
// it took the single-GPU fallback, and skips the peer table, which needs
// a second device to mean anything. Both paths print which one ran.
//
// Correctness is exact, not approximate. a[i] is i % 1024 and b[i] is
// twice that, so the answer is 3 * (i % 1024), a whole number below 3072
// that float represents exactly, and one float add of two exact small
// integers is exact. A tolerance here would only hide a shard that
// landed at the wrong offset. The output buffer is filled with a
// sentinel first, so an element no device ever wrote fails too.
//
// What this program does not do: overlap the copies with the compute
// (day 54 does that), use more than two devices, or move any data
// between the shards. A vector add needs no communication at all, which
// is exactly why it is the honest place to start: whatever the second
// GPU wins here, it wins with zero communication, and every real
// workload pays more.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o multi_gpu multi_gpu.cu
// Run: ./multi_gpu
//
// PARTIALLY VERIFIED: the single-GPU fallback ran on a Tesla T4 on
// 2026-09-02. The required two-GPU split and peer-copy path have not run.
#include <cstdio>
#include <cstdlib>
#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)
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// Two is the number this day is about, and it is the number Kaggle's free
// T4 x2 gives. A machine with more devices still uses two, so the printed
// table means the same thing everywhere.
constexpr int kMaxDevices = 2;
// 2^24 + 611. 611 is 13 x 47, so no block size divides the count and the
// kernel's bounds check runs on every launch. The count is also odd, so
// the two shards differ in size by one element.
constexpr size_t kElems = 16777216ull + 611ull;
// The period of the input pattern. 1024 x 3 is 3072, well under 2^24, so
// every expected value is an exactly representable float.
constexpr size_t kPeriod = 1024ull;
constexpr float kSentinel = -1.0f;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kElems % 2 == 1,
"the split must carry a tail, or the lesson's bookkeeping "
"point never runs");
// out[i] = a[i] + b[i]. One thread owns one element. This is day 5's
// vector add with __restrict__ added, unchanged in every other way,
// because the point of today is where the array lives, not what the
// arithmetic is.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes, four 32-byte sectors.
//
// Launch assumption: gridDim.x * blockDim.x >= n, where n is this shard's
// element count, not the whole array's. Nothing in the kernel knows it is
// looking at a shard; the host does the arithmetic and hands over a base
// pointer.
__global__ void vectorAdd(const float* __restrict__ a,
const float* __restrict__ b, float* __restrict__ out,
size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = a[i] + b[i];
}
}
// One device's share of the work: which device, which slice of the global
// array, and the stream, events and allocations that belong to it. Nothing
// here is shared between devices, which is the whole point.
struct Shard {
int device;
size_t begin;
size_t count;
cudaStream_t stream;
cudaEvent_t start;
cudaEvent_t stop;
float* d_a;
float* d_b;
float* d_c;
};
// Cuts [0, n) into `used` shards and gives each one its own stream, event
// pair and device allocations.
//
// The tail is the part people get wrong. n / used loses the remainder, so
// the first `n % used` shards take one extra element each. With
// 16,777,827 elements over two devices that is 8,388,914 and 8,388,913,
// and the two shards are never the same size.
static void createShards(Shard* shards, int used, size_t n) {
// snippet: split
const size_t base = n / static_cast<size_t>(used);
const size_t rem = n % static_cast<size_t>(used);
size_t begin = 0;
for (int d = 0; d < used; ++d) {
shards[d].device = d;
shards[d].begin = begin;
shards[d].count = base + (static_cast<size_t>(d) < rem ? 1ull : 0ull);
begin += shards[d].count;
}
// end snippet
for (int d = 0; d < used; ++d) {
Shard& s = shards[d];
const size_t bytes = s.count * sizeof(float);
// snippet: per-device
// Everything below is created while device d is current, and
// belongs to device d for good. A stream, an event and an
// allocation are not portable objects: "a kernel launch will fail
// if it is issued to a stream that is not associated to the
// current device".
CUDA_CHECK(cudaSetDevice(s.device));
CUDA_CHECK(cudaStreamCreate(&s.stream));
CUDA_CHECK(cudaEventCreate(&s.start));
CUDA_CHECK(cudaEventCreate(&s.stop));
CUDA_CHECK(cudaMalloc(&s.d_a, bytes));
CUDA_CHECK(cudaMalloc(&s.d_b, bytes));
CUDA_CHECK(cudaMalloc(&s.d_c, bytes));
// end snippet
}
}
static void destroyShards(Shard* shards, int used) {
for (int d = 0; d < used; ++d) {
Shard& s = shards[d];
CUDA_CHECK(cudaSetDevice(s.device));
CUDA_CHECK(cudaFree(s.d_a));
CUDA_CHECK(cudaFree(s.d_b));
CUDA_CHECK(cudaFree(s.d_c));
CUDA_CHECK(cudaEventDestroy(s.start));
CUDA_CHECK(cudaEventDestroy(s.stop));
CUDA_CHECK(cudaStreamDestroy(s.stream));
}
}
// Puts one shard's work on its own stream: two host-to-device copies, the
// kernel, and one device-to-host copy back into the shard's own window of
// the host array. Nothing here synchronises, so the caller can issue the
// other device's work immediately afterwards and both run at once.
//
// The third launch argument is 0 because this kernel uses no dynamic
// shared memory; the fourth is the shard's stream, which is the argument
// that matters today.
static void enqueueShard(const Shard& s, const float* h_a, const float* h_b,
float* h_c, bool withCopies) {
const size_t bytes = s.count * sizeof(float);
const int blocks =
static_cast<int>((s.count + kThreadsPerBlock - 1) / kThreadsPerBlock);
CUDA_CHECK(cudaSetDevice(s.device));
if (withCopies) {
CUDA_CHECK(cudaMemcpyAsync(s.d_a, h_a + s.begin, bytes,
cudaMemcpyHostToDevice, s.stream));
CUDA_CHECK(cudaMemcpyAsync(s.d_b, h_b + s.begin, bytes,
cudaMemcpyHostToDevice, s.stream));
}
vectorAdd<<<blocks, kThreadsPerBlock, 0, s.stream>>>(s.d_a, s.d_b, s.d_c,
s.count);
CUDA_CHECK(cudaGetLastError());
if (withCopies) {
CUDA_CHECK(cudaMemcpyAsync(h_c + s.begin, s.d_c, bytes,
cudaMemcpyDeviceToHost, s.stream));
}
}
// Runs one path once and waits for every device to finish it.
static void runPathOnce(Shard* shards, int used, const float* h_a,
const float* h_b, float* h_c) {
for (int d = 0; d < used; ++d) {
enqueueShard(shards[d], h_a, h_b, h_c, true);
}
for (int d = 0; d < used; ++d) {
CUDA_CHECK(cudaSetDevice(shards[d].device));
CUDA_CHECK(cudaStreamSynchronize(shards[d].stream));
}
}
// Times a path and writes each device's own mean milliseconds into
// perDeviceMs. Returns the largest of them.
//
// The largest, not the sum: the devices run at the same time, so the path
// is finished when the slowest one is. It is a lower bound on wall time
// and the program says so, because each device's events only see that
// device's timeline. The host issues device 0's work first, and the few
// microseconds of skew before device 1 starts are invisible here.
//
// This is not the course's timeKernel helper. That one records on the
// legacy default stream under whatever device happens to be current, and
// both halves of that are wrong with two devices in play, since
// "cudaEventRecord() will fail if the input event and input stream are
// associated to different devices".
static float timePath(Shard* shards, int used, const float* h_a,
const float* h_b, float* h_c, bool withCopies,
float* perDeviceMs) {
for (int r = 0; r < kWarmupRuns; ++r) {
for (int d = 0; d < used; ++d) {
enqueueShard(shards[d], h_a, h_b, h_c, withCopies);
}
}
for (int d = 0; d < used; ++d) {
CUDA_CHECK(cudaSetDevice(shards[d].device));
CUDA_CHECK(cudaStreamSynchronize(shards[d].stream));
}
CUDA_CHECK(cudaGetLastError());
for (int d = 0; d < used; ++d) {
CUDA_CHECK(cudaSetDevice(shards[d].device));
CUDA_CHECK(cudaEventRecord(shards[d].start, shards[d].stream));
}
// Interleaved, so neither device waits on the host to finish issuing
// the other's whole loop.
for (int r = 0; r < kTimedRuns; ++r) {
for (int d = 0; d < used; ++d) {
enqueueShard(shards[d], h_a, h_b, h_c, withCopies);
}
}
for (int d = 0; d < used; ++d) {
CUDA_CHECK(cudaSetDevice(shards[d].device));
CUDA_CHECK(cudaEventRecord(shards[d].stop, shards[d].stream));
}
float worst = 0.0f;
for (int d = 0; d < used; ++d) {
CUDA_CHECK(cudaSetDevice(shards[d].device));
CUDA_CHECK(cudaEventSynchronize(shards[d].stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, shards[d].start, shards[d].stop));
perDeviceMs[d] = ms / kTimedRuns;
if (perDeviceMs[d] > worst) {
worst = perDeviceMs[d];
}
}
return worst;
}
// Times work issued to one stream, with that stream's own device current.
// Used only for the single-stream copies in the peer-to-peer table; the
// two-device path above cannot use it, because this helper waits for its
// own stream before returning and that would serialise the devices.
template <typename IssueFn>
static float timeStream(int device, cudaStream_t stream, cudaEvent_t start,
cudaEvent_t stop, IssueFn issue) {
CUDA_CHECK(cudaSetDevice(device));
for (int i = 0; i < kWarmupRuns; ++i) {
issue();
}
CUDA_CHECK(cudaStreamSynchronize(stream));
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaEventRecord(start, stream));
for (int i = 0; i < kTimedRuns; ++i) {
issue();
}
CUDA_CHECK(cudaEventRecord(stop, stream));
CUDA_CHECK(cudaEventSynchronize(stop));
CUDA_CHECK(cudaGetLastError());
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
return ms / kTimedRuns;
}
// GB/s for `bytes` moved in `ms` milliseconds, decimal GB as everywhere
// else in this course.
static double bandwidth(size_t bytes, float ms) {
return static_cast<double>(bytes) / (static_cast<double>(ms) * 1.0e-3) /
1.0e9;
}
// Checks the whole output array against the closed form 3 * (i % 1024).
// The expected value comes from the index, never from a[] and b[], so a
// shard written at the wrong offset fails instead of matching itself.
// Returns the number of wrong elements and reports the first one.
static size_t checkOutput(const float* h_c, size_t n, const char* label) {
size_t bad = 0;
size_t firstBad = 0;
for (size_t i = 0; i < n; ++i) {
const float want = 3.0f * static_cast<float>(i % kPeriod);
if (h_c[i] != want) {
if (bad == 0) {
firstBad = i;
}
++bad;
}
}
if (bad != 0) {
const float want = 3.0f * static_cast<float>(firstBad % kPeriod);
std::fprintf(stderr,
"%s: %zu of %zu elements wrong; first at index %zu, "
"got %.9g, want %.9g\n",
label, bad, n, firstBad, h_c[firstBad], want);
}
return bad;
}
int main() {
int deviceCount = 0;
CUDA_CHECK(cudaGetDeviceCount(&deviceCount));
if (deviceCount < 1) {
std::fprintf(stderr, "no CUDA device visible\n");
return EXIT_FAILURE;
}
const int used = (deviceCount < kMaxDevices) ? deviceCount : kMaxDevices;
for (int d = 0; d < used; ++d) {
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, d));
// The PCI address is here because the peer table below often turns
// on topology: two cards under one switch can usually reach each
// other, two under different root complexes often cannot.
std::printf(
"device %d: %s (compute capability %d.%d), %d SMs, "
"PCI %04x:%02x:%02x\n",
d, prop.name, prop.major, prop.minor, prop.multiProcessorCount,
static_cast<unsigned int>(prop.pciDomainID),
static_cast<unsigned int>(prop.pciBusID),
static_cast<unsigned int>(prop.pciDeviceID));
}
if (used >= 2) {
std::printf("path: two-GPU split (%d devices visible)\n", deviceCount);
} else {
std::printf(
"path: single-GPU fallback (1 device visible). The split runs "
"with one shard, so it does the same work as the baseline and "
"prices the extra stream and event machinery. The "
"peer-to-peer table needs a second device and is skipped.\n");
}
// Pinned host memory, because every copy below is a cudaMemcpyAsync
// and day 53 measured what happens when the source is pageable: the
// call blocks and the copy is not asynchronous at all.
const size_t bytes = kElems * sizeof(float);
float* h_a = nullptr;
float* h_b = nullptr;
float* h_c = nullptr;
CUDA_CHECK(cudaMallocHost(&h_a, bytes));
CUDA_CHECK(cudaMallocHost(&h_b, bytes));
CUDA_CHECK(cudaMallocHost(&h_c, bytes));
for (size_t i = 0; i < kElems; ++i) {
h_a[i] = static_cast<float>(i % kPeriod);
h_b[i] = 2.0f * static_cast<float>(i % kPeriod);
}
Shard whole[1];
Shard split[kMaxDevices];
createShards(whole, 1, kElems);
createShards(split, used, kElems);
std::printf("\n%zu elements, %.1f MiB per array\n", kElems,
static_cast<double>(bytes) / (1024.0 * 1024.0));
std::printf("%-14s %10s %14s %12s\n", "shard", "device", "first element",
"elements");
for (int d = 0; d < used; ++d) {
std::printf("%-14d %10d %14zu %12zu\n", d, split[d].device,
split[d].begin, split[d].count);
}
// Every failure from here on records itself and falls through to the
// one cleanup block at the bottom, so no path can return with device
// or pinned host memory still allocated.
size_t failures = 0;
for (size_t i = 0; i < kElems; ++i) {
h_c[i] = kSentinel;
}
runPathOnce(whole, 1, h_a, h_b, h_c);
failures += checkOutput(h_c, kElems, "whole array");
for (size_t i = 0; i < kElems; ++i) {
h_c[i] = kSentinel;
}
runPathOnce(split, used, h_a, h_b, h_c);
failures += checkOutput(h_c, kElems, "split");
// Two measurements per path. With copies is what the program actually
// costs; kernels only is what the arithmetic costs, and the gap
// between the two ratios is the argument of this whole lesson.
float wholeMs[1] = {0.0f};
float splitMs[kMaxDevices] = {0.0f};
const float wholeCopy = timePath(whole, 1, h_a, h_b, h_c, true, wholeMs);
const float splitCopy = timePath(split, used, h_a, h_b, h_c, true, splitMs);
float splitPer[kMaxDevices];
for (int d = 0; d < used; ++d) {
splitPer[d] = splitMs[d];
}
const float wholeKernel = timePath(whole, 1, h_a, h_b, h_c, false, wholeMs);
const float splitKernel =
timePath(split, used, h_a, h_b, h_c, false, splitMs);
// 3 arrays cross PCIe per run (two in, one out); the kernel itself
// reads two and writes one, so both rows move the same 3 x 64 MiB.
const size_t movedBytes = 3ull * bytes;
std::printf("\nmean of %d runs after %d warm-ups, slowest device shown\n",
kTimedRuns, kWarmupRuns);
std::printf("%-24s %12s %10s %10s\n", "path", "time (ms)", "GB/s",
"vs whole");
std::printf("%-24s %12.3f %10.1f %10.2f\n", "whole array, with copies",
wholeCopy, bandwidth(movedBytes, wholeCopy), 1.0);
std::printf("%-24s %12.3f %10.1f %10.2f\n", "split, with copies", splitCopy,
bandwidth(movedBytes, splitCopy), wholeCopy / splitCopy);
std::printf("%-24s %12.3f %10.1f %10.2f\n", "whole array, kernel only",
wholeKernel, bandwidth(movedBytes, wholeKernel), 1.0);
std::printf("%-24s %12.3f %10.1f %10.2f\n", "split, kernel only",
splitKernel, bandwidth(movedBytes, splitKernel),
wholeKernel / splitKernel);
for (int d = 0; d < used; ++d) {
std::printf(
" device %d alone, with copies: %.3f ms over %zu "
"elements\n",
d, splitPer[d], split[d].count);
}
if (used >= 2) {
// What the two devices can do for each other, if anything. The
// query is per ordered pair and the answer is not symmetric by
// definition: "access granted by this call is unidirectional".
// snippet: peer-query
std::printf("\npeer access, as reported by cudaDeviceCanAccessPeer\n");
int canPeer = 0;
for (int i = 0; i < used; ++i) {
for (int j = 0; j < used; ++j) {
if (i == j) {
continue;
}
int can = 0;
CUDA_CHECK(cudaDeviceCanAccessPeer(&can, i, j));
std::printf(" device %d -> device %d: %s\n", i, j,
can ? "yes" : "no");
if (i == 0 && j == 1) {
canPeer = can;
}
}
}
// end snippet
// Move shard 0's output to device 1 three ways. Every row moves
// the same bytes, so the only thing that changes is the route.
const size_t peerBytes = split[0].count * sizeof(float);
float* d_peer = nullptr;
float* h_stage = nullptr;
CUDA_CHECK(cudaSetDevice(1));
CUDA_CHECK(cudaMalloc(&d_peer, peerBytes));
CUDA_CHECK(cudaMallocHost(&h_stage, peerBytes));
// Leg 1 and leg 2 of the host-staged route are timed separately,
// each on its own device's events, because an event pair cannot
// span two devices. Their sum is what the staged copy costs.
const float d2hMs =
timeStream(0, split[0].stream, split[0].start, split[0].stop, [&] {
CUDA_CHECK(cudaMemcpyAsync(h_stage, split[0].d_c, peerBytes,
cudaMemcpyDeviceToHost,
split[0].stream));
});
const float h2dMs =
timeStream(1, split[1].stream, split[1].start, split[1].stop, [&] {
CUDA_CHECK(cudaMemcpyAsync(d_peer, h_stage, peerBytes,
cudaMemcpyHostToDevice,
split[1].stream));
});
// The peer copy before anything is enabled. It is legal either
// way; what changes is the route it takes underneath.
const float peerOffMs =
timeStream(0, split[0].stream, split[0].start, split[0].stop, [&] {
CUDA_CHECK(cudaMemcpyPeerAsync(d_peer, 1, split[0].d_c, 0,
peerBytes, split[0].stream));
});
float peerOnMs = 0.0f;
if (canPeer != 0) {
// snippet: peer-enable
CUDA_CHECK(cudaSetDevice(0));
CUDA_CHECK(cudaDeviceEnablePeerAccess(1, 0));
// end snippet
peerOnMs = timeStream(
0, split[0].stream, split[0].start, split[0].stop, [&] {
CUDA_CHECK(cudaMemcpyPeerAsync(d_peer, 1, split[0].d_c, 0,
peerBytes, split[0].stream));
});
CUDA_CHECK(cudaSetDevice(0));
CUDA_CHECK(cudaDeviceDisablePeerAccess(1));
}
std::printf("\nmoving %.1f MiB from device 0 to device 1\n",
static_cast<double>(peerBytes) / (1024.0 * 1024.0));
std::printf("%-34s %12s %10s\n", "route", "time (ms)", "GB/s");
std::printf("%-34s %12.3f %10.1f\n", "leg 1: device 0 -> pinned host",
d2hMs, bandwidth(peerBytes, d2hMs));
std::printf("%-34s %12.3f %10.1f\n", "leg 2: pinned host -> device 1",
h2dMs, bandwidth(peerBytes, h2dMs));
std::printf("%-34s %12.3f %10.1f\n", "staged total (leg 1 + leg 2)",
d2hMs + h2dMs, bandwidth(peerBytes, d2hMs + h2dMs));
std::printf("%-34s %12.3f %10.1f\n", "cudaMemcpyPeer, access off",
peerOffMs, bandwidth(peerBytes, peerOffMs));
if (canPeer != 0) {
std::printf("%-34s %12.3f %10.1f\n", "cudaMemcpyPeer, access on",
peerOnMs, bandwidth(peerBytes, peerOnMs));
} else {
std::printf("%-34s %12s %10s\n", "cudaMemcpyPeer, access on", "n/a",
"n/a");
std::printf(
" cudaDeviceCanAccessPeer said no for 0 -> 1, so "
"peer access was never enabled and every "
"cross-device byte went through host memory.\n");
}
CUDA_CHECK(cudaSetDevice(1));
CUDA_CHECK(cudaFree(d_peer));
CUDA_CHECK(cudaFreeHost(h_stage));
}
destroyShards(whole, 1);
destroyShards(split, used);
CUDA_CHECK(cudaFreeHost(h_a));
CUDA_CHECK(cudaFreeHost(h_b));
CUDA_CHECK(cudaFreeHost(h_c));
if (failures != 0) {
std::fprintf(stderr, "%zu wrong element(s) in total\n", failures);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}