code/day10-block-size/block_size.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 10: choosing threads per block, measured on one kernel.
//
// The kernel is day 7's grayscale conversion written as day 8's grid-stride
// loop. That form is the reason this program can exist: the loop covers every
// pixel for any grid and any block size, so a sweep over block sizes changes
// one thing. A kernel sized to one thread per pixel changes its grid every
// time the block size moves, and no row of the table could then be attributed
// to either number on its own.
//
// What it does: for each block size, asks the driver how many blocks of that
// size fit on one SM, launches exactly one full wave of them, times the kernel
// with events, and checks the result against a CPU reference. Then it asks
// cudaOccupancyMaxPotentialBlockSize what it would have chosen.
//
// What it does not do. It never measures the tail effect: every row runs one
// full wave, so no row pays for a grid that is 1.1 waves deep. It sweeps the
// size of a 1D block and not the shape of a 2D one, because a 2D block changes
// the access pattern as well as the occupancy, and that is day 12's subject.
// And it says nothing about a kernel with different register pressure; the
// blocks/SM column below is this kernel's answer, not your kernel's.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o block_size block_size.cu
// Run: ./block_size
//
// Verified 2026-08-30 on a Tesla T4 (compute capability 7.5), driver
// 595.84, CUDA 12.6 (V12.6.85). Transcript: evidence/run-2026-08-30.txt
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#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)
// One 4K frame, 8,294,400 pixels. Big enough that a launch is far above the
// half-microsecond event resolution and far above launch overhead, small
// enough that four planes fit in 128 MiB of device memory.
constexpr int kWidth = 3840;
constexpr int kHeight = 2160;
constexpr size_t kPixels = static_cast<size_t>(kWidth) * kHeight;
// Rec. 601 luma weights, the ones day 7's kernel uses.
constexpr float kRedWeight = 0.299f;
constexpr float kGreenWeight = 0.587f;
constexpr float kBlueWeight = 0.114f;
// 32 on every GPU this course targets. warpSize is a run-time built-in, so it
// cannot appear in a static_assert; the compile-time copy lives here and
// main() checks the two agree on the card you ran on.
constexpr int kWarpSize = 32;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kRelTolerance = 1e-5f;
// The sweep. Four groups, and each one is here to answer a question.
//
// 16 below a warp, so half of every warp slot it takes is idle
// 32 64 the sizes the per-SM block cap reaches first on a small card
// 96 160 192 not divisors of a 32-warp SM, which is where the sawtooth is
// 100 not a multiple of 32 at all: 4 warp slots, 28 lanes idle
// 128 to 1024 the range NVIDIA's best practices guide suggests starting in
constexpr int kBlockSizes[] = {16, 32, 64, 96, 100, 128, 160,
192, 256, 384, 512, 768, 1024};
constexpr int kNumBlockSizes =
static_cast<int>(sizeof(kBlockSizes) / sizeof(kBlockSizes[0]));
// Warp slots a block of t threads occupies, which is ceil(t / 32). The guide's
// own formula, and it rounds up: 100 threads take four slots and run the
// fourth with 28 of its 32 lanes switched off.
//
// constexpr so the two claims below are checked at compile time. They follow
// only from the constants in this file, so a run-time check would look like a
// measurement of your GPU and be nothing of the kind.
constexpr int warpsPerBlock(int t) {
return (t + kWarpSize - 1) / kWarpSize;
}
constexpr int largestBlockSize() {
int largest = 0;
for (int i = 0; i < kNumBlockSizes; ++i) {
if (kBlockSizes[i] > largest) {
largest = kBlockSizes[i];
}
}
return largest;
}
static_assert(largestBlockSize() <= 1024,
"1024 threads per block is the ceiling on every compute "
"capability this course covers");
static_assert(warpsPerBlock(100) == 4,
"100 threads take four warp slots and idle 28 lanes, which is "
"the row that separates threads from warp slots");
// gray[i] = 0.299 r[i] + 0.587 g[i] + 0.114 b[i]. One thread does one pixel
// per iteration and then jumps a whole grid forward.
//
// Memory: consecutive threads take consecutive pixels, so one warp's 32
// addresses cover 128 contiguous bytes in each plane, which is four 32-byte
// sectors. That holds at every block size in the sweep, which is what makes
// the block size the only variable in the table.
//
// Launch assumption: none, and that is the point. Any grid and any block size
// covers every pixel exactly once. Day 8 built this loop; day 7 wrote the
// arithmetic inside it against a 2D grid, which taught 2D indexing but buys
// nothing here, because a grayscale conversion never looks at a neighbour.
// snippet: kernel
__global__ void grayscale(const float* r, const float* g, const float* b,
float* out, size_t n) {
const size_t step = gridDim.x * static_cast<size_t>(blockDim.x);
for (size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
i < n; i += step) {
out[i] = kRedWeight * r[i] + kGreenWeight * g[i] + kBlueWeight * b[i];
}
}
// end snippet
// CPU reference. Written for obvious correctness, not speed: plain loop, no
// OpenMP, no intrinsics. It never allocates; the caller owns every buffer.
static void grayscaleCpu(const float* r, const float* g, const float* b,
float* out, size_t n) {
for (size_t i = 0; i < n; ++i) {
out[i] = kRedWeight * r[i] + kGreenWeight * g[i] + kBlueWeight * b[i];
}
}
// 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 whole point: "wrong at 40960" names the block, "wrong" does not.
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;
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// Copy it verbatim; the alternative is six copies of the event boilerplate,
// which is how a warm-up goes missing from one of them.
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;
}
// Every row reads three planes and writes one, so the figure is computed once
// from one constant rather than per row, where the two could drift apart.
static double bandwidthGBs(float ms) {
const double bytes = 4.0 * static_cast<double>(kPixels) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// One measured configuration. Collected first and printed afterwards, because
// the "fraction of the best row" column cannot be filled until every row has
// been measured.
struct Row {
int threads;
int warps;
int idleLanes;
int blocksPerSM;
int gridBlocks;
float ms;
};
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
// Counts checks that failed, so one bad row still reports every other.
int wrong = 0;
// cudaDeviceProp carries no maxWarpsPerSM field. It reports resident
// threads per SM, and the warp count is that over the warp size.
const int maxWarpsPerSM = prop.maxThreadsPerMultiProcessor / prop.warpSize;
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf(" SMs %d, resident warps/SM %d, resident blocks/SM %d\n",
prop.multiProcessorCount, maxWarpsPerSM,
prop.maxBlocksPerMultiProcessor);
std::printf(
"Image: %d x %d, %zu pixels, %.1f MiB moved per pass\n", kWidth,
kHeight, kPixels,
4.0 * static_cast<double>(kPixels) * sizeof(float) / (1024.0 * 1024.0));
std::printf("Timing: mean of %d runs after %d warm-ups, CUDA events\n\n",
kTimedRuns, kWarmupRuns);
const size_t bytes = kPixels * sizeof(float);
std::vector<float> h_r(kPixels);
std::vector<float> h_g(kPixels);
std::vector<float> h_b(kPixels);
std::vector<float> h_out(kPixels);
std::vector<float> h_want(kPixels);
// A synthetic frame rather than a file, so the program has no data
// dependency and CI can run it anywhere. The pattern is arbitrary; the
// only property that matters is that the three planes differ, so a kernel
// that read the wrong one would be caught.
for (size_t i = 0; i < kPixels; ++i) {
h_r[i] = static_cast<float>(i % 251) / 251.0f;
h_g[i] = static_cast<float>(i % 241) / 241.0f;
h_b[i] = static_cast<float>(i % 239) / 239.0f;
}
grayscaleCpu(h_r.data(), h_g.data(), h_b.data(), h_want.data(), kPixels);
float* d_r = nullptr;
float* d_g = nullptr;
float* d_b = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_r, bytes));
CUDA_CHECK(cudaMalloc(&d_g, bytes));
CUDA_CHECK(cudaMalloc(&d_b, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMemcpy(d_r, h_r.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_g, h_g.data(), bytes, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_b, h_b.data(), bytes, cudaMemcpyHostToDevice));
Row rows[kNumBlockSizes];
for (int c = 0; c < kNumBlockSizes; ++c) {
const int t = kBlockSizes[c];
const int warps = warpsPerBlock(t);
// Blocks of this size that fit on one SM, from the driver rather than
// from arithmetic on the lesson page. It answers about the compiled
// kernel, so it accounts for a register count the source never shows.
//
// One full wave: exactly the blocks the device holds at once. The grid
// is therefore a property of the machine, not of the image, which is
// what the grid-stride loop bought. Capped at the blocks a
// one-thread-per-pixel launch would need, so a wave larger than the
// image cannot spend its extra blocks doing nothing, and floored at
// one, because a grid of zero blocks is not a launch.
// snippet: grid-policy
int blocksPerSM = 0;
CUDA_CHECK(cudaOccupancyMaxActiveBlocksPerMultiprocessor(
&blocksPerSM, grayscale, t, 0));
size_t wave = static_cast<size_t>(blocksPerSM) *
static_cast<size_t>(prop.multiProcessorCount);
const size_t needed =
(kPixels + static_cast<size_t>(t) - 1) / static_cast<size_t>(t);
if (wave > needed) {
wave = needed;
}
if (wave < 1) {
wave = 1;
}
const int gridBlocks = static_cast<int>(wave);
// end snippet
// A residency of zero means no block of this size fits at all, which
// on a real kernel means registers or shared memory, not threads. The
// floor above keeps the program running to the end so every other row
// still reports, and this branch makes sure the run still fails.
if (blocksPerSM < 1) {
std::fprintf(stderr,
"no block of %d threads fits on an SM for this "
"kernel, so its row is not a measurement\n",
t);
++wrong;
}
// The two per-SM caps the lesson states in words, checked against the
// driver's own answer rather than asserted at the reader. These are
// real branches, not asserts: CI builds Release, Release defines
// NDEBUG, and an assert under NDEBUG is deleted, so a check written
// that way would vanish in exactly the build that matters.
if (blocksPerSM > prop.maxBlocksPerMultiProcessor) {
std::fprintf(stderr,
"residency %d exceeds the per-SM block cap %d at %d "
"threads per block\n",
blocksPerSM, prop.maxBlocksPerMultiProcessor, t);
++wrong;
}
if (blocksPerSM * warps > maxWarpsPerSM) {
std::fprintf(stderr,
"resident warps %d exceeds the per-SM warp cap %d at "
"%d threads per block\n",
blocksPerSM * warps, maxWarpsPerSM, t);
++wrong;
}
// Zeroed before the run so that a configuration which covered only
// part of the image is caught by the comparison below instead of
// passing on the previous row's answer.
CUDA_CHECK(cudaMemset(d_out, 0, bytes));
const float ms = timeKernel([&] {
grayscale<<<gridBlocks, t>>>(d_r, d_g, d_b, d_out, kPixels);
});
// Correctness is re-checked at every configuration, not once. The
// claim this program rests on is that the kernel is
// launch-configuration independent, and a single check at one block
// size would not test it.
CUDA_CHECK(
cudaMemcpy(h_out.data(), d_out, bytes, cudaMemcpyDeviceToHost));
const size_t bad =
firstMismatch(h_out.data(), h_want.data(), kPixels, kRelTolerance);
if (bad != kPixels) {
std::fprintf(stderr,
"grayscale wrong at %zu with %d threads per block: "
"got %.9g, want %.9g\n",
bad, t, h_out[bad], h_want[bad]);
++wrong;
}
rows[c].threads = t;
rows[c].warps = warps;
rows[c].idleLanes = warps * kWarpSize - t;
rows[c].blocksPerSM = blocksPerSM;
rows[c].gridBlocks = gridBlocks;
rows[c].ms = ms;
}
int best = 0;
for (int c = 1; c < kNumBlockSizes; ++c) {
if (rows[c].ms < rows[best].ms) {
best = c;
}
}
std::printf("%8s %6s %5s %7s %9s %6s %7s %9s %8s %8s\n", "thr/blk", "warps",
"idle", "blk/SM", "warps/SM", "occ", "grid", "ms", "GB/s",
"of best");
int withinFive = 0;
for (int c = 0; c < kNumBlockSizes; ++c) {
const Row& row = rows[c];
const int residentWarps = row.blocksPerSM * row.warps;
const double fraction = rows[best].ms / row.ms;
if (fraction >= 0.95) {
++withinFive;
}
std::printf("%8d %6d %5d %7d %9d %5.0f%% %7d %9.3f %8.1f %8.2f\n",
row.threads, row.warps, row.idleLanes, row.blocksPerSM,
residentWarps, 100.0 * residentWarps / maxWarpsPerSM,
row.gridBlocks, row.ms, bandwidthGBs(row.ms), fraction);
}
// The API that answers the question for you. It returns the block size
// that reaches the highest occupancy this kernel can reach, and the
// smallest grid that fills the device at that size. What it maximises is
// occupancy; what the table above measures is time, and the two agree only
// when the kernel is latency-bound on warp supply.
// snippet: suggestion
int suggestedGrid = 0;
int suggestedBlock = 0;
CUDA_CHECK(cudaOccupancyMaxPotentialBlockSize(&suggestedGrid,
&suggestedBlock, grayscale));
// end snippet
std::printf(
"\ncudaOccupancyMaxPotentialBlockSize suggests %d threads per block "
"and a grid of %d blocks.\n",
suggestedBlock, suggestedGrid);
std::printf("Fastest measured: %d threads per block at %.3f ms.\n",
rows[best].threads, rows[best].ms);
std::printf(
"%d of %d configurations reach at least 95 percent of that row.\n",
withinFive, kNumBlockSizes);
// This one consults the device, so it is a runtime check.
if (prop.warpSize != kWarpSize) {
std::fprintf(stderr,
"this GPU reports warpSize %d; every number on the page "
"assumes %d\n",
prop.warpSize, kWarpSize);
++wrong;
}
CUDA_CHECK(cudaFree(d_r));
CUDA_CHECK(cudaFree(d_g));
CUDA_CHECK(cudaFree(d_b));
CUDA_CHECK(cudaFree(d_out));
if (wrong > 0) {
std::fprintf(stderr,
"\n%d check(s) failed. The table above is not a "
"measurement of what the lesson claims.\n",
wrong);
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}