code/day12-transpose/transpose.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 12: why transposing a matrix is slow.
//
// A transpose reads a row and writes a column. Consecutive threads take
// consecutive columns, so on one of the two sides the 32 lanes of a warp are
// `width` elements apart. At 4096 columns that is 16 KiB per lane, which is
// one 32-byte sector per lane instead of four sectors for the whole warp.
//
// Three kernels, one launch configuration, one byte count:
//
// copyRows coalesced read, coalesced write. The ceiling.
// transposeNaive coalesced read, strided write.
// transposeStridedRead strided read, coalesced write.
//
// The third exists because swapping which side pays is the first thing
// everyone tries, and it is not a fix. Nothing here is a fix. Day 13 fixes
// it with shared memory, which is also what NVIDIA's own guide says to do:
// "To alleviate these uncoalesced writes, the use of shared memory can be
// employed, which will be described in the next section." CUDA Programming
// Guide section 2.3.4.1, "Coalesced Global Memory Access", checked
// 2026-08-30. The full URL is in this directory's README and in the lesson.
//
// The first draft had no copy row and compared against day 11's coalesced
// copy, 232.9 GB/s on this project's verification node, which came from a 1D
// kernel at 256 threads a block over a 256 MiB buffer on another day. That
// would have folded the block shape, the grid shape, the buffer size and the
// card's clock state into a figure the page then calls "the transpose
// penalty". The ceiling is measured here, in this launch, or it is not one.
//
// Build: nvcc -std=c++17 -O3 -arch=sm_75 -o transpose transpose.cu
// Run: ./transpose
//
// 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)
// 4096 x 4096 floats, 64 MiB per buffer, 128 MiB moved per launch. Square, so
// the input and the output have the same shape and one buffer pair serves all
// three kernels.
//
// A power-of-two width is the case people actually write and the case that
// hurts most: the write stride is exactly 16 KiB, so consecutive lanes land in
// different sectors, different 128-byte lines and different DRAM pages at
// once. Day 15 comes back to what powers of two do to an address pattern.
constexpr size_t kWidth = 4096;
constexpr size_t kElems = kWidth * kWidth;
// A block is 32 by 32, which is 1024 threads, not the course default of 256.
// The tile shape forces it: 32 in x is exactly one warp per row of the block,
// which is what makes one side of every kernel here perfectly coalesced, and
// 32 in y makes the block square so one block covers a square patch of both
// matrices. 1024 is the maximum block on every GPU this course targets and it
// is 32 whole warps.
constexpr int kBlockDimX = 32;
constexpr int kBlockDimY = 32;
constexpr int kWarmupRuns = 3;
constexpr int kTimedRuns = 10;
constexpr int kNumRows = 3;
// Blocks along one side of the grid, by integer ceiling division.
//
// constexpr because the static_assert below calls it. A fact that follows only
// from the constants in this file is checked by the compiler, never at run
// time, and never with assert(): CI builds Release, Release defines NDEBUG,
// and assert() under NDEBUG expands to nothing.
constexpr int blocksPerSide(size_t width, int blockDim) {
return static_cast<int>((width + blockDim - 1) / blockDim);
}
static_assert(kBlockDimX == 32,
"one block row must be exactly one warp, or neither side of a "
"transpose is coalesced and the measurement means nothing");
static_assert(kBlockDimX * kBlockDimY % 32 == 0,
"block size must be a whole number of warps");
static_assert(kBlockDimX * kBlockDimY <= 1024,
"1024 threads per block is the maximum on every GPU this course "
"targets");
static_assert(static_cast<size_t>(blocksPerSide(kWidth, kBlockDimX)) *
kBlockDimX ==
kWidth,
"the grid must tile the matrix exactly on x, which is why no "
"kernel below carries a bounds check");
static_assert(static_cast<size_t>(blocksPerSide(kWidth, kBlockDimY)) *
kBlockDimY ==
kWidth,
"the grid must tile the matrix exactly on y as well");
static_assert(kElems <= (1ull << 24),
"every element index must be exactly representable as a float, "
"because the correctness check below compares exactly");
// out[row][col] = in[row][col]. One thread owns one element.
//
// Memory: threadIdx.x is the fastest-varying index, so the 32 lanes of a warp
// hold 32 consecutive values of `col` on one row. Both sides cover 128
// contiguous bytes, which is four 32-byte sectors read and four written.
//
// This is the ceiling row. It moves the same bytes as the two transposes
// below, from the same launch configuration, with the same index arithmetic.
// The only thing that differs is which subscript is swapped.
//
// Launch assumption: grid.x * blockDim.x == width and grid.y * blockDim.y ==
// width, checked by the static_asserts above, so there is no bounds check.
// snippet: copy-kernel
__global__ void copyRows(const float* __restrict__ in, float* __restrict__ out,
size_t width) {
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;
out[row * width + col] = in[row * width + col];
}
// end snippet: copy-kernel
// out[col][row] = in[row][col]. One thread owns one element.
//
// Memory: the read is the same coalesced 128 bytes as copyRows. The write is
// not. The 32 lanes hold consecutive `col`, so their write addresses are
// `width` elements apart, which is 16 KiB here. Every lane lands in its own
// 32-byte sector, so the warp asks for 32 sectors to deliver 128 useful
// bytes. NVIDIA's guide names the same thing in its own transpose example:
// "the writing of c is not coalesced, because consecutive values of
// threadIdx.x ... are writing elements to c that are ld (leading dimension)
// elements apart from each other."
//
// Launch assumption: as copyRows.
// snippet: naive-kernel
__global__ void transposeNaive(const float* __restrict__ in,
float* __restrict__ out, size_t width) {
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;
out[col * width + row] = in[row * width + col];
}
// end snippet: naive-kernel
// The same transpose with the two subscripts swapped, so the strided side is
// the read and the write is coalesced. Same output as transposeNaive, checked
// against the same expected values.
//
// Memory: 32 sectors read, four written. The sector count per warp is
// unchanged at 36. This kernel is here to show that moving the strided side
// is not a fix, so do not read it as one.
//
// Launch assumption: as copyRows.
// snippet: swapped-kernel
__global__ void transposeStridedRead(const float* __restrict__ in,
float* __restrict__ out, size_t width) {
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;
out[row * width + col] = in[col * width + row];
}
// end snippet: swapped-kernel
// Returns the first index where `got` is not a copy of the input, or n if it
// is a copy everywhere.
//
// The comparison is exact, and that is deliberate. These kernels move floats,
// they do not compute with them, so there is no rounding to allow for and a
// tolerance would hide a real index bug. main() fills every element with its
// own linear index, which the last static_assert above proves is exactly
// representable in float here.
static size_t firstBadCopy(const float* got, size_t n) {
for (size_t i = 0; i < n; ++i) {
if (got[i] != static_cast<float>(i)) {
return i;
}
}
return n;
}
// Returns the first index where `got` is not the transpose of the input, or
// width * width if the whole matrix is right. Same exactness argument.
static size_t firstBadTranspose(const float* got, size_t width) {
for (size_t row = 0; row < width; ++row) {
for (size_t col = 0; col < width; ++col) {
const float want = static_cast<float>(col * width + row);
if (got[row * width + col] != want) {
return row * width + col;
}
}
}
return width * width;
}
// Times a launch with CUDA events and returns the mean milliseconds per run.
// A host-side clock around a launch measures the launch, not the kernel,
// because launches are asynchronous. Day 9 takes that apart.
//
// This is the one template and the one lambda allowed in module 1 to 3 code.
// Copy it verbatim.
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 kernel here reads kElems floats and writes kElems floats. That
// equality is the whole point of the comparison, so the figure is computed
// once from one constant rather than per kernel, where the three could drift
// apart. Copies between host and device are not included and never enter a
// timed region.
static double bandwidthGBs(float ms) {
const double bytes = 2.0 * static_cast<double>(kElems) * sizeof(float);
return bytes / (static_cast<double>(ms) * 1.0e-3) / 1.0e9;
}
// Runs each kernel once and checks what it produced. Returns EXIT_SUCCESS only
// when all three are right.
//
// Nothing here is timed. The copies back would sit inside the measurement, and
// the first launch of each kernel still pays its own module load.
//
// d_out is filled with 0xFF bytes before every launch, which is a bit pattern
// that is not a number. An element the kernel never writes therefore fails the
// comparison instead of passing on whatever the previous kernel left there.
static int checkAll(const float* d_in, float* d_out, float* h_out, dim3 grid,
dim3 block, size_t bytes) {
int status = EXIT_SUCCESS;
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
copyRows<<<grid, block>>>(d_in, d_out, kWidth);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));
size_t bad = firstBadCopy(h_out, kElems);
if (bad != kElems) {
std::fprintf(stderr, "copyRows wrong at %zu: got %.9g, want %.9g\n",
bad, static_cast<double>(h_out[bad]),
static_cast<double>(bad));
status = EXIT_FAILURE;
}
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
transposeNaive<<<grid, block>>>(d_in, d_out, kWidth);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstBadTranspose(h_out, kWidth);
if (bad != kElems) {
std::fprintf(stderr, "transposeNaive wrong at %zu: got %.9g\n", bad,
static_cast<double>(h_out[bad]));
status = EXIT_FAILURE;
}
CUDA_CHECK(cudaMemset(d_out, 0xFF, bytes));
transposeStridedRead<<<grid, block>>>(d_in, d_out, kWidth);
CUDA_CHECK(cudaGetLastError());
CUDA_CHECK(cudaDeviceSynchronize());
CUDA_CHECK(cudaMemcpy(h_out, d_out, bytes, cudaMemcpyDeviceToHost));
bad = firstBadTranspose(h_out, kWidth);
if (bad != kElems) {
std::fprintf(stderr, "transposeStridedRead wrong at %zu: got %.9g\n",
bad, static_cast<double>(h_out[bad]));
status = EXIT_FAILURE;
}
return status;
}
int main() {
const int device = 0;
CUDA_CHECK(cudaSetDevice(device));
cudaDeviceProp prop;
CUDA_CHECK(cudaGetDeviceProperties(&prop, device));
const size_t bytes = kElems * sizeof(float);
const dim3 block(kBlockDimX, kBlockDimY);
const dim3 grid(blocksPerSide(kWidth, kBlockDimX),
blocksPerSide(kWidth, kBlockDimY));
std::vector<float> h_in(kElems);
std::vector<float> h_out(kElems);
for (size_t i = 0; i < kElems; ++i) {
h_in[i] = static_cast<float>(i);
}
float* d_in = nullptr;
float* d_out = nullptr;
CUDA_CHECK(cudaMalloc(&d_in, bytes));
CUDA_CHECK(cudaMalloc(&d_out, bytes));
CUDA_CHECK(cudaMemcpy(d_in, h_in.data(), bytes, cudaMemcpyHostToDevice));
int status = checkAll(d_in, d_out, h_out.data(), grid, block, bytes);
float copyMs = 0.0f;
float naiveMs = 0.0f;
float swappedMs = 0.0f;
if (status == EXIT_SUCCESS) {
copyMs =
timeKernel([&] { copyRows<<<grid, block>>>(d_in, d_out, kWidth); });
naiveMs = timeKernel(
[&] { transposeNaive<<<grid, block>>>(d_in, d_out, kWidth); });
swappedMs = timeKernel([&] {
transposeStridedRead<<<grid, block>>>(d_in, d_out, kWidth);
});
// A real branch returning EXIT_FAILURE, not an assert(), because CI
// builds Release and NDEBUG deletes assert(). This is arithmetic and
// not a claim about your card: a zero or negative elapsed time means
// the event pair never separated, and every GB/s figure derived from
// it would be infinite or nonsense.
if (copyMs <= 0.0f || naiveMs <= 0.0f || swappedMs <= 0.0f) {
std::fprintf(stderr,
"timing broken: elapsed times were %.6f, %.6f and "
"%.6f ms; none of them can be zero\n",
static_cast<double>(copyMs),
static_cast<double>(naiveMs),
static_cast<double>(swappedMs));
status = EXIT_FAILURE;
}
}
if (status == EXIT_SUCCESS) {
std::printf("GPU: %s (compute capability %d.%d)\n", prop.name,
prop.major, prop.minor);
std::printf("Matrix: %zu x %zu floats, %.0f MiB per buffer\n", kWidth,
kWidth, static_cast<double>(bytes) / (1024.0 * 1024.0));
std::printf(
"Block: %d x %d threads. Mean of %d runs after %d "
"warm-ups.\n\n",
kBlockDimX, kBlockDimY, kTimedRuns, kWarmupRuns);
std::printf("%-34s %10s %10s %9s\n", "kernel", "time (ms)", "GB/s",
"of copy");
std::printf("%-34s %10s %10s %9s\n",
"---------------------------------", "---------",
"---------", "--------");
std::printf("%-34s %10.3f %10.1f %9.2f\n", "copy, both sides coalesced",
copyMs, bandwidthGBs(copyMs), 1.0);
std::printf("%-34s %10.3f %10.1f %9.2f\n", "transpose, strided write",
naiveMs, bandwidthGBs(naiveMs),
static_cast<double>(copyMs / naiveMs));
std::printf("%-34s %10.3f %10.1f %9.2f\n", "transpose, strided read",
swappedMs, bandwidthGBs(swappedMs),
static_cast<double>(copyMs / swappedMs));
std::printf(
"\nAll %d rows moved exactly %.0f MiB, from the same grid and the "
"same\nblock, with the same index arithmetic. Only which "
"subscript is\nswapped differs, so the gap is the address pattern "
"and nothing else.\n",
kNumRows, 2.0 * static_cast<double>(bytes) / (1024.0 * 1024.0));
}
CUDA_CHECK(cudaFree(d_in));
CUDA_CHECK(cudaFree(d_out));
return status;
}