code/day59-host-threads/host_threads.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 59: host threads and the GPU.
//
// What it does: pushes 24 chunks of 8 MiB through one CUDA stream twice.
// First serially, one host thread doing fill, copy, kernel, wait per chunk.
// Then as a producer-consumer pipeline: a producer std::thread fills two
// pinned staging slots while the main thread, the consumer, owns the stream
// and enqueues the async copy and the kernel. A cudaLaunchHostFunc callback
// hands each slot back to the producer once the stream is done reading it,
// through a condition variable. Both runs must produce bit-identical output.
//
// What it does not measure: two GPUs, or two streams. One producer, one
// consumer, one stream, so the only concurrency on trial is host-side.
// Day 52 built the cross-stream dependencies; day 60 combines the two.
//
// Timing is CUDA events throughout, including the CPU-only fill leg: two
// events recorded on an idle stream complete as soon as the GPU reaches
// them, so the elapsed time between them is the host wall time between the
// two record calls, give or take a launch overhead measured in microseconds.
// That keeps every number on one clock and keeps std::chrono out.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -lineinfo -o host_threads \
// host_threads.cu
// Run: ./host_threads
#include <cmath>
#include <cstdint>
#include <cstdio>
#include <cstdlib>
#include <condition_variable>
#include <mutex>
#include <thread>
#include <vector>
#include <cuda_runtime.h>
#include <nvtx3/nvToolsExt.h>
// The one error macro. This file is standalone, so it carries its own
// verbatim copy. `err_` has a trailing underscore so it cannot collide with a
// variable at the call site.
#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)
// 24 chunks of 2^21 floats is 192 MiB end to end, with 8 MiB in flight per
// copy: big enough that the fill and the copy both take real milliseconds,
// small enough to leave a 16 GB card mostly empty. Two slots is the floor
// for overlap, same as day 54's double buffer, just with the two halves in
// two threads instead of two halves of one loop.
constexpr size_t kChunkElems = size_t{1} << 21;
constexpr int kChunks = 24;
constexpr int kSlots = 2;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kInnerIters = 128;
// x converges toward kAdd / (1 - kMul) = 1.0f, so 128 iterations neither
// overflow nor underflow whatever the fill produced.
constexpr float kMul = 0.999f;
constexpr float kAdd = 0.001f;
// The GPU may contract x * kMul + kAdd into one fma where the host compiler
// keeps two rounded steps, and over 128 iterations the difference drifts a
// few hundred ulp. 1e-4 relative absorbs that and still catches any indexing
// mistake, which lands whole chunks apart, not ulps apart.
constexpr float kRelTolerance = 1e-4f;
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kChunkElems % kThreadsPerBlock == 0,
"chunk size must fill its last block");
// What the producer writes into a staging slot. A hash of the global element
// index, so any chunk can be regenerated anywhere: the CPU reference calls
// the same function and never sees the staging buffers at all.
static float fillValue(size_t globalIndex) {
uint32_t h = static_cast<uint32_t>(globalIndex) * 2654435761u;
h ^= h >> 16;
h *= 2246822519u;
h ^= h >> 13;
return static_cast<float>(h & 0xffffffu) * (1.0f / 16777216.0f);
}
static void fillChunk(float* dst, int chunk) {
const size_t base = static_cast<size_t>(chunk) * kChunkElems;
for (size_t i = 0; i < kChunkElems; ++i) {
dst[i] = fillValue(base + i);
}
}
// snippet: kernel
__global__ void processChunk(const float* __restrict__ in,
float* __restrict__ out, size_t n) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
float x = in[i];
for (int k = 0; k < kInnerIters; ++k) {
x = x * kMul + kAdd;
}
out[i] = x;
}
}
// end snippet
static void referenceChunk(float* out, int chunk) {
const size_t base = static_cast<size_t>(chunk) * kChunkElems;
for (size_t i = 0; i < kChunkElems; ++i) {
float x = fillValue(base + i);
for (int k = 0; k < kInnerIters; ++k) {
x = x * kMul + kAdd;
}
out[i] = x;
}
}
static size_t firstMismatch(const float* got, const float* want, size_t n,
float relTolerance) {
for (size_t i = 0; i < n; ++i) {
const float scale = (want[i] == 0.0f) ? 1.0f : std::fabs(want[i]);
if (std::fabs(got[i] - want[i]) > relTolerance * scale) {
return i;
}
}
return n;
}
// The handoff. One mutex, one condition variable, one state per slot. The
// producer owns a slot while it is kFree, the consumer while it is kReady,
// the stream while it is kInFlight. Every transition happens under the lock
// and notifies, so neither thread ever spins.
// snippet: handoff
enum class SlotState { kFree, kReady, kInFlight };
// end snippet
struct Handoff {
std::mutex m;
std::condition_variable cv;
SlotState state[kSlots] = {SlotState::kFree, SlotState::kFree};
};
struct SlotDone {
Handoff* handoff;
int slot;
};
// snippet: callback
// Runs on CUDA's internal callback thread when the stream reaches it. It may
// not call any CUDA API, and it must return fast: work queued after it in
// this stream waits for it, and one callback thread can serve every stream
// in the process. Lock, flip the state, notify, get out. notify_all, not
// notify_one: the producer and the consumer wait on the same condition
// variable, and a single wakeup delivered to the wrong one is a deadlock.
static void CUDART_CB slotFree(void* userData) {
SlotDone* done = static_cast<SlotDone*>(userData);
{
std::lock_guard<std::mutex> lock(done->handoff->m);
done->handoff->state[done->slot] = SlotState::kFree;
}
done->handoff->cv.notify_all();
}
// end snippet
// The producer thread. Not one CUDA call in the loop: it fills pinned
// memory and moves slot states, nothing else. The cudaSetDevice up front is
// the per-thread habit this lesson is about. It is cheap and it binds this
// thread to device 0's process-wide context; a fresh thread on a one-GPU box
// gets device 0 anyway, and on a multi-GPU box the habit is load-bearing.
// snippet: producer
static void producerLoop(Handoff* handoff, float** h_staging) {
CUDA_CHECK(cudaSetDevice(0));
for (int chunk = 0; chunk < kChunks; ++chunk) {
const int slot = chunk % kSlots;
{
std::unique_lock<std::mutex> lock(handoff->m);
handoff->cv.wait(
lock, [&] { return handoff->state[slot] == SlotState::kFree; });
}
nvtxRangePushA("fill");
fillChunk(h_staging[slot], chunk);
nvtxRangePop();
{
std::lock_guard<std::mutex> lock(handoff->m);
handoff->state[slot] = SlotState::kReady;
}
handoff->cv.notify_all();
}
}
// end snippet
// Milliseconds between two events that bracket whatever ran on `stream`,
// including host time the stream spent waiting to be fed.
struct EventClock {
cudaEvent_t start = nullptr;
cudaEvent_t stop = nullptr;
void begin(cudaStream_t stream) {
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start, stream));
}
float endMs(cudaStream_t stream) {
CUDA_CHECK(cudaEventRecord(stop, stream));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
CUDA_CHECK(cudaEventDestroy(start));
CUDA_CHECK(cudaEventDestroy(stop));
return ms;
}
};
static void enqueueChunk(cudaStream_t stream, float* d_in, const float* h_src,
float* d_out, int chunk) {
const size_t bytes = kChunkElems * sizeof(float);
const int blocks = static_cast<int>(kChunkElems / kThreadsPerBlock);
CUDA_CHECK(
cudaMemcpyAsync(d_in, h_src, bytes, cudaMemcpyHostToDevice, stream));
processChunk<<<blocks, kThreadsPerBlock, 0, stream>>>(
d_in, d_out + static_cast<size_t>(chunk) * kChunkElems, kChunkElems);
CUDA_CHECK(cudaGetLastError());
}
int main() {
CUDA_CHECK(cudaSetDevice(0));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, 0));
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("%d chunks of %zu floats, %d staging slots\n", kChunks,
kChunkElems, kSlots);
const size_t chunkBytes = kChunkElems * sizeof(float);
const size_t totalElems = static_cast<size_t>(kChunks) * kChunkElems;
const size_t totalBytes = totalElems * sizeof(float);
float* h_staging[kSlots] = {nullptr, nullptr};
for (int s = 0; s < kSlots; ++s) {
CUDA_CHECK(cudaMallocHost(&h_staging[s], chunkBytes));
}
float* d_in[kSlots] = {nullptr, nullptr};
for (int s = 0; s < kSlots; ++s) {
CUDA_CHECK(cudaMalloc(&d_in[s], chunkBytes));
}
float* d_outSerial = nullptr;
float* d_outPipe = nullptr;
CUDA_CHECK(cudaMalloc(&d_outSerial, totalBytes));
CUDA_CHECK(cudaMalloc(&d_outPipe, totalBytes));
cudaStream_t stream;
CUDA_CHECK(cudaStreamCreate(&stream));
// Warm-up: the one kernel this program times, plus an H2D copy, so no
// timed leg below pays first-launch module loading.
fillChunk(h_staging[0], 0);
enqueueChunk(stream, d_in[0], h_staging[0], d_outSerial, 0);
CUDA_CHECK(cudaStreamSynchronize(stream));
// Leg 1, fill only: what the producer's work costs with the GPU idle.
nvtxRangePushA("fill-leg");
EventClock fillClock;
fillClock.begin(stream);
for (int chunk = 0; chunk < kChunks; ++chunk) {
fillChunk(h_staging[chunk % kSlots], chunk);
}
const float fillMs = fillClock.endMs(stream);
nvtxRangePop();
// Leg 2, GPU only: copies and kernels back to back from already-filled
// staging. Wrong data in all but two chunks, so it writes to d_outPipe,
// which the pipeline run overwrites before anything reads it.
nvtxRangePushA("gpu-leg");
EventClock gpuClock;
gpuClock.begin(stream);
for (int chunk = 0; chunk < kChunks; ++chunk) {
enqueueChunk(stream, d_in[chunk % kSlots], h_staging[chunk % kSlots],
d_outPipe, chunk);
}
const float gpuMs = gpuClock.endMs(stream);
nvtxRangePop();
// The serial baseline: one thread does everything, and the stream drains
// before each next fill because the fill reuses the slot the copy reads.
// snippet: serial
nvtxRangePushA("serial");
EventClock serialClock;
serialClock.begin(stream);
for (int chunk = 0; chunk < kChunks; ++chunk) {
const int slot = chunk % kSlots;
fillChunk(h_staging[slot], chunk);
enqueueChunk(stream, d_in[slot], h_staging[slot], d_outSerial, chunk);
CUDA_CHECK(cudaStreamSynchronize(stream));
}
const float serialMs = serialClock.endMs(stream);
nvtxRangePop();
// end snippet
// The pipeline: producer thread fills, this thread owns the stream. The
// consumer never waits on the GPU inside the loop; the callback frees
// each slot when the stream has finished with it, and the producer's
// condition-variable wait is the only backpressure.
Handoff handoff;
SlotDone done[kSlots];
for (int s = 0; s < kSlots; ++s) {
done[s] = {&handoff, s};
}
nvtxRangePushA("pipeline");
EventClock pipeClock;
pipeClock.begin(stream);
std::thread producer(producerLoop, &handoff, h_staging);
// snippet: consumer
for (int chunk = 0; chunk < kChunks; ++chunk) {
const int slot = chunk % kSlots;
{
std::unique_lock<std::mutex> lock(handoff.m);
handoff.cv.wait(
lock, [&] { return handoff.state[slot] == SlotState::kReady; });
handoff.state[slot] = SlotState::kInFlight;
}
nvtxRangePushA("enqueue");
enqueueChunk(stream, d_in[slot], h_staging[slot], d_outPipe, chunk);
CUDA_CHECK(cudaLaunchHostFunc(stream, slotFree, &done[slot]));
nvtxRangePop();
}
// end snippet
producer.join();
const float pipeMs = pipeClock.endMs(stream);
nvtxRangePop();
std::printf("\n%-22s %10s %14s\n", "leg", "total ms", "ms per chunk");
std::printf("%-22s %10.3f %14.3f\n", "fill only (CPU)", fillMs,
fillMs / kChunks);
std::printf("%-22s %10.3f %14.3f\n", "copy+kernel only", gpuMs,
gpuMs / kChunks);
std::printf("%-22s %10.3f %14.3f\n", "serial, one thread", serialMs,
serialMs / kChunks);
std::printf("%-22s %10.3f %14.3f\n", "pipeline, two threads", pipeMs,
pipeMs / kChunks);
std::printf("\nspeedup over serial: %.2fx\n", serialMs / pipeMs);
std::printf("floor (longer leg): %.3f ms\n",
fillMs > gpuMs ? fillMs : gpuMs);
std::vector<float> h_serial(totalElems);
std::vector<float> h_pipe(totalElems);
std::vector<float> h_want(kChunkElems);
CUDA_CHECK(cudaMemcpy(h_serial.data(), d_outSerial, totalBytes,
cudaMemcpyDeviceToHost));
CUDA_CHECK(cudaMemcpy(h_pipe.data(), d_outPipe, totalBytes,
cudaMemcpyDeviceToHost));
// Free before the gates, so the failure paths free too.
for (int s = 0; s < kSlots; ++s) {
CUDA_CHECK(cudaFreeHost(h_staging[s]));
CUDA_CHECK(cudaFree(d_in[s]));
}
CUDA_CHECK(cudaFree(d_outSerial));
CUDA_CHECK(cudaFree(d_outPipe));
CUDA_CHECK(cudaStreamDestroy(stream));
// Gate 1: the pipeline agrees with the CPU reference, chunk by chunk.
for (int chunk = 0; chunk < kChunks; ++chunk) {
referenceChunk(h_want.data(), chunk);
const size_t base = static_cast<size_t>(chunk) * kChunkElems;
const size_t bad = firstMismatch(h_pipe.data() + base, h_want.data(),
kChunkElems, kRelTolerance);
if (bad != kChunkElems) {
std::fprintf(stderr,
"pipeline wrong: chunk %d element %zu: got %.9g, "
"want %.9g\n",
chunk, bad, h_pipe[base + bad], h_want[bad]);
return EXIT_FAILURE;
}
}
// Gate 2: threading changed nothing. Same kernel, same inputs, so the
// serial and pipelined outputs must match to the bit, not to a tolerance.
for (size_t i = 0; i < totalElems; ++i) {
if (h_serial[i] != h_pipe[i]) {
std::fprintf(stderr,
"serial and pipeline differ at %zu: %.9g vs %.9g\n", i,
h_serial[i], h_pipe[i]);
return EXIT_FAILURE;
}
}
std::printf("\nall %d chunks match the CPU reference\n", kChunks);
std::printf("serial and pipelined outputs are bit-identical\n");
return EXIT_SUCCESS;
}