code/day20-convolution/convolution.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 20, capstone 1: the image convolution pipeline.
//
// Blur a grayscale image with a 5x5 Gaussian, take the Sobel gradient
// magnitude, threshold it, and write out the coastlines as white lines on
// black. Then time every stage against this card's own copy bandwidth and
// against a single-threaded CPU doing the same three steps.
//
// It uses every day of module 2 at once: coalescing (11), a shared-memory
// tile with a halo (13), the barrier (14), the shared-memory bank map (15)
// and the filter in constant memory (18).
//
// The image is generated here and round-tripped through the day 7 PGM reader,
// so the repo carries no binary and the program needs no data directory. The
// graded capstone reads a NASA Blue Marble crop instead; see
// research/CAPSTONE-SPECS.md section 0.6.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o convolution convolution.cu
// Run: ./convolution
//
// 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`.
// snippet: check-macro
#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)
// end snippet
// ---------------------------------------------------------------------------
// Sizes and shapes.
//
// Four cases, chosen so each one can only fail for its own reason.
//
// `odd` is 611 x 613. 611 is 13 x 47 and 613 is prime, so neither side divides
// the tile and both halves of every bounds check run on every launch.
// `impulse` is a 256 x 256 field of zeros with one 1.0 in it, whose blur
// output is the filter itself. `base` is the picture the program writes out.
// `perf` is `base` tiled 4 by 4, so it is a whole number of tiles and no guard
// ever fires there. Passing `perf` and failing `odd` is the single most common
// outcome for a first tiled kernel, and that is why both are in the run.
constexpr size_t kBaseWidth = 1024;
constexpr size_t kBaseHeight = 1024;
constexpr int kPerfTiles = 4;
constexpr size_t kPerfWidth = kBaseWidth * kPerfTiles;
constexpr size_t kPerfHeight = kBaseHeight * kPerfTiles;
constexpr size_t kOddWidth = 611;
constexpr size_t kOddHeight = 613;
constexpr size_t kImpulseDim = 256;
// A 32 x 32 output tile is 1024 threads, the maximum block on every compute
// capability this course targets. On a T4 that is 32 resident warps in one
// block, which is full occupancy and also means no second block is resident to
// cover the barrier. Day 10's sweep is the tool for arguing about that; the
// exercise at the bottom of the lesson asks you to.
constexpr int kTileDim = 32;
constexpr int kBlurRadius = 2;
constexpr int kSobelRadius = 1;
constexpr int kBlurDim = 2 * kBlurRadius + 1; // 5
constexpr int kSobelDim = 2 * kSobelRadius + 1; // 3
constexpr int kBlurShared = kTileDim + 2 * kBlurRadius; // 36
constexpr int kSobelShared = kTileDim + 2 * kSobelRadius; // 34
// Two row pitches for the same shared tile, so one kernel can measure both.
// 36 is the tile as it comes; 37 is the day 15 padding fix applied to it.
constexpr int kBlurPitch = kBlurShared;
constexpr int kBlurPitchPadded = kBlurShared + 1;
constexpr int kThreadsPerBlock = 256; // 8 warps, for the 1D kernels
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
// The edge threshold, in the same units as the blurred image, which is 0 to 1.
// A coastline in the generated image is a step of about 0.45, and a 3x3 Sobel
// turns a step of d into a magnitude near 4d, so 0.30 sits well below the
// coastlines and well above the texture inside a region.
constexpr float kEdgeThreshold = 0.30f;
// Tolerances. The blur takes the pair from research/CAPSTONE-SPECS.md section
// 1.5 unchanged: 25 terms is well under the 64 the length-dependent rule
// starts at, so nothing widens.
constexpr float kRelTolerance = 1e-5f;
constexpr float kAbsTolerance = 1e-6f;
// The gradient does not, and the reason is worth four lines.
//
// gx is (a + 2b + c) - (d + 2e + f): a difference of two sums, each of order 4
// on an image that runs 0 to 1. Inside a smooth region those two sums agree to
// five figures and the answer is the small change left over, so the float
// error in the operands survives into the result at full size while the result
// itself is near zero. The error floor of a gradient is set by its inputs, not
// by its output, and grading it on a relative tolerance would fail correct
// code everywhere the picture is flat. 1e-4 is about two hundred times that
// floor and three thousand times smaller than the threshold it feeds.
constexpr float kSobelAbsTolerance = 1e-4f;
// A thresholded image is a bad comparison target: a pixel whose magnitude sits
// exactly on the threshold can land either side of it when one implementation
// accumulates in float and the other in double. So a byte mismatch is only a
// failure when the reference magnitude is further than this from the
// threshold. The float-versus-double gap on these magnitudes is around 1e-7,
// so this band is a thousand times wider than the noise it forgives.
constexpr float kEdgeBand = 1e-4f;
// ---------------------------------------------------------------------------
// The filters.
//
// The 5x5 Gaussian is the outer product of {1, 4, 6, 4, 1} with itself. Held
// as integers so the weights sum to exactly 256 and the normalisation is one
// division by a power of two, which is exact in float and leaves the
// separable and the 2D versions nothing to disagree about beyond the order of
// the adds.
constexpr int kGauss1dInt[kBlurDim] = {1, 4, 6, 4, 1};
constexpr int kGauss1dSum = 16;
constexpr int kGauss2dSum = kGauss1dSum * kGauss1dSum;
// Sobel, the usual 3x3 pair. gy is gx transposed, which is why a transposed
// filter is invisible in hypotf(gx, gy) and why the impulse case below tests
// the blur rather than the gradient.
constexpr float kSobelXHost[kSobelDim * kSobelDim] = {
-1.0f, 0.0f, 1.0f, -2.0f, 0.0f, 2.0f, -1.0f, 0.0f, 1.0f};
constexpr float kSobelYHost[kSobelDim * kSobelDim] = {
-1.0f, -2.0f, -1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 2.0f, 1.0f};
// Constant memory. Every lane of a warp reads the same coefficient at the same
// instant, which is the one access pattern constant memory is built for: the
// broadcast is served from the constant cache and costs one fetch for the
// whole warp. Day 18 measures the case where lanes read different addresses
// and constant memory loses.
__constant__ float cGauss2d[kBlurDim * kBlurDim];
__constant__ float cGauss1d[kBlurDim];
__constant__ float cSobelX[kSobelDim * kSobelDim];
__constant__ float cSobelY[kSobelDim * kSobelDim];
// Integer ceiling division. constexpr so one function can size the grid at run
// time and appear in a static_assert. The reason these are static_asserts and
// not asserts: CI builds Release, Release defines NDEBUG, and NDEBUG deletes
// assert(), so a check that can be made at compile time has to be.
static constexpr size_t ceilDiv(size_t a, size_t b) {
return (a + b - 1) / b;
}
static constexpr int gauss1dSum() {
int total = 0;
for (int k = 0; k < kBlurDim; ++k) {
total += kGauss1dInt[k];
}
return total;
}
static_assert(gauss1dSum() == kGauss1dSum,
"the 1D weights are divided by kGauss1dSum, so they must add to "
"it or the blur changes the image's brightness");
static_assert(kTileDim * kTileDim <= 1024,
"a thread block holds at most 1024 threads on every compute "
"capability this course targets");
static_assert((kTileDim * kTileDim) % 32 == 0,
"block size must be a whole number of warps");
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
// The traffic budget that rules out the naive kernel, written as a check
// rather than as prose. A 32 x 32 output tile with a 2-pixel apron reads a
// 36 x 36 square from global memory, which is 1296 / 1024 = 1.27 reads per
// output pixel. The capstone's budget is 1.6. The naive kernel reads 25. Pick
// a tile that breaks the budget and the build stops instead of the grader.
static_assert(kBlurShared * kBlurShared * 10 <= kTileDim * kTileDim * 16,
"the apron must stay under 1.6 global reads per output pixel");
// A T4 gives a block 48 KiB of shared memory by default and 64 KiB only after
// an explicit cudaFuncSetAttribute opt-in. Everything here stays under the
// default, so nothing in this file needs the opt-in.
static_assert(kBlurShared * kBlurPitchPadded * sizeof(float) <= 48 * 1024,
"the padded apron must fit the 48 KiB a T4 gives a block "
"without cudaFuncSetAttribute");
static_assert((kBlurShared * kBlurShared + kBlurShared * kTileDim) *
sizeof(float) <=
48 * 1024,
"the separable kernel's two shared arrays must fit the same "
"48 KiB");
static_assert(kSobelShared * kSobelShared * sizeof(float) <= 48 * 1024,
"the Sobel apron must fit the same 48 KiB");
static_assert(kOddWidth % static_cast<size_t>(kTileDim) != 0 &&
kOddHeight % static_cast<size_t>(kTileDim) != 0,
"the odd case must divide unevenly on both axes, or the grid "
"does not round up and neither half of the bounds check runs");
static_assert(kPerfWidth % static_cast<size_t>(kTileDim) == 0 &&
kPerfHeight % static_cast<size_t>(kTileDim) == 0,
"the timed case is deliberately a whole number of tiles, so the "
"two cases are 'guard fires' against 'guard never fires'");
static_assert(ceilDiv(kPerfHeight, static_cast<size_t>(kTileDim)) <= 65535,
"gridDim.y and gridDim.z stop at 65535, unlike gridDim.x. A "
"taller image than that needs its rows folded into x.");
// ---------------------------------------------------------------------------
// Kernels.
// out[i] = in[i]. One thread owns one element.
//
// Memory: consecutive threads take consecutive elements, so one warp's 32
// addresses cover 128 contiguous bytes and cost four 32-byte sectors, the best
// case day 11 measured.
//
// Launch assumption: gridDim.x * blockDim.x >= n. This is the baseline every
// tier is a fraction of, so it has to stay the simplest fully coalesced kernel
// in the file and nothing else.
__global__ void copyFloats(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) {
out[i] = in[i];
}
}
// out[row][col] = the 5x5 Gaussian of in around (row, col), zero outside the
// image. One thread owns one output pixel and reads 25 inputs from global
// memory.
//
// Memory: one warp is 32 consecutive values of col in one row, so each of the
// 25 loads is a contiguous 128-byte span and coalesces into four sectors. The
// pattern is not the problem here. The count is: 25 loads per output pixel
// against the 1.6 the capstone allows. L1 and L2 recover some of it, because
// neighbouring threads ask for overlapping windows, and the measurement below
// says how much.
//
// Launch assumption: a 2D grid that covers the image. It rounds up on both
// axes, so both halves of the guard have to be there.
__global__ void blurNaive(const float* __restrict__ in, float* __restrict__ out,
size_t rows, size_t cols) {
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 < rows && col < cols) {
float sum = 0.0f;
for (int p = 0; p < kBlurDim; ++p) {
for (int q = 0; q < kBlurDim; ++q) {
const long long gy =
static_cast<long long>(row) + p - kBlurRadius;
const long long gx =
static_cast<long long>(col) + q - kBlurRadius;
if (gy >= 0 && gy < static_cast<long long>(rows) && gx >= 0 &&
gx < static_cast<long long>(cols)) {
sum += cGauss2d[p * kBlurDim + q] *
in[static_cast<size_t>(gy) * cols +
static_cast<size_t>(gx)];
}
}
}
out[row * cols + col] = sum;
}
}
// Same output as blurNaive, from a shared tile with a halo.
//
// The block loads a (kTileDim + 4) square apron into shared memory, waits at
// one barrier, then every thread reads its 25 taps from there. Global traffic
// drops from 25 reads per output pixel to 1296 / 1024 = 1.27.
//
// `pitch` is the row stride of the shared tile in floats, passed in rather
// than fixed, so one kernel measures 36 and 37 with an identical instruction
// mix and only the address arithmetic differing. That is the day 11 discipline
// applied to shared memory: if you want to compare two layouts, compare them
// in one kernel.
//
// Memory: the apron load is 36 contiguous floats per row, so a warp's 32 lanes
// cover 128 contiguous bytes and coalesce. The row starts two floats before
// the aligned boundary the image row starts on, so each row's load crosses one
// more 32-byte sector than a 32-float aligned row would. On the shared side,
// lanes of a warp differ by one float in every access, so they land in 32
// different banks whatever the pitch is, which is the thing the two pitches
// exist to test.
//
// Launch assumptions: blockDim is exactly kTileDim x kTileDim, the grid covers
// the image, and the dynamic shared memory is kBlurShared * pitch floats.
__global__ void blurTiled(const float* __restrict__ in, float* __restrict__ out,
size_t rows, size_t cols, int pitch) {
extern __shared__ float tile[];
const int tx = static_cast<int>(threadIdx.x);
const int ty = static_cast<int>(threadIdx.y);
// The apron corner is two pixels above and left of the tile corner, so it
// is signed and can be negative. size_t cannot hold that, which is the one
// place in this file where an index is not a size_t.
const long long y0 =
static_cast<long long>(blockIdx.y) * kTileDim - kBlurRadius;
const long long x0 =
static_cast<long long>(blockIdx.x) * kTileDim - kBlurRadius;
// snippet: apron-load
// 1024 threads fill 1296 apron slots, so 272 of the loads are second
// helpings, and they fall on the 240 threads whose tx or ty is under 4.
// The loops step by kTileDim rather than dividing by the pitch, because a
// runtime integer division inside the load would cost more than the load,
// and would land differently on the padded row than on the unpadded one.
for (int sy = ty; sy < kBlurShared; sy += kTileDim) {
for (int sx = tx; sx < kBlurShared; sx += kTileDim) {
const long long gy = y0 + sy;
const long long gx = x0 + sx;
const bool inside = gy >= 0 && gy < static_cast<long long>(rows) &&
gx >= 0 && gx < static_cast<long long>(cols);
tile[sy * pitch + sx] = inside ? in[static_cast<size_t>(gy) * cols +
static_cast<size_t>(gx)]
: 0.0f;
}
}
__syncthreads();
// end snippet
const size_t row =
blockIdx.y * static_cast<size_t>(kTileDim) + static_cast<size_t>(ty);
const size_t col =
blockIdx.x * static_cast<size_t>(kTileDim) + static_cast<size_t>(tx);
if (row < rows && col < cols) {
float sum = 0.0f;
for (int p = 0; p < kBlurDim; ++p) {
for (int q = 0; q < kBlurDim; ++q) {
sum += cGauss2d[p * kBlurDim + q] *
tile[(ty + p) * pitch + (tx + q)];
}
}
out[row * cols + col] = sum;
}
}
// The same 5x5 Gaussian again, taken as two 5-tap passes inside one kernel.
//
// A 5x5 Gaussian is the outer product of two 5-tap filters, so blurring rows
// and then columns gives the same answer for 10 multiply-adds per output pixel
// instead of 25. Doing that as two kernels writes and re-reads a whole
// intermediate image, which on a memory-bound kernel costs more than the 15
// multiply-adds it saves. Doing it in one kernel keeps the intermediate in
// shared memory, so the global traffic is identical to blurTiled and only the
// arithmetic changes. That is the whole point, and it is the comparison the
// exercise asks you to make.
//
// Memory: identical to blurTiled, one 36 x 36 apron read per block. The second
// shared array holds the row pass, 36 rows by 32 columns, and both are read
// with lanes one float apart, so both are conflict-free.
//
// Launch assumptions: blockDim is exactly kTileDim x kTileDim and the grid
// covers the image. Two barriers, and every thread reaches both, including the
// threads whose output pixel is off the end of the image.
__global__ void blurSeparable(const float* __restrict__ in,
float* __restrict__ out, size_t rows,
size_t cols) {
__shared__ float apron[kBlurShared][kBlurShared];
__shared__ float rowPass[kBlurShared][kTileDim];
const int tx = static_cast<int>(threadIdx.x);
const int ty = static_cast<int>(threadIdx.y);
const long long y0 =
static_cast<long long>(blockIdx.y) * kTileDim - kBlurRadius;
const long long x0 =
static_cast<long long>(blockIdx.x) * kTileDim - kBlurRadius;
for (int sy = ty; sy < kBlurShared; sy += kTileDim) {
for (int sx = tx; sx < kBlurShared; sx += kTileDim) {
const long long gy = y0 + sy;
const long long gx = x0 + sx;
const bool inside = gy >= 0 && gy < static_cast<long long>(rows) &&
gx >= 0 && gx < static_cast<long long>(cols);
apron[sy][sx] = inside ? in[static_cast<size_t>(gy) * cols +
static_cast<size_t>(gx)]
: 0.0f;
}
}
__syncthreads();
// snippet: separable-passes
// The row pass runs over all 36 apron rows, because the column pass below
// needs two rows above and two below its own. It produces only 32 columns,
// which is the width of the output tile.
for (int sy = ty; sy < kBlurShared; sy += kTileDim) {
float sum = 0.0f;
for (int q = 0; q < kBlurDim; ++q) {
sum += cGauss1d[q] * apron[sy][tx + q];
}
rowPass[sy][tx] = sum;
}
__syncthreads();
const size_t row =
blockIdx.y * static_cast<size_t>(kTileDim) + static_cast<size_t>(ty);
const size_t col =
blockIdx.x * static_cast<size_t>(kTileDim) + static_cast<size_t>(tx);
if (row < rows && col < cols) {
float sum = 0.0f;
for (int p = 0; p < kBlurDim; ++p) {
sum += cGauss1d[p] * rowPass[ty + p][tx];
}
out[row * cols + col] = sum;
}
// end snippet
}
// Loads a (kTileDim + 2) square apron of `in` into `tile`, zero outside the
// image, with one thread per output pixel taking one or two slots.
//
// The caller owns the barrier that has to follow this, not the callee. A
// barrier hidden inside a helper is a barrier a reader stops accounting for,
// and the two kernels below want it in different places relative to their
// guard.
__device__ inline void loadSobelApron(const float* in,
float tile[kSobelShared][kSobelShared],
size_t rows, size_t cols) {
const int tx = static_cast<int>(threadIdx.x);
const int ty = static_cast<int>(threadIdx.y);
const long long y0 =
static_cast<long long>(blockIdx.y) * kTileDim - kSobelRadius;
const long long x0 =
static_cast<long long>(blockIdx.x) * kTileDim - kSobelRadius;
for (int sy = ty; sy < kSobelShared; sy += kTileDim) {
for (int sx = tx; sx < kSobelShared; sx += kTileDim) {
const long long gy = y0 + sy;
const long long gx = x0 + sx;
const bool inside = gy >= 0 && gy < static_cast<long long>(rows) &&
gx >= 0 && gx < static_cast<long long>(cols);
tile[sy][sx] = inside ? in[static_cast<size_t>(gy) * cols +
static_cast<size_t>(gx)]
: 0.0f;
}
}
}
// The Sobel gradient magnitude at one thread's pixel, read from a loaded
// apron. Both kernels below call this, so the two paths cannot compute a
// different number and the fused-versus-staged comparison can be exact rather
// than tolerant.
__device__ inline float sobelMagnitudeAt(float tile[kSobelShared][kSobelShared],
int ty, int tx) {
float gx = 0.0f;
float gy = 0.0f;
for (int p = 0; p < kSobelDim; ++p) {
for (int q = 0; q < kSobelDim; ++q) {
const float v = tile[ty + p][tx + q];
gx += cSobelX[p * kSobelDim + q] * v;
gy += cSobelY[p * kSobelDim + q] * v;
}
}
return hypotf(gx, gy);
}
// out[row][col] = hypotf(gx, gy) of the 3x3 Sobel pair. One thread owns one
// pixel, reading a 34 x 34 apron once per block.
//
// Memory: 34 contiguous floats per apron row, so a warp's lanes are one float
// apart on both the global read and the shared read. Reads 4 bytes and writes
// 4 bytes per output pixel, the same traffic as a copy.
//
// Launch assumptions: blockDim is exactly kTileDim x kTileDim and the grid
// covers the image. Every thread reaches the barrier.
__global__ void sobelMagnitude(const float* __restrict__ in,
float* __restrict__ out, size_t rows,
size_t cols) {
__shared__ float tile[kSobelShared][kSobelShared];
loadSobelApron(in, tile, rows, cols);
__syncthreads();
const size_t row = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
const size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
if (row < rows && col < cols) {
out[row * cols + col] = sobelMagnitudeAt(
tile, static_cast<int>(threadIdx.y), static_cast<int>(threadIdx.x));
}
}
// out[i] = in[i] > t ? 255 : 0. One thread owns one pixel.
//
// Memory: consecutive threads take consecutive elements. It reads 4 bytes and
// writes 1, so its ceiling is five bytes per pixel rather than the eight a
// copy moves, and its GB/s should be read against its own byte count and not
// against the copy row's.
//
// Launch assumption: gridDim.x * blockDim.x >= n.
__global__ void thresholdBytes(const float* __restrict__ in,
unsigned char* __restrict__ out, size_t n,
float t) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
out[i] = (in[i] > t) ? 255 : 0;
}
}
// The two stages above in one kernel: Sobel magnitude, compared to the
// threshold, one byte out.
//
// Fusing them saves a full write and a full read of a float image. The staged
// pair moves 8 bytes per pixel for the Sobel and 5 for the threshold; this
// moves 5 in total. The measurement below says whether the ratio the byte
// counts predict is the ratio the card delivers.
//
// Memory and launch assumptions are sobelMagnitude's, except that the output
// is one byte per pixel, so a warp writes 32 contiguous bytes: one sector for
// the whole warp instead of four.
// snippet: sobel-threshold-fused
__global__ void sobelThreshold(const float* __restrict__ in,
unsigned char* __restrict__ out, size_t rows,
size_t cols, float t) {
__shared__ float tile[kSobelShared][kSobelShared];
loadSobelApron(in, tile, rows, cols);
__syncthreads();
const size_t row = blockIdx.y * static_cast<size_t>(kTileDim) + threadIdx.y;
const size_t col = blockIdx.x * static_cast<size_t>(kTileDim) + threadIdx.x;
if (row < rows && col < cols) {
const float mag = sobelMagnitudeAt(tile, static_cast<int>(threadIdx.y),
static_cast<int>(threadIdx.x));
out[row * cols + col] = (mag > t) ? 255 : 0;
}
}
// end snippet
// ---------------------------------------------------------------------------
// CPU references. Written for obvious correctness, not speed: plain loops, no
// OpenMP, no blocking, no intrinsics. They never allocate; the caller owns
// every buffer. They accumulate in double where the kernels accumulate in
// float, because a reference exists to be right rather than to match bit for
// bit, and the tolerance covers the difference.
//
// These are also the denominator of the speedup the lesson prints, which makes
// them a weak baseline on purpose and says so.
static void blurCpu(const float* in, float* out, size_t rows, size_t cols) {
for (size_t row = 0; row < rows; ++row) {
for (size_t col = 0; col < cols; ++col) {
double sum = 0.0;
for (int p = 0; p < kBlurDim; ++p) {
for (int q = 0; q < kBlurDim; ++q) {
const long long gy =
static_cast<long long>(row) + p - kBlurRadius;
const long long gx =
static_cast<long long>(col) + q - kBlurRadius;
if (gy >= 0 && gy < static_cast<long long>(rows) &&
gx >= 0 && gx < static_cast<long long>(cols)) {
// The weights stay integers inside the loop and the
// one division happens after it. That is 25 fewer
// double divisions per pixel, which on 16.7 million
// pixels is the difference between a reference that
// takes seconds and one that takes a minute.
sum += static_cast<double>(kGauss1dInt[p] *
kGauss1dInt[q]) *
static_cast<double>(
in[static_cast<size_t>(gy) * cols +
static_cast<size_t>(gx)]);
}
}
}
out[row * cols + col] =
static_cast<float>(sum / static_cast<double>(kGauss2dSum));
}
}
}
static void sobelCpu(const float* in, float* out, size_t rows, size_t cols) {
for (size_t row = 0; row < rows; ++row) {
for (size_t col = 0; col < cols; ++col) {
double gx = 0.0;
double gy = 0.0;
for (int p = 0; p < kSobelDim; ++p) {
for (int q = 0; q < kSobelDim; ++q) {
const long long sy =
static_cast<long long>(row) + p - kSobelRadius;
const long long sx =
static_cast<long long>(col) + q - kSobelRadius;
if (sy >= 0 && sy < static_cast<long long>(rows) &&
sx >= 0 && sx < static_cast<long long>(cols)) {
const double v = static_cast<double>(
in[static_cast<size_t>(sy) * cols +
static_cast<size_t>(sx)]);
gx += static_cast<double>(
kSobelXHost[p * kSobelDim + q]) *
v;
gy += static_cast<double>(
kSobelYHost[p * kSobelDim + q]) *
v;
}
}
}
out[row * cols + col] = static_cast<float>(std::hypot(gx, gy));
}
}
}
static void thresholdCpu(const float* in, unsigned char* out, size_t n,
float t) {
for (size_t i = 0; i < n; ++i) {
out[i] = (in[i] > t) ? 255 : 0;
}
}
// ---------------------------------------------------------------------------
// Comparison helpers. Every one of them returns an index rather than a bool,
// because "wrong at row 610, col 0" names the tile and "wrong" does not.
static size_t firstMismatch(const float* got, const float* want, size_t n,
float absTolerance) {
for (size_t i = 0; i < n; ++i) {
const float diff = std::fabs(got[i] - want[i]);
if (diff > absTolerance + kRelTolerance * std::fabs(want[i])) {
return i;
}
}
return n;
}
// A byte mismatch counts only where the reference magnitude is further than
// kEdgeBand from the threshold. Inside the band the two implementations are
// allowed to disagree, and the caller prints how many pixels that covered so
// the allowance is visible rather than silent.
static size_t firstEdgeMismatch(const unsigned char* got,
const unsigned char* want, const float* mag,
size_t n, size_t* forgiven) {
*forgiven = 0;
size_t bad = n;
for (size_t i = 0; i < n; ++i) {
if (got[i] == want[i]) {
continue;
}
if (std::fabs(mag[i] - kEdgeThreshold) <= kEdgeBand) {
++(*forgiven);
continue;
}
if (bad == n) {
bad = i;
}
}
return bad;
}
static size_t firstByteMismatch(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;
}
// ---------------------------------------------------------------------------
// Timing.
//
// Copied verbatim from day 9, which is the point: it is the one template and
// the one lambda this course allows in module 1 to 3 code, and the warm-up
// lives inside it so it cannot go missing from the fourth kernel somebody adds
// later. Lazy module loading has been the default since CUDA 12.2 on Linux, so
// the first launch of each kernel pays its own load and every kernel timed
// here is warmed by this function before it is timed.
// snippet: time-kernel
template <typename LaunchFn>
static float timeKernel(LaunchFn launch) {
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
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;
}
// end snippet
// Effective bandwidth in GB/s from the bytes a stage must move and its time.
// The byte count is passed in and printed, so every row of the table can be
// checked against its own arithmetic rather than trusted. Comparing two rows
// that moved different byte counts is the bug day 11 exists to correct.
static double gigabytesPerSecond(double movedBytes, float ms) {
return movedBytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
static void printRow(const char* label, float ms, double movedBytes,
double copyGbs) {
const double gbs = gigabytesPerSecond(movedBytes, ms);
std::printf("%-34s %9.3f %9.1f %13.0f %8.2f\n", label, ms, gbs, movedBytes,
gbs / copyGbs);
}
// ---------------------------------------------------------------------------
// Netpbm, all of it. Taken unchanged from code/day07-2d-grids/grayscale.cu,
// which is the whole reason it is here: the reader already exists, it has run
// on this hardware, and writing a second one would be twenty more lines of
// file format and no CUDA.
//
// 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.
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;
}
// 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;
}
// ---------------------------------------------------------------------------
// The test image, built here so the repo carries no binary file.
//
// A few sinusoids added together make blobby regions with long curved
// boundaries; the sign of that sum is land or sea and a finer term gives each
// region its own texture. It is not a photograph and it does not need to be.
// What it needs is boundaries that are long, curved and not axis-aligned, so
// that a transposed filter or an apron that is off by one changes the picture
// in a way you can see without a diff.
//
// The graded capstone uses a 1024 x 1024 crop of NASA's Blue Marble Next
// Generation composite instead, shipped under research/CAPSTONE-SPECS.md
// section 0.6. This program generates its input so that it runs anywhere with
// no data directory and no download.
static void makeEarthImage(std::vector<unsigned char>* image, size_t width,
size_t height) {
image->resize(width * height);
for (size_t row = 0; row < height; ++row) {
for (size_t col = 0; col < width; ++col) {
const double y =
static_cast<double>(row) / static_cast<double>(height);
const double x =
static_cast<double>(col) / static_cast<double>(width);
const double land =
std::sin(6.1 * x + 0.7) * std::cos(4.3 * y - 1.1) +
0.62 * std::sin(11.7 * x - 2.0) * std::cos(9.1 * y + 0.4) +
0.31 * std::sin(19.3 * (x + 0.6 * y)) +
0.18 * std::cos(27.5 * (y - 0.4 * x));
const double grain =
0.5 + 0.5 * std::sin(83.0 * x + 51.0 * y) * std::cos(67.0 * y);
const double value =
(land > 0.0) ? (0.62 + 0.16 * grain) : (0.10 + 0.10 * grain);
(*image)[row * width + col] =
static_cast<unsigned char>(value * 255.0 + 0.5);
}
}
}
// Repeats a base image `times` by `times`. The seams are edges too, and the
// capstone's timed case is defined this way so that one megabyte on disk
// becomes a 67 megabyte buffer in memory.
static void tileImage(const unsigned char* base, size_t baseWidth,
size_t baseHeight, int times,
std::vector<unsigned char>* out) {
const size_t width = baseWidth * static_cast<size_t>(times);
const size_t height = baseHeight * static_cast<size_t>(times);
out->resize(width * height);
for (size_t row = 0; row < height; ++row) {
const size_t srcRow = row % baseHeight;
for (size_t col = 0; col < width; ++col) {
(*out)[row * width + col] =
base[srcRow * baseWidth + (col % baseWidth)];
}
}
}
static void toFloat(const unsigned char* bytes, std::vector<float>* values,
size_t n) {
values->resize(n);
for (size_t i = 0; i < n; ++i) {
(*values)[i] = static_cast<float>(bytes[i]) / 255.0f;
}
}
static void toBytes(const float* values, std::vector<unsigned char>* bytes,
size_t n) {
bytes->resize(n);
for (size_t i = 0; i < n; ++i) {
const float clamped =
(values[i] < 0.0f) ? 0.0f : ((values[i] > 1.0f) ? 1.0f : values[i]);
(*bytes)[i] = static_cast<unsigned char>(clamped * 255.0f + 0.5f);
}
}
// ---------------------------------------------------------------------------
// Launch helpers, so the grid arithmetic appears once.
static dim3 tileGrid(size_t rows, size_t cols) {
return dim3(
static_cast<unsigned int>(ceilDiv(cols, static_cast<size_t>(kTileDim))),
static_cast<unsigned int>(
ceilDiv(rows, static_cast<size_t>(kTileDim))));
}
static int linearBlocks(size_t n) {
return static_cast<int>(ceilDiv(n, static_cast<size_t>(kThreadsPerBlock)));
}
// The whole pipeline, as one call, the way the capstone's `solve_pipeline`
// contract asks for it: blur into a scratch buffer, then Sobel and threshold
// fused into the output. Two launches, not three, and no host round trip
// between them.
static void runPipeline(const float* d_in, float* d_scratch,
unsigned char* d_out, size_t rows, size_t cols) {
const dim3 block(kTileDim, kTileDim);
const dim3 grid = tileGrid(rows, cols);
blurSeparable<<<grid, block>>>(d_in, d_scratch, rows, cols);
sobelThreshold<<<grid, block>>>(d_scratch, d_out, rows, cols,
kEdgeThreshold);
}
// ---------------------------------------------------------------------------
// Correctness.
// Runs every kernel over one image and compares it to the CPU reference.
// Returns true when all of them agree.
//
// The blurred and thresholded reference buffers are also handed back, because
// the timed case reuses them and re-running a 16.7 megapixel double-precision
// reference to get them twice would double the program's run time for nothing.
static bool checkCase(const char* name, const float* h_in, size_t rows,
size_t cols, float* d_in, float* d_blur, float* d_mag,
unsigned char* d_bytes, std::vector<float>* h_blurWant,
std::vector<float>* h_magWant,
std::vector<unsigned char>* h_edgeWant) {
const size_t pixels = rows * cols;
const size_t floatBytes = pixels * sizeof(float);
const dim3 block(kTileDim, kTileDim);
const dim3 grid = tileGrid(rows, cols);
std::printf("case %-8s %zu x %zu, grid (%u, %u), block (32, 32)\n", name,
rows, cols, grid.x, grid.y);
h_blurWant->resize(pixels);
h_magWant->resize(pixels);
h_edgeWant->resize(pixels);
blurCpu(h_in, h_blurWant->data(), rows, cols);
sobelCpu(h_blurWant->data(), h_magWant->data(), rows, cols);
thresholdCpu(h_magWant->data(), h_edgeWant->data(), pixels, kEdgeThreshold);
CUDA_CHECK(cudaMemcpy(d_in, h_in, floatBytes, cudaMemcpyHostToDevice));
std::vector<float> h_got(pixels);
std::vector<unsigned char> h_gotBytes(pixels);
// Four spellings of the same blur, all compared against the one
// reference. Two of them differ only in the pitch of the shared tile.
const char* blurNames[4] = {"blur naive", "blur tiled, pitch 36",
"blur tiled, pitch 37", "blur separable"};
for (int variant = 0; variant < 4; ++variant) {
if (variant == 0) {
blurNaive<<<grid, block>>>(d_in, d_blur, rows, cols);
} else if (variant == 3) {
blurSeparable<<<grid, block>>>(d_in, d_blur, rows, cols);
} else {
const int pitch = (variant == 1) ? kBlurPitch : kBlurPitchPadded;
const size_t sharedBytes =
static_cast<size_t>(kBlurShared) * pitch * sizeof(float);
blurTiled<<<grid, block, sharedBytes>>>(d_in, d_blur, rows, cols,
pitch);
}
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_got.data(), d_blur, floatBytes,
cudaMemcpyDeviceToHost));
const size_t bad = firstMismatch(h_got.data(), h_blurWant->data(),
pixels, kAbsTolerance);
if (bad != pixels) {
std::fprintf(stderr,
"%s wrong on %s at row %zu, col %zu: got %.9g, "
"want %.9g\n",
blurNames[variant], name, bad / cols, bad % cols,
static_cast<double>(h_got[bad]),
static_cast<double>((*h_blurWant)[bad]));
return false;
}
}
// From here on the blur output on the device is the separable one, which
// is what the pipeline uses, so the Sobel below is fed the same bytes the
// pipeline feeds it.
sobelMagnitude<<<grid, block>>>(d_blur, d_mag, rows, cols);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_got.data(), d_mag, floatBytes, cudaMemcpyDeviceToHost));
const size_t badMag = firstMismatch(h_got.data(), h_magWant->data(), pixels,
kSobelAbsTolerance);
if (badMag != pixels) {
std::fprintf(stderr,
"sobel wrong on %s at row %zu, col %zu: got %.9g, "
"want %.9g\n",
name, badMag / cols, badMag % cols,
static_cast<double>(h_got[badMag]),
static_cast<double>((*h_magWant)[badMag]));
return false;
}
thresholdBytes<<<linearBlocks(pixels), kThreadsPerBlock>>>(
d_mag, d_bytes, pixels, kEdgeThreshold);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_gotBytes.data(), d_bytes, pixels, cudaMemcpyDeviceToHost));
size_t forgiven = 0;
const size_t badEdge =
firstEdgeMismatch(h_gotBytes.data(), h_edgeWant->data(),
h_magWant->data(), pixels, &forgiven);
if (badEdge != pixels) {
std::fprintf(stderr,
"threshold wrong on %s at row %zu, col %zu: got %u, "
"want %u, reference magnitude %.9g\n",
name, badEdge / cols, badEdge % cols,
static_cast<unsigned int>(h_gotBytes[badEdge]),
static_cast<unsigned int>((*h_edgeWant)[badEdge]),
static_cast<double>((*h_magWant)[badEdge]));
return false;
}
// The fused kernel and the staged pair run the same device function over
// the same inputs, so this comparison is exact rather than tolerant. A
// difference here is a bug in the fusion and can be nothing else.
std::vector<unsigned char> h_fused(pixels);
sobelThreshold<<<grid, block>>>(d_blur, d_bytes, rows, cols,
kEdgeThreshold);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_fused.data(), d_bytes, pixels, cudaMemcpyDeviceToHost));
const size_t badFused =
firstByteMismatch(h_fused.data(), h_gotBytes.data(), pixels);
if (badFused != pixels) {
std::fprintf(stderr,
"fused sobel+threshold differs from the staged pair on "
"%s at row %zu, col %zu: got %u, want %u\n",
name, badFused / cols, badFused % cols,
static_cast<unsigned int>(h_fused[badFused]),
static_cast<unsigned int>(h_gotBytes[badFused]));
return false;
}
std::printf(
" every kernel matches the CPU reference at %zu pixels, "
"%zu pixels forgiven inside the threshold band\n",
pixels, forgiven);
return true;
}
// The impulse response, which is the one check the CPU reference cannot give
// you, because it compares against the filter constants rather than against
// another implementation of the same arithmetic.
//
// Blur an image that is zero everywhere except one pixel and the output is the
// filter itself, centred on that pixel. Three things are checked: the peak
// sits where the impulse sat, the 25 values around it equal the 25
// coefficients, and the whole image sums to 1. A dropped apron tap fails the
// sum, an apron loaded from the wrong corner moves the peak, and a filter
// indexed the wrong way round fails the 25.
//
// What it cannot catch: a transposed Gaussian. This filter is symmetric, so it
// is its own transpose, and the same is true of the Sobel pair once the two
// gradients go through hypotf. Say that rather than claim the check is
// stronger than it is.
static bool checkImpulse(float* d_in, float* d_blur) {
const size_t rows = kImpulseDim;
const size_t cols = kImpulseDim;
const size_t pixels = rows * cols;
const size_t floatBytes = pixels * sizeof(float);
const size_t centreRow = rows / 2;
const size_t centreCol = cols / 2;
std::vector<float> h_in(pixels, 0.0f);
h_in[centreRow * cols + centreCol] = 1.0f;
CUDA_CHECK(
cudaMemcpy(d_in, h_in.data(), floatBytes, cudaMemcpyHostToDevice));
const dim3 block(kTileDim, kTileDim);
const dim3 grid = tileGrid(rows, cols);
const size_t sharedBytes =
static_cast<size_t>(kBlurShared) * kBlurPitch * sizeof(float);
blurTiled<<<grid, block, sharedBytes>>>(d_in, d_blur, rows, cols,
kBlurPitch);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<float> h_got(pixels);
CUDA_CHECK(
cudaMemcpy(h_got.data(), d_blur, floatBytes, cudaMemcpyDeviceToHost));
double total = 0.0;
size_t peak = 0;
for (size_t i = 0; i < pixels; ++i) {
total += static_cast<double>(h_got[i]);
if (h_got[i] > h_got[peak]) {
peak = i;
}
}
if (peak != centreRow * cols + centreCol) {
std::fprintf(stderr,
"impulse: the blur peak is at row %zu, col %zu, not at "
"row %zu, col %zu. The apron is loaded from the wrong "
"corner.\n",
peak / cols, peak % cols, centreRow, centreCol);
return false;
}
if (std::fabs(total - 1.0) > 1e-4) {
std::fprintf(stderr,
"impulse: the blurred image sums to %.9g, not 1. A tap "
"is missing from the apron or counted twice.\n",
total);
return false;
}
const size_t top = centreRow - static_cast<size_t>(kBlurRadius);
const size_t left = centreCol - static_cast<size_t>(kBlurRadius);
for (int p = 0; p < kBlurDim; ++p) {
for (int q = 0; q < kBlurDim; ++q) {
const size_t row = top + static_cast<size_t>(p);
const size_t col = left + static_cast<size_t>(q);
const float want =
static_cast<float>(kGauss1dInt[p] * kGauss1dInt[q]) /
static_cast<float>(kGauss2dSum);
const float got = h_got[row * cols + col];
if (std::fabs(got - want) > kAbsTolerance) {
std::fprintf(stderr,
"impulse: response at offset (%d, %d) is %.9g, "
"want %.9g. The filter is indexed wrongly.\n",
p - kBlurRadius, q - kBlurRadius,
static_cast<double>(got),
static_cast<double>(want));
return false;
}
}
}
std::printf(
"case impulse %zu x %zu, peak in place, sums to %.6f, all "
"25 taps equal the filter\n",
rows, cols, total);
return true;
}
// ---------------------------------------------------------------------------
// The timed table.
static bool timeStages(const float* h_in, const std::vector<float>& h_magWant,
const std::vector<unsigned char>& h_edgeWant,
size_t rows, size_t cols, float* d_in, float* d_blur,
float* d_mag, unsigned char* d_bytes) {
const size_t pixels = rows * cols;
const size_t floatBytes = pixels * sizeof(float);
const dim3 block(kTileDim, kTileDim);
const dim3 grid = tileGrid(rows, cols);
const int blocks1d = linearBlocks(pixels);
CUDA_CHECK(cudaMemcpy(d_in, h_in, floatBytes, cudaMemcpyHostToDevice));
// Byte counts, stated once, printed with every row. A stage's effective
// bandwidth is the bytes it must touch over its time, and two rows with
// different byte counts are not comparable in milliseconds.
const double fourBytes = static_cast<double>(pixels) * sizeof(float);
const double oneByte = static_cast<double>(pixels);
const double copyBytes = 2.0 * fourBytes;
const double blurBytes = 2.0 * fourBytes;
const double sobelBytes = 2.0 * fourBytes;
const double thresholdBytesMoved = fourBytes + oneByte;
const double stagedBytes = sobelBytes + thresholdBytesMoved;
const double fusedBytes = fourBytes + oneByte;
const double pipelineBytes = blurBytes + fusedBytes;
const size_t sharedPlain =
static_cast<size_t>(kBlurShared) * kBlurPitch * sizeof(float);
const size_t sharedPadded =
static_cast<size_t>(kBlurShared) * kBlurPitchPadded * sizeof(float);
// The copy baseline writes into d_mag, which nothing has filled yet at
// this point in the run. Every tier below is a fraction of this row, so it
// is measured first, on the same buffers, seconds before the kernels it
// grades.
const float msCopy = timeKernel([&] {
copyFloats<<<blocks1d, kThreadsPerBlock>>>(d_in, d_mag, pixels);
});
const double copyGbs = gigabytesPerSecond(copyBytes, msCopy);
const float msNaive = timeKernel(
[&] { blurNaive<<<grid, block>>>(d_in, d_blur, rows, cols); });
const float msTiled = timeKernel([&] {
blurTiled<<<grid, block, sharedPlain>>>(d_in, d_blur, rows, cols,
kBlurPitch);
});
const float msPadded = timeKernel([&] {
blurTiled<<<grid, block, sharedPadded>>>(d_in, d_blur, rows, cols,
kBlurPitchPadded);
});
const float msSeparable = timeKernel(
[&] { blurSeparable<<<grid, block>>>(d_in, d_blur, rows, cols); });
// d_blur now holds the separable blur of the input, which is what the
// Sobel rows below read.
const float msSobel = timeKernel(
[&] { sobelMagnitude<<<grid, block>>>(d_blur, d_mag, rows, cols); });
const float msThreshold = timeKernel([&] {
thresholdBytes<<<blocks1d, kThreadsPerBlock>>>(d_mag, d_bytes, pixels,
kEdgeThreshold);
});
const float msFused = timeKernel([&] {
sobelThreshold<<<grid, block>>>(d_blur, d_bytes, rows, cols,
kEdgeThreshold);
});
const float msPipeline =
timeKernel([&] { runPipeline(d_in, d_blur, d_bytes, rows, cols); });
// The CPU reference, timed with the same clock as everything above. The
// two events sit on an idle stream, so each record completes as soon as
// the host submits it and the gap between them is the host work. It is a
// host measurement taken on the device clock, which is what keeps every
// row of this table on one clock and is why the site's ban on std::chrono
// in a .cu file costs nothing. The reference has already run once during
// the correctness pass, so this pass is warm.
std::vector<float> h_cpuBlur(pixels);
std::vector<float> h_cpuMag(pixels);
std::vector<unsigned char> h_cpuEdge(pixels);
cudaEvent_t cpuStart, cpuStop;
CUDA_CHECK(cudaEventCreate(&cpuStart));
CUDA_CHECK(cudaEventCreate(&cpuStop));
CUDA_CHECK(cudaEventRecord(cpuStart));
blurCpu(h_in, h_cpuBlur.data(), rows, cols);
sobelCpu(h_cpuBlur.data(), h_cpuMag.data(), rows, cols);
thresholdCpu(h_cpuMag.data(), h_cpuEdge.data(), pixels, kEdgeThreshold);
CUDA_CHECK(cudaEventRecord(cpuStop));
CUDA_CHECK(cudaEventSynchronize(cpuStop));
float msCpu = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&msCpu, cpuStart, cpuStop));
CUDA_CHECK(cudaEventDestroy(cpuStart));
CUDA_CHECK(cudaEventDestroy(cpuStop));
std::printf(
"\ntimed on %zu x %zu, mean of %d runs after %d warm-ups, "
"copies not included\n",
rows, cols, kTimedRuns, kWarmupRuns);
std::printf("%-34s %9s %9s %13s %8s\n", "stage", "ms", "GB/s", "bytes",
"x copy");
printRow("copy (the baseline)", msCopy, copyBytes, copyGbs);
printRow("blur, naive", msNaive, blurBytes, copyGbs);
printRow("blur, tiled, pitch 36", msTiled, blurBytes, copyGbs);
printRow("blur, tiled, pitch 37 (padded)", msPadded, blurBytes, copyGbs);
printRow("blur, separable, one kernel", msSeparable, blurBytes, copyGbs);
printRow("sobel magnitude", msSobel, sobelBytes, copyGbs);
printRow("threshold", msThreshold, thresholdBytesMoved, copyGbs);
printRow("sobel + threshold, two kernels", msSobel + msThreshold,
stagedBytes, copyGbs);
printRow("sobel + threshold, fused", msFused, fusedBytes, copyGbs);
printRow("pipeline, blur + fused", msPipeline, pipelineBytes, copyGbs);
std::printf("%-34s %9.3f %9s %13s %8s\n", "CPU, one thread, whole pipeline",
msCpu, "n/a", "n/a", "n/a");
std::printf("\nCPU / GPU pipeline speedup: %.1fx\n",
static_cast<double>(msCpu) / static_cast<double>(msPipeline));
std::printf(
"pipeline at %.1f%% of this card's copy bandwidth\n",
100.0 * gigabytesPerSecond(pipelineBytes, msPipeline) / copyGbs);
std::printf(
"padding the shared tile changed the blur by %.2f%%\n",
100.0 * (static_cast<double>(msTiled) - static_cast<double>(msPadded)) /
static_cast<double>(msTiled));
// The last launch left the fused pipeline's bytes on the device. Check
// them against the reference that the correctness pass already built, so
// the timed path is not an untested path. This is a real branch returning
// false, not an assert: CI builds Release, Release defines NDEBUG, and an
// assert under NDEBUG is deleted.
std::vector<unsigned char> h_gotBytes(pixels);
CUDA_CHECK(
cudaMemcpy(h_gotBytes.data(), d_bytes, pixels, cudaMemcpyDeviceToHost));
size_t forgiven = 0;
const size_t bad = firstEdgeMismatch(h_gotBytes.data(), h_edgeWant.data(),
h_magWant.data(), pixels, &forgiven);
if (bad != pixels) {
std::fprintf(stderr,
"the timed pipeline is wrong at row %zu, col %zu: got "
"%u, want %u\n",
bad / cols, bad % cols,
static_cast<unsigned int>(h_gotBytes[bad]),
static_cast<unsigned int>(h_edgeWant[bad]));
return false;
}
// Two containment gates, and neither is a claim about your card. The
// pipeline launches blurSeparable and then sobelThreshold on one stream,
// so it cannot be faster than either of them measured alone. If one of
// these fires, the timing is broken and every number above it is
// worthless. Whether fusing beats the staged pair, and by how much, is a
// prediction the table answers rather than a gate; see README.md.
if (msPipeline < msFused) {
std::fprintf(stderr,
"the pipeline (%.3f ms) is faster than the fused edge "
"pass it contains (%.3f ms). The timing is broken.\n",
static_cast<double>(msPipeline),
static_cast<double>(msFused));
return false;
}
if (msPipeline < msSeparable) {
std::fprintf(stderr,
"the pipeline (%.3f ms) is faster than the blur it "
"contains (%.3f ms). The timing is broken.\n",
static_cast<double>(msPipeline),
static_cast<double>(msSeparable));
return false;
}
return true;
}
// ---------------------------------------------------------------------------
// Everything between the allocation and the free, so main() has exactly one
// cleanup block and every path, including every failure, falls through it.
static int runEverything(const std::vector<unsigned char>& h_baseBytes,
float* d_in, float* d_blur, float* d_mag,
unsigned char* d_bytes) {
std::vector<float> h_blurWant;
std::vector<float> h_magWant;
std::vector<unsigned char> h_edgeWant;
// Case 1, the odd size. Neither side divides the tile, so both halves of
// every bounds check run and the last tile in each direction is partly
// outside the image.
std::vector<unsigned char> h_oddBytes;
std::vector<float> h_odd;
makeEarthImage(&h_oddBytes, kOddWidth, kOddHeight);
toFloat(h_oddBytes.data(), &h_odd, kOddWidth * kOddHeight);
if (!checkCase("odd", h_odd.data(), kOddHeight, kOddWidth, d_in, d_blur,
d_mag, d_bytes, &h_blurWant, &h_magWant, &h_edgeWant)) {
return EXIT_FAILURE;
}
// Case 2, the impulse.
if (!checkImpulse(d_in, d_blur)) {
return EXIT_FAILURE;
}
// Case 3, the picture. This is the deliverable: run the pipeline over the
// base image and write the edges out as a PGM you can open.
const size_t basePixels = kBaseWidth * kBaseHeight;
std::vector<float> h_base;
toFloat(h_baseBytes.data(), &h_base, basePixels);
CUDA_CHECK(cudaMemcpy(d_in, h_base.data(), basePixels * sizeof(float),
cudaMemcpyHostToDevice));
runPipeline(d_in, d_blur, d_bytes, kBaseHeight, kBaseWidth);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
std::vector<unsigned char> h_edges(basePixels);
CUDA_CHECK(cudaMemcpy(h_edges.data(), d_bytes, basePixels,
cudaMemcpyDeviceToHost));
std::vector<float> h_blurred(basePixels);
CUDA_CHECK(cudaMemcpy(h_blurred.data(), d_blur, basePixels * sizeof(float),
cudaMemcpyDeviceToHost));
std::vector<unsigned char> h_blurredBytes;
toBytes(h_blurred.data(), &h_blurredBytes, basePixels);
if (!writeNetpbm("day20_blur.pgm", h_blurredBytes.data(), kBaseWidth,
kBaseHeight, 1) ||
!writeNetpbm("day20_edges.pgm", h_edges.data(), kBaseWidth, kBaseHeight,
1)) {
std::fprintf(stderr, "could not write the output images\n");
return EXIT_FAILURE;
}
size_t white = 0;
for (size_t i = 0; i < basePixels; ++i) {
white += (h_edges[i] != 0) ? 1 : 0;
}
std::printf(
"wrote day20_blur.pgm and day20_edges.pgm, %.2f%% of the "
"edge image is white\n",
100.0 * static_cast<double>(white) / static_cast<double>(basePixels));
// Case 4, the timed size. Correctness first, on the same buffers, and the
// reference it builds is handed to the timing pass so the CPU work is not
// repeated a third time.
std::vector<unsigned char> h_perfBytes;
std::vector<float> h_perf;
tileImage(h_baseBytes.data(), kBaseWidth, kBaseHeight, kPerfTiles,
&h_perfBytes);
toFloat(h_perfBytes.data(), &h_perf, kPerfWidth * kPerfHeight);
if (!checkCase("perf", h_perf.data(), kPerfHeight, kPerfWidth, d_in, d_blur,
d_mag, d_bytes, &h_blurWant, &h_magWant, &h_edgeWant)) {
return EXIT_FAILURE;
}
if (!timeStages(h_perf.data(), h_magWant, h_edgeWant, kPerfHeight,
kPerfWidth, d_in, d_blur, d_mag, d_bytes)) {
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
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);
// The 2D Gaussian, built from the 1D weights so there is one set of
// numbers in this file and not two that can drift apart.
float h_gauss2d[kBlurDim * kBlurDim];
float h_gauss1d[kBlurDim];
for (int p = 0; p < kBlurDim; ++p) {
h_gauss1d[p] = static_cast<float>(kGauss1dInt[p]) /
static_cast<float>(kGauss1dSum);
for (int q = 0; q < kBlurDim; ++q) {
h_gauss2d[p * kBlurDim + q] =
static_cast<float>(kGauss1dInt[p] * kGauss1dInt[q]) /
static_cast<float>(kGauss2dSum);
}
}
CUDA_CHECK(cudaMemcpyToSymbol(cGauss2d, h_gauss2d, sizeof(h_gauss2d)));
CUDA_CHECK(cudaMemcpyToSymbol(cGauss1d, h_gauss1d, sizeof(h_gauss1d)));
CUDA_CHECK(cudaMemcpyToSymbol(cSobelX, kSobelXHost, sizeof(kSobelXHost)));
CUDA_CHECK(cudaMemcpyToSymbol(cSobelY, kSobelYHost, sizeof(kSobelYHost)));
// The image, written and read back through the day 7 Netpbm code. Reading
// it back rather than reusing the buffer already in memory is the point:
// 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 rows that are not the length it thinks they are.
std::vector<unsigned char> h_madeBytes;
makeEarthImage(&h_madeBytes, kBaseWidth, kBaseHeight);
if (!writeNetpbm("day20_input.pgm", h_madeBytes.data(), kBaseWidth,
kBaseHeight, 1)) {
std::fprintf(stderr, "could not write day20_input.pgm\n");
return EXIT_FAILURE;
}
std::vector<unsigned char> h_baseBytes;
size_t readWidth = 0;
size_t readHeight = 0;
size_t readChannels = 0;
if (!readNetpbm("day20_input.pgm", &h_baseBytes, &readWidth, &readHeight,
&readChannels)) {
std::fprintf(stderr, "could not read day20_input.pgm back\n");
return EXIT_FAILURE;
}
if (readWidth != kBaseWidth || readHeight != kBaseHeight ||
readChannels != 1 || h_baseBytes != h_madeBytes) {
std::fprintf(stderr,
"PGM round trip failed: wrote %zu x %zu x 1, read "
"%zu x %zu x %zu\n",
kBaseWidth, kBaseHeight, readWidth, readHeight,
readChannels);
return EXIT_FAILURE;
}
// One allocation at the largest size every case uses, so the small cases
// reuse it and nothing is allocated inside a loop or a timed region.
const size_t maxPixels = kPerfWidth * kPerfHeight;
const size_t maxFloatBytes = maxPixels * sizeof(float);
std::printf("timed case %zu x %zu, %.1f MiB per float buffer\n", kPerfWidth,
kPerfHeight,
static_cast<double>(maxFloatBytes) / (1024.0 * 1024.0));
float* d_in = nullptr;
float* d_blur = nullptr;
float* d_mag = nullptr;
unsigned char* d_bytes = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, maxFloatBytes));
CUDA_CHECK(cudaMalloc(&d_blur, maxFloatBytes));
CUDA_CHECK(cudaMalloc(&d_mag, maxFloatBytes));
CUDA_CHECK(cudaMalloc(&d_bytes, maxPixels));
const int status = runEverything(h_baseBytes, d_in, d_blur, d_mag, d_bytes);
// The one cleanup block. Every path above returns into it, including
// every failure, which is the pairing day 6 is about.
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_blur));
CUDA_CHECK(cudaFree(d_mag));
CUDA_CHECK(cudaFree(d_bytes));
return status;
}