code/day07-2d-grids/grayscale.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 7: a 2D grid over an image, and the two bounds checks it needs.
//
// The program builds its own test image, writes it as a binary PPM, reads it
// back with the reader below, converts it to grayscale on the GPU, checks
// every pixel against a CPU reference, and writes a binary PGM. No image
// library, and no binary committed to the repo: everything on disk is made
// here and can be deleted afterwards.
//
// Nothing is timed. Day 9 is where timing arrives, because a host clock
// around a first CUDA program measures context creation and two PCIe copies
// rather than the kernel.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o grayscale grayscale.cu
// Run: ./grayscale
//
// 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 <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)
// 611 x 397. Neither side divides its block dimension: 611 = 19 * 32 + 3 and
// 397 = 49 * 8 + 5. So the grid rounds up on both axes and both halves of the
// kernel's guard run on every launch. Pick 1024 x 1024 and neither ever
// fires, which is how a missing guard reaches somebody else's machine.
constexpr size_t kWidth = 611;
constexpr size_t kHeight = 397;
constexpr size_t kChannels = 3;
// threadIdx.x is the fastest-varying index inside a block, so 32 in x means
// one warp covers 32 consecutive columns of one row. Move the 32 to y and one
// warp spans 32 different rows instead, which is a stride of `width`.
constexpr int kBlockDimX = 32;
constexpr int kBlockDimY = 8; // 32 * 8 = 256 threads, the course default
// Fixed-point luma weights: the 8-bit form of 0.299, 0.587 and 0.114, the
// coefficients in ITU-R BT.601. https://www.itu.int/rec/R-REC-BT.601 is the
// recommendation's index page (checked 2026-08-30); the coefficients are in
// the recommendation text rather than on that page.
//
// Integers rather than floats, for one reason: a GPU may fuse a multiply and
// an add where the host compiler does not, which moves the last bit and can
// round a byte differently on the two sides. With integers the kernel and the
// CPU reference produce the same byte every time, so a mismatch in main() is
// an indexing bug and can be nothing else.
constexpr unsigned int kRedWeight = 77;
constexpr unsigned int kGreenWeight = 150;
constexpr unsigned int kBlueWeight = 29;
static_assert(kRedWeight + kGreenWeight + kBlueWeight == 256,
"the weighted sum is divided by 256 with a shift, so the "
"weights must add to 256 or the image changes brightness");
static_assert(kBlockDimX * kBlockDimY <= 1024,
"a thread block holds at most 1024 threads on every compute "
"capability this course targets");
static_assert((kBlockDimX * kBlockDimY) % 32 == 0,
"block size must be a whole number of warps");
static_assert(kWidth % kBlockDimX != 0 && kHeight % kBlockDimY != 0,
"both image sides must divide unevenly by their block "
"dimension, or the grid does not round up and neither half of "
"the bounds check ever runs");
// Integer ceiling division. constexpr so one function can size the grid at
// run time and appear in a static_assert. Day 2 shipped a static_assert over
// a function that was not constexpr and stopped compiling. The reason it is a
// static_assert and not an assert: CI builds Release, Release defines NDEBUG,
// and NDEBUG deletes assert(), so a check you can make at compile time has to
// be made at compile time.
static constexpr size_t ceilDiv(size_t a, size_t b) {
return (a + b - 1) / b;
}
static_assert(ceilDiv(kHeight, static_cast<size_t>(kBlockDimY)) <= 65535,
"gridDim.y and gridDim.z stop at 65535, unlike gridDim.x which "
"goes to 2^31 - 1. A taller image than that needs its rows "
"folded into x.");
// The luma formula, in the one place both processors can see it. __host__
// __device__ compiles it twice, once for each, so the kernel and the CPU
// reference run the same expression and the comparison in main() tests the
// indexing rather than the arithmetic.
//
// The arithmetic is exact at every input: 255 * 256 >> 8 is 255, so there is
// no rounding for the two sides to disagree about.
__host__ __device__ inline unsigned char luma(unsigned char r, unsigned char g,
unsigned char b) {
return static_cast<unsigned char>(
(kRedWeight * r + kGreenWeight * g + kBlueWeight * b) >> 8);
}
// gray[row][col] = luma of rgb[row][col]. One thread owns one pixel.
//
// Memory: threadIdx.x is the fastest-varying index inside a block, so one
// warp is 32 consecutive values of `col` inside one row. Its 32 output
// addresses are 32 consecutive bytes and its 96 input bytes are contiguous.
// Take `col` from y instead and one warp spans 32 rows, 32 addresses `width`
// bytes apart, which is the pattern day 11 measures.
//
// Launch assumption: the grid covers the image. It rounds up on both axes, so
// both halves of the guard have to be there. Dropping `col < width` does not
// fault on most runs, because row * width + col with col >= width is the
// address of pixel (row + 1, col - width), a real pixel, right up to the last
// row where it runs off the end of the buffer. README.md has the edit.
// snippet: kernel
__global__ void grayscale(const unsigned char* rgb, unsigned char* gray,
size_t width, size_t height) {
const size_t col =
blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
const size_t row =
blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
if (row < height && col < width) {
const size_t i = row * width + col;
gray[i] = luma(rgb[kChannels * i + 0], rgb[kChannels * i + 1],
rgb[kChannels * i + 2]);
}
}
// end snippet
// CPU reference. Two nested loops in row-major order, written for obvious
// correctness rather than speed: plain loops, no OpenMP, no intrinsics. It
// never allocates; the caller owns every buffer.
//
// The loop nest and the kernel's grid are two independent spellings of one
// mapping, which is the only reason comparing them proves anything.
static void grayscaleCpu(const unsigned char* rgb, unsigned char* gray,
size_t width, size_t height) {
for (size_t row = 0; row < height; ++row) {
for (size_t col = 0; col < width; ++col) {
const size_t i = row * width + col;
gray[i] = luma(rgb[kChannels * i + 0], rgb[kChannels * i + 1],
rgb[kChannels * i + 2]);
}
}
}
// Returns the first index where got and want differ, or n if they agree
// everywhere. main() turns that index back into a row and a column with / and
// %, which is the inverse of row * width + col and the only place in this
// file where the flattening runs backwards.
static size_t firstMismatch(const unsigned char* got, const unsigned char* want,
size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != want[i]) {
return i;
}
}
return n;
}
// The test image, built here so the repo carries no binary file. The three
// channels disagree with each other on purpose: red ramps along the row,
// green ramps down the column, blue is a 16-pixel checkerboard. A kernel that
// swaps row and col, or reads the wrong channel, then produces a different
// image rather than a plausible one.
static void makeTestImage(std::vector<unsigned char>* rgb, size_t width,
size_t height) {
rgb->resize(width * height * kChannels);
for (size_t row = 0; row < height; ++row) {
for (size_t col = 0; col < width; ++col) {
const size_t i = (row * width + col) * kChannels;
const size_t checker = ((row / 16) + (col / 16)) % 2;
(*rgb)[i + 0] = static_cast<unsigned char>(col % 256);
(*rgb)[i + 1] = static_cast<unsigned char>(row % 256);
(*rgb)[i + 2] = static_cast<unsigned char>(checker * 255);
}
}
}
// Netpbm, all of it. A binary PGM starts "P5" and a binary PPM starts "P6",
// then width, height and the maximum sample value as ASCII decimal separated
// by whitespace, then one whitespace character, then the raster, one byte per
// sample while the maximum is under 256. A '#' comment runs to the end of the
// line and may sit between any two header fields. Spec:
// https://netpbm.sourceforge.net/doc/ppm.html and
// https://netpbm.sourceforge.net/doc/pgm.html , both checked 2026-08-30.
//
// This is C tax, not CUDA. Read it once and never again.
// snippet: netpbm-write
static bool writeNetpbm(const char* path, const unsigned char* pixels,
size_t width, size_t height, size_t channels) {
std::FILE* file = std::fopen(path, "wb");
if (file == nullptr) {
return false;
}
std::fprintf(file, "P%d\n%zu %zu\n255\n", (channels == 1) ? 5 : 6, width,
height);
const size_t n = width * height * channels;
const bool wrote = std::fwrite(pixels, 1, n, file) == n;
return (std::fclose(file) == 0) && wrote;
}
// end snippet
// Reads one ASCII decimal header field, skipping whitespace and any '#'
// comment in front of it. It consumes one character past the digits, which on
// the last field is exactly the single whitespace the spec puts before the
// raster.
static bool readHeaderField(std::FILE* file, size_t* value) {
int c = std::fgetc(file);
while (c == '#' || c == ' ' || c == '\t' || c == '\n' || c == '\r') {
if (c == '#') {
while (c != '\n' && c != EOF) {
c = std::fgetc(file);
}
}
c = std::fgetc(file);
}
if (c < '0' || c > '9') {
return false;
}
*value = 0;
while (c >= '0' && c <= '9') {
*value = *value * 10 + static_cast<size_t>(c - '0');
c = std::fgetc(file);
}
return true;
}
// Reads a binary P5 or P6 file. Eight bits per sample only, which is what the
// writer above produces. The spec allows a maximum up to 65535 with two bytes
// per sample; this reader rejects those rather than reading them wrong.
static bool readNetpbm(const char* path, std::vector<unsigned char>* pixels,
size_t* width, size_t* height, size_t* channels) {
std::FILE* file = std::fopen(path, "rb");
if (file == nullptr) {
return false;
}
const int magic0 = std::fgetc(file);
const int magic1 = std::fgetc(file);
size_t maxValue = 0;
bool ok = magic0 == 'P' && (magic1 == '5' || magic1 == '6') &&
readHeaderField(file, width) && readHeaderField(file, height) &&
readHeaderField(file, &maxValue) && maxValue == 255;
if (ok) {
*channels = (magic1 == '5') ? 1 : 3;
pixels->resize((*width) * (*height) * (*channels));
ok = std::fread(pixels->data(), 1, pixels->size(), file) ==
pixels->size();
}
std::fclose(file);
return ok;
}
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)\n", prop.name, prop.major,
prop.minor);
const size_t pixels = kWidth * kHeight;
std::vector<unsigned char> h_rgb;
makeTestImage(&h_rgb, kWidth, kHeight);
if (!writeNetpbm("day07_input.ppm", h_rgb.data(), kWidth, kHeight,
kChannels)) {
std::fprintf(stderr, "could not write day07_input.ppm\n");
return EXIT_FAILURE;
}
// Read the file back rather than reusing the buffer already in memory.
// The reader is the half of this I/O that nothing else exercises, and a
// reader that quietly returns the wrong shape is how a 2D program ends up
// indexing a buffer whose rows are not the length it thinks they are.
std::vector<unsigned char> h_read;
size_t readWidth = 0;
size_t readHeight = 0;
size_t readChannels = 0;
if (!readNetpbm("day07_input.ppm", &h_read, &readWidth, &readHeight,
&readChannels)) {
std::fprintf(stderr, "could not read day07_input.ppm back\n");
return EXIT_FAILURE;
}
if (readWidth != kWidth || readHeight != kHeight ||
readChannels != kChannels || h_read != h_rgb) {
std::fprintf(stderr,
"PPM round trip failed: wrote %zu x %zu x %zu, read "
"%zu x %zu x %zu\n",
kWidth, kHeight, kChannels, readWidth, readHeight,
readChannels);
return EXIT_FAILURE;
}
// Every host pointer carries h_ and every device pointer d_. The image is
// one flat row-major array on both sides: nothing about `unsigned char*`
// knows it is two-dimensional, and row * width + col is the only thing
// that makes it so.
std::vector<unsigned char> h_gray(pixels);
std::vector<unsigned char> h_want(pixels);
unsigned char* d_rgb = nullptr;
unsigned char* d_gray = nullptr;
CUDA_CHECK(cudaMalloc(&d_rgb, pixels * kChannels));
CUDA_CHECK(cudaMalloc(&d_gray, pixels));
CUDA_CHECK(cudaMemcpy(d_rgb, h_read.data(), pixels * kChannels,
cudaMemcpyHostToDevice));
// dim3 holds three unsigned values and fills the ones you leave out with
// 1, so this block is 32 x 8 x 1 and this grid is 20 x 50 x 1. x carries
// the column because x is the fastest-varying dimension.
// snippet: launch
const dim3 block(kBlockDimX, kBlockDimY);
const dim3 grid(static_cast<unsigned int>(
ceilDiv(kWidth, static_cast<size_t>(kBlockDimX))),
static_cast<unsigned int>(
ceilDiv(kHeight, static_cast<size_t>(kBlockDimY))));
grayscale<<<grid, block>>>(d_rgb, d_gray, kWidth, kHeight);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
// end snippet
const size_t launched = static_cast<size_t>(grid.x) * block.x *
static_cast<size_t>(grid.y) * block.y;
std::printf("image %zu x %zu, %zu pixels\n", kWidth, kHeight, pixels);
std::printf("block (%u, %u), grid (%u, %u), %zu threads\n", block.x,
block.y, grid.x, grid.y, launched);
std::printf("%zu threads have no pixel\n", launched - pixels);
CUDA_CHECK(
cudaMemcpy(h_gray.data(), d_gray, pixels, cudaMemcpyDeviceToHost));
grayscaleCpu(h_read.data(), h_want.data(), kWidth, kHeight);
const size_t bad = firstMismatch(h_gray.data(), h_want.data(), pixels);
// Free before reporting, so the failure path frees too. Every cudaMalloc
// has a matching cudaFree before every return, including the one taken
// when the answer is wrong.
CUDA_CHECK(cudaFree(d_rgb));
CUDA_CHECK(cudaFree(d_gray));
if (bad != pixels) {
std::fprintf(stderr, "wrong at row %zu, col %zu: got %u, want %u\n",
bad / kWidth, bad % kWidth,
static_cast<unsigned int>(h_gray[bad]),
static_cast<unsigned int>(h_want[bad]));
return EXIT_FAILURE;
}
if (!writeNetpbm("day07_gray.pgm", h_gray.data(), kWidth, kHeight, 1)) {
std::fprintf(stderr, "could not write day07_gray.pgm\n");
return EXIT_FAILURE;
}
std::printf("wrote day07_input.ppm and day07_gray.pgm\n");
std::printf("all %zu pixels match the CPU reference\n", pixels);
return EXIT_SUCCESS;
}