code/day46-ptx-sass/ptx_sass.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 46: reading PTX and SASS.
//
// The interesting output of this program is not its stdout. It is the two
// disassemblies you take out of the binary once it is built, one from each
// half of the compiler:
//
// cuobjdump -ptx ./ptx_sass the virtual ISA nvcc handed to ptxas
// cuobjdump -sass ./ptx_sass the Turing machine code ptxas produced
//
// Three kernels, each chosen to put one specific thing in that second dump.
//
// decayRolled the fold's trip count is a kernel argument, so no stage
// of the compiler knows it. The loop survives into SASS,
// with a backward branch and a comparison per iteration.
// decayUnrolled the same fold over a compile-time constant, with
// #pragma unroll. The branch is gone and the body is a run
// of fused multiply-adds.
// decayStaged a two-pass window that keeps all kTaps values live at
// once, with __launch_bounds__ capping the register budget
// below what it wants. ptxas fits it by spilling, and a
// spill is an STL and an LDL instruction you can point at.
//
// decayRolled is handed taps = kTaps, so it must produce the same numbers as
// decayUnrolled, element for element: same operations, same order, different
// instruction count. The program checks that rather than promising it.
//
// decayStaged folds a running window rather than the raw input, so it
// computes a different value on purpose and has its own reference. Reading
// the window backwards is what keeps every tap live between the two loops,
// which is what creates the register pressure the launch bound then refuses
// to pay for.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o ptx_sass ptx_sass.cu
// Run: ./ptx_sass
//
#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)
// kElems is not a multiple of the block size, so the bounds check runs on
// every launch instead of never. The input buffer carries kTaps - 1 extra
// elements so the last window reads in[i + 63] without wrapping, which keeps
// a 64-bit modulo out of the inner loop. Day 11 paid for one of those and it
// moved the numbers it was there to measure.
constexpr size_t kElems = 1024ull * 1024ull + 611ull;
constexpr int kTaps = 64;
constexpr size_t kPaddedElems = kElems + kTaps - 1;
constexpr int kThreadsPerBlock = 256; // 8 warps
constexpr int kMinBlocksPerSm = 4;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr float kDecay = 0.5f;
constexpr float kRelTolerance = 1e-5f;
// Everything a compile-time constant already settles is a static_assert, not
// a run-time branch, because the compiler can answer it and the GPU should
// not have to. All four hold at any block size that is a whole number of
// warps, including 32, which is what the README invites you to try.
static_assert(kThreadsPerBlock % 32 == 0,
"block size must be a whole number of warps");
static_assert(kElems % kThreadsPerBlock != 0,
"pick an element count the block size does not divide, so the "
"bounds check runs on every launch");
static_assert(kTaps > 1, "a one-tap fold has no loop to unroll");
static_assert(kThreadsPerBlock * kMinBlocksPerSm <= 1024,
"sm_75 holds at most 1024 resident threads per SM (compute "
"capabilities appendix, Table 30), so a larger launch bound is "
"one the hardware cannot honour");
// out[i] folds a decaying window over in[i .. i + taps - 1], newest tap
// first. `taps` is a kernel argument, so neither the front end nor ptxas can
// see its value and neither can unroll the loop away. What survives into SASS
// is a loop: a comparison, a body, and a branch back to the top.
//
// One thread owns one output element. On every tap a warp's 32 lanes read 32
// consecutive floats, so each tap is four 32-byte sectors: the coalesced
// pattern from day 11.
//
// Launch assumption: gridDim.x * blockDim.x >= n, `in` holds at least
// n + kTaps - 1 elements, and the caller passes taps = kTaps so this kernel
// and decayUnrolled compute the same thing.
// snippet: rolled
__global__ void decayRolled(const float* __restrict__ in,
float* __restrict__ out, size_t n, int taps) {
const size_t i = blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
if (i < n) {
float acc = 0.0f;
for (int k = taps - 1; k >= 0; --k) {
acc = acc * kDecay + in[i + static_cast<size_t>(k)];
}
out[i] = acc;
}
}
// end snippet: rolled
// The same fold with the trip count nailed down at compile time and the
// pragma that asks for the loop to be replaced by its body, repeated. Same
// arithmetic in the same order as decayRolled, so the same answer, out of a
// different number of instructions.
//
// Memory and launch assumption: identical to decayRolled.
// snippet: unrolled
__global__ void decayUnrolled(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 acc = 0.0f;
#pragma unroll
for (int k = kTaps - 1; k >= 0; --k) {
acc = acc * kDecay + in[i + static_cast<size_t>(k)];
}
out[i] = acc;
}
}
// end snippet: unrolled
// A two-pass window with a register budget it cannot meet.
//
// The forward pass carries a running value and keeps every step, so all
// kTaps floats are live when the backward pass starts. __launch_bounds__(T, B)
// then promises ptxas that T threads per block and B blocks per SM will be
// resident, which is a register ceiling of the file divided by T times B.
// ptxas has to fit the kernel under that ceiling, and the values that no
// longer fit are written out and read back: spill stores and spill loads,
// STL and LDL in the disassembly.
//
// Memory: the reads of `in` are unchanged and still coalesced. The spill
// traffic is per thread and off-chip, because local memory is device DRAM
// with a per-thread address. Day 17 measured what that costs.
//
// Launch assumption: identical to decayRolled, plus the bound. A launch with
// more than kThreadsPerBlock threads per block now fails outright with
// "too many resources requested for launch" rather than running slowly.
// snippet: staged
__global__ __launch_bounds__(kThreadsPerBlock, kMinBlocksPerSm) void
decayStaged(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 window[kTaps];
float carry = 0.0f;
#pragma unroll
for (int k = 0; k < kTaps; ++k) {
carry = carry * kDecay + in[i + static_cast<size_t>(k)];
window[k] = carry;
}
float acc = 0.0f;
#pragma unroll
for (int k = kTaps - 1; k >= 0; --k) {
acc = acc * kDecay + window[k];
}
out[i] = acc;
}
}
// end snippet: staged
// CPU reference for decayRolled and decayUnrolled. Written for obvious
// correctness, not speed: plain loops, no OpenMP, no intrinsics. It
// accumulates in double where the kernels accumulate in float, because a
// reference exists to be right rather than to match bit for bit, and
// kRelTolerance covers the difference. It never allocates; the caller owns
// every buffer.
static void decayFoldCpu(const float* in, float* out, size_t n, int taps) {
for (size_t i = 0; i < n; ++i) {
double acc = 0.0;
for (int k = taps - 1; k >= 0; --k) {
acc = acc * kDecay + static_cast<double>(in[i + k]);
}
out[i] = static_cast<float>(acc);
}
}
// CPU reference for decayStaged, which folds the running window rather than
// the input and so is a different function of the same data.
static void decayStagedCpu(const float* in, float* out, size_t n, int taps) {
double window[kTaps];
for (size_t i = 0; i < n; ++i) {
double carry = 0.0;
for (int k = 0; k < taps; ++k) {
carry = carry * kDecay + static_cast<double>(in[i + k]);
window[k] = carry;
}
double acc = 0.0;
for (int k = taps - 1; k >= 0; --k) {
acc = acc * kDecay + window[k];
}
out[i] = static_cast<float>(acc);
}
}
// 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 point: "wrong at 512" 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.
// Copy it verbatim; it is the same helper day 9 handed you. The warm-up lives
// inside so it cannot go missing from the fourth kernel someone adds later.
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;
}
// One row of the table. The register and local-memory figures come from the
// driver's view of the compiled kernel, so the page cannot quote a number the
// binary does not have. localSizeBytes is a size_t and prints with %zu.
static void printRow(const char* name, const cudaFuncAttributes& attr,
float ms) {
std::printf("%-14s %5d %9zu %9.3f\n", name, attr.numRegs,
attr.localSizeBytes, ms);
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
// The ceiling the launch bound is measured against, read from the device
// rather than copied from a table. A card reporting no resident-thread
// limit would make the division meaningless, so this is a branch and not
// an assumption.
if (prop.maxThreadsPerMultiProcessor <= 0) {
std::fprintf(stderr, "device reports %d resident threads per SM\n",
prop.maxThreadsPerMultiProcessor);
return EXIT_FAILURE;
}
const int regsUnderBound =
prop.regsPerMultiprocessor / (kThreadsPerBlock * kMinBlocksPerSm);
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name, prop.major,
prop.minor);
std::printf("32-bit registers per SM: %d\n", prop.regsPerMultiprocessor);
std::printf("resident threads per SM: %d\n",
prop.maxThreadsPerMultiProcessor);
std::printf(
"__launch_bounds__(%d, %d) therefore caps decayStaged at %d "
"registers per thread\n\n",
kThreadsPerBlock, kMinBlocksPerSm, regsUnderBound);
std::printf("n = %zu elements, %d threads per block, %d taps\n\n", kElems,
kThreadsPerBlock, kTaps);
const size_t inBytes = kPaddedElems * sizeof(float);
const size_t outBytes = kElems * sizeof(float);
const int blocks =
static_cast<int>((kElems + kThreadsPerBlock - 1) / kThreadsPerBlock);
// Every host pointer carries h_ and every device pointer d_, because the
// prefix is the only thing standing between you and passing one where the
// other belongs.
std::vector<float> h_in(kPaddedElems);
std::vector<float> h_rolled(kElems);
std::vector<float> h_unrolled(kElems);
std::vector<float> h_staged(kElems);
std::vector<float> h_want(kElems);
for (size_t i = 0; i < kPaddedElems; ++i) {
h_in[i] = static_cast<float>(i % 17) * 0.25f;
}
float* d_in = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, inBytes));
CUDA_CHECK(cudaMalloc(&d_out, outBytes));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), inBytes, cudaMemcpyHostToDevice));
// One correctness launch per kernel, each followed by both checks in this
// order: cudaGetLastError() reports a launch the driver refused, and
// cudaDeviceSynchronize() reports what the kernel did.
decayRolled<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems, kTaps);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_rolled.data(), d_out, outBytes, cudaMemcpyDeviceToHost));
decayUnrolled<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_unrolled.data(), d_out, outBytes, cudaMemcpyDeviceToHost));
decayStaged<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(
cudaMemcpy(h_staged.data(), d_out, outBytes, cudaMemcpyDeviceToHost));
// What the driver thinks each compiled kernel needs. These are the ptxas
// report read from the other end, and the two have to agree.
cudaFuncAttributes attrRolled;
cudaFuncAttributes attrUnrolled;
cudaFuncAttributes attrStaged;
CUDA_CHECK(cudaFuncGetAttributes(&attrRolled, decayRolled));
CUDA_CHECK(cudaFuncGetAttributes(&attrUnrolled, decayUnrolled));
CUDA_CHECK(cudaFuncGetAttributes(&attrStaged, decayStaged));
const float msRolled = timeKernel([&] {
decayRolled<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems, kTaps);
});
const float msUnrolled = timeKernel([&] {
decayUnrolled<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
});
const float msStaged = timeKernel([&] {
decayStaged<<<blocks, kThreadsPerBlock>>>(d_in, d_out, kElems);
});
// Freed before anything is reported, so the failure paths below free too.
// Every cudaMalloc has a matching cudaFree before every return, including
// the return you take when the answer is wrong.
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
decayFoldCpu(h_in.data(), h_want.data(), kElems, kTaps);
size_t bad =
firstMismatch(h_rolled.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "decayRolled wrong at %zu: got %.9g, want %.9g\n",
bad, h_rolled[bad], h_want[bad]);
return EXIT_FAILURE;
}
bad =
firstMismatch(h_unrolled.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr,
"decayUnrolled wrong at %zu: got %.9g, want %.9g\n", bad,
h_unrolled[bad], h_want[bad]);
return EXIT_FAILURE;
}
decayStagedCpu(h_in.data(), h_want.data(), kElems, kTaps);
bad = firstMismatch(h_staged.data(), h_want.data(), kElems, kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr, "decayStaged wrong at %zu: got %.9g, want %.9g\n",
bad, h_staged[bad], h_want[bad]);
return EXIT_FAILURE;
}
// The claim the unrolling half of the page rests on: the two kernels are
// one program compiled two ways, not two programs. A real branch rather
// than an assert, because CI builds Release, Release defines NDEBUG, and
// NDEBUG deletes assert() from exactly the build that matters.
bad = firstMismatch(h_rolled.data(), h_unrolled.data(), kElems,
kRelTolerance);
if (bad != kElems) {
std::fprintf(stderr,
"decayRolled and decayUnrolled disagree at %zu (%.9g "
"against %.9g); they are meant to differ only in how "
"many instructions the fold becomes\n",
bad, h_rolled[bad], h_unrolled[bad]);
return EXIT_FAILURE;
}
// Reported, not gated. Both kernels issue the same operations in the same
// order, so every element should also be equal bit for bit. A handful of
// last-bit differences would mean ptxas contracted a multiply and an add
// in one of them and not the other, which is worth knowing and is not a
// failure.
size_t exactDiffs = 0;
for (size_t i = 0; i < kElems; ++i) {
if (h_rolled[i] != h_unrolled[i]) {
++exactDiffs;
}
}
std::printf("%-14s %5s %9s %9s\n", "kernel", "regs", "local B", "ms");
std::printf("%-14s %5s %9s %9s\n", "--------------", "-----", "---------",
"---------");
printRow("decayRolled", attrRolled, msRolled);
printRow("decayUnrolled", attrUnrolled, msUnrolled);
printRow("decayStaged", attrStaged, msStaged);
std::printf(
"\nall three kernels match their CPU reference; decayRolled and "
"decayUnrolled\ndiffer at %zu of %zu elements, so the two are one "
"program compiled twice\n",
exactDiffs, kElems);
return EXIT_SUCCESS;
}