Appendix K: Matrix Traversal, From First Principles — Row-Major, Column-Major, Transpose, Spiral, Diagonal, and the Grids That Read Them¶
"Chapter 6 and 7 gave a matrix's storage a name — row-major, the last axis contiguous — and Chapter 13 built the two operations, multiplication and transpose, that actually move through that storage. Chapter 18 gave a CUDA kernel its own grid of threads and the ceiling-division-plus-bounds-check pattern that keeps every one of them from stepping off the edge. This appendix is the place all four of those ideas meet on one object at once: a 2D matrix is nothing more than a flat buffer and an addressing formula, a CUDA grid is nothing more than that same formula computed once per thread instead of once per loop iteration, and a boundary check is nothing more than making sure the formula was only ever asked to answer for a position that actually exists. That same lens turns out to reach two traversal orders that look, at first, nothing like a flat formula at all — spiral order, whose rank at a position depends on a shrinking ring boundary, and diagonal order, whose whole point is that some positions depend on others. Every claim below — which traversal order is fast, which axis a real bank-conflict fix survives, which boundary check is only half-written, which layout a real cuBLAS call expects, which of two intricate-looking traversals still fits one kernel launch and which one genuinely cannot — is checked by a genuine compile and, everywhere a real answer doesn't require a GPU to give one honestly, a genuine run."
What you will understand by the end of this appendix¶
- That row-major STORAGE and the order a loop nest TRAVERSES that storage are two independent decisions — proven with a real, deterministic cache simulator (not a wall-clock timing, which this book's own build environment already notes will jitter run to run) that shows a matching traversal order hitting a 93.75% hit rate and a crossing order hitting exactly zero, on the identical data, in the identical memory.
- A complete, non-square, boundary-checked shared-memory tiled transpose kernel — genuinely compiled for a real architecture and, since launching it needs a device this sandbox does not have, checked instead by a host-side simulator that walks the same grid/block/tile loop nest a real launch would, exactly the way this book's own Appendix E.3 checks a memory-traffic count without a GPU. Appendix C.3's bank-conflict fix (
[32][32]padded to[32][33]) is confirmed, by the identical arithmetic, to survive completely unchanged once boundary checks are added. - Why CUDA's
blockIdx.y * blockDim.y + threadIdx.y/blockIdx.x * blockDim.x + threadIdx.xconvention puts columns onxand rows onyspecifically for coalescing, exhaustively verified (every position in several differently-shaped matrices touched exactly once, not zero, not twice) rather than merely asserted — and why that assignment flips completely for column-major data such as cuBLAS expects. - A genuine, scheduling-dependent data race produced by a boundary check that guards only one of a 2D kernel's two axes — shown to be silently harmless for element-wise addition specifically (a real, worked-out reason why, involving flat-index aliasing) and genuinely corrupting for a position-dependent fill kernel, with the corruption's very presence shown to depend on an execution order CUDA never promises.
- A working row-major/column-major converter for handing this book's own default
Tensorlayout to a column-major API, a genuine round-trip check that the conversion loses nothing, and Chapter 13.3's naive-reshape trap shown to be exactly what happens when a column-major buffer is read as row-major without converting it at all. - Spiral order's genuine O(1)-per-position closed form — derived from nothing but distance to the nearest edge, cross-checked exhaustively against a sequential simulation across ten shapes, and then handed to a real CUDA kernel where every thread computes its own rank with no dependency on any other thread, despite spiral order's reputation as an inherently sequential traversal.
- Diagonal order's opposite lesson: every position on one anti-diagonal is independent of every other position on that same diagonal, but a genuine wavefront dependency exists ACROSS diagonals — demonstrated with a real lattice-path dynamic-programming kernel that provably cannot be launched once for the whole grid the way this book's ordinary elementwise kernels are, and is instead genuinely simulated as the exact sequence of per-diagonal launches it would actually need.
What you need to know first¶
- Chapter 6.2 and Chapter 7.2 in full —
TensorStrides::row_majorandcol_major, and the offset formula both layouts share. This appendix's ownRowMajorMatrix/ColMajorMatrixstructs specialize those exact formulas down to two dimensions, the same way Chapter 13'sMatrix::getalready did for row-major alone. - Chapter 13.1 and 13.2 — the shape-check gap's sentinel technique (reused directly in Section K.4) and the SIMD-scatter transpose (Section K.2's GPU tiled transpose is the same operation, moved onto a device and given real shared memory).
- Chapter 18.1 and 18.3 in full — ceiling-division launch configuration, the two-part "compute an index, then guard it" template every kernel in this book follows, and the tiled convolution kernel Section K.2's transpose kernel structurally mirrors.
- Appendix C.3 (shared-memory bank conflicts) and Appendix E.3's own
[COMMON TRAP], which named the exact gap this appendix's Section K.2 closes: "a real kernel handling a matrix size that doesn't divide evenly needs the same bounds-checking and padding logic." - Appendix H.7's glossary entry for a cache line (
commonly 64 bytes, or 8 doubles) — Section K.1's cache simulator uses this book's own stated figure, not a new assumption.
K.1 Row-Major versus Column-Major: Storage Is Not Traversal [FOUNDATIONAL]¶
Intuition¶
A row-major matrix's storage answers exactly one question: which flat slot does logical position (row, col) live in. It says nothing about the order any particular piece of code is required to visit those slots in. A loop nest that visits row outer and col inner walks the buffer front to back, one cache line at a time, each line fully consumed before the next is touched. A loop nest that visits col outer and row inner walks the identical buffer in enormous strides, abandoning most of every line it fetches before using more than one value from it. Both loops compute the exact same sum over the exact same data — this section makes both facts, sameness of the answer and enormous difference in the cost, genuinely observable at once.
Background¶
The file below reuses Chapter 13's exact Matrix struct (float* data; int rows, cols;, with get/set computing row * cols + col) as its row-major baseline, alongside a ColMajorMatrix specializing Chapter 7.2's col_major offset formula the same way. Measuring the cost honestly means not measuring wall-clock time: Getting Started's own account of this book's build environment already notes that "any chapter's runtime numbers, timing comparisons... will jitter somewhat from run to run and machine to machine" (Chapter 19's own GFLOPS figures). A cache-line-touch count is different — it is pure index arithmetic, reproduces bit-for-bit on every rerun, on any machine, and is exactly the same kind of genuinely-counted-not-timed evidence Appendix E.3 already used for a memory-traffic claim. This file implements a small, fully deterministic LRU (least-recently-used) cache simulator, capacity given in lines rather than bytes, and feeds it the exact sequence of cache lines two traversal orders touch.
File: 01_row_major_col_major_traversal.cpp¶
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <list>
#include <unordered_map>
// Appendix K.1 -- row-major storage (Chapter 3's convention, Chapter 13's
// Matrix::get) fixes WHERE every element lives. It says nothing at all
// about the order a loop nest is allowed to VISIT those elements in. This
// file keeps storage fixed and varies only the traversal order, then
// measures the consequence with a real, deterministic cache simulator --
// not a wall-clock timing, which Getting Started's own account of this
// book's build environment already notes will "jitter somewhat from run
// to run and machine to machine" (Chapter 19's own GFLOPS figures). A
// cache-line-touch count, computed from pure index arithmetic, reproduces
// bit-for-bit on every rerun, on any machine, exactly like Appendix E.3's
// genuinely counted (not timed) global-memory reads.
struct Matrix {
float* data;
int rows, cols;
Matrix(int rows_, int cols_) : rows(rows_), cols(cols_) {
data = (float*)malloc((size_t)rows * cols * sizeof(float));
}
~Matrix() { free(data); }
float get(int row, int col) const { return data[row * cols + col]; } // row-major, Chapter 3
void set(int row, int col, float value) { data[row * cols + col] = value; }
};
// The same 2x3 buffer, but stored COLUMN-major instead -- Chapter 7.2's
// col_major offset formula (offset = col*rows + row) specialized to two
// dimensions, the same way Matrix::get specialized Chapter 6's general
// row-major formula down to two.
struct ColMajorMatrix {
float* data;
int rows, cols;
ColMajorMatrix(int rows_, int cols_) : rows(rows_), cols(cols_) {
data = (float*)malloc((size_t)rows * cols * sizeof(float));
}
~ColMajorMatrix() { free(data); }
float get(int row, int col) const { return data[col * rows + row]; } // column-major, Chapter 7.2
void set(int row, int col, float value) { data[col * rows + row] = value; }
};
// A CPU cache line is a fixed-size chunk (this book's own glossary, Appendix
// H.7, states "commonly 64 bytes, or 8 doubles" -- for 4-byte float this is
// 16 elements per line) fetched as one unit on a miss. A tiny, fully
// deterministic LRU (least-recently-used) simulator, capacity given in
// LINES rather than bytes, reproduces the real replacement decision a
// hardware cache makes without needing to touch actual memory or a clock.
struct LRUCacheSim {
int capacity_lines;
std::list<long long> lru_order; // front = most recently used
std::unordered_map<long long, std::list<long long>::iterator> present;
long long hits = 0, misses = 0;
explicit LRUCacheSim(int capacity_lines_) : capacity_lines(capacity_lines_) {}
void access(long long line_id) {
auto it = present.find(line_id);
if (it != present.end()) {
hits++;
lru_order.erase(it->second);
lru_order.push_front(line_id);
present[line_id] = lru_order.begin();
return;
}
misses++;
if ((int)lru_order.size() >= capacity_lines) {
long long victim = lru_order.back();
lru_order.pop_back();
present.erase(victim);
}
lru_order.push_front(line_id);
present[line_id] = lru_order.begin();
}
};
static const int FLOATS_PER_LINE = 16; // 64-byte line / 4-byte float, Appendix H.7's glossary
long long line_of(int row, int col, int cols) {
long long flat = (long long)row * cols + col;
return flat / FLOATS_PER_LINE;
}
// Walk a row-major-stored matrix outer-row/inner-col ("natural" order for
// row-major data) and feed every touched line, in visit order, to a fresh
// cache simulator of the given capacity.
void simulate_row_outer(int rows, int cols, int capacity_lines, long long& hits, long long& misses) {
LRUCacheSim sim(capacity_lines);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
sim.access(line_of(r, c, cols));
hits = sim.hits;
misses = sim.misses;
}
// Walk the SAME row-major-stored matrix outer-col/inner-row instead --
// storage never changes, only which index the loop nest holds fixed.
void simulate_col_outer(int rows, int cols, int capacity_lines, long long& hits, long long& misses) {
LRUCacheSim sim(capacity_lines);
for (int c = 0; c < cols; c++)
for (int r = 0; r < rows; r++)
sim.access(line_of(r, c, cols));
hits = sim.hits;
misses = sim.misses;
}
void run_case(const char* label, int rows, int cols, int capacity_lines) {
long long row_hits, row_misses, col_hits, col_misses;
simulate_row_outer(rows, cols, capacity_lines, row_hits, row_misses);
simulate_col_outer(rows, cols, capacity_lines, col_hits, col_misses);
long long total = (long long)rows * cols;
printf("%s: %dx%d matrix, %d-line cache (%d floats/line -> %d floats resident)\n",
label, rows, cols, capacity_lines, FLOATS_PER_LINE, capacity_lines * FLOATS_PER_LINE);
printf(" row-outer (matches row-major storage): %lld hits, %lld misses / %lld accesses (%.2f%% hit rate)\n",
row_hits, row_misses, total, 100.0 * row_hits / total);
printf(" col-outer (crosses row-major storage): %lld hits, %lld misses / %lld accesses (%.2f%% hit rate)\n",
col_hits, col_misses, total, 100.0 * col_hits / total);
}
int main() {
printf("=== Section K.1: row-major vs column-major TRAVERSAL of row-major-stored data ===\n\n");
run_case("Case 1", 8, 64, 4);
long long case1_row_hits, case1_row_misses, case1_col_hits, case1_col_misses;
simulate_row_outer(8, 64, 4, case1_row_hits, case1_row_misses);
simulate_col_outer(8, 64, 4, case1_col_hits, case1_col_misses);
printf("\n");
run_case("Case 2 (larger, same shape ratio)", 32, 128, 4);
printf("\n=== Correctness check: traversal order never changes the answer, only the cost ===\n");
Matrix m(8, 64);
for (int r = 0; r < m.rows; r++)
for (int c = 0; c < m.cols; c++)
m.set(r, c, (float)(r * 64 + c));
double sum_row_outer = 0.0;
for (int r = 0; r < m.rows; r++)
for (int c = 0; c < m.cols; c++)
sum_row_outer += m.get(r, c);
double sum_col_outer = 0.0;
for (int c = 0; c < m.cols; c++)
for (int r = 0; r < m.rows; r++)
sum_col_outer += m.get(r, c);
printf("sum via row-outer traversal: %.1f\n", sum_row_outer);
printf("sum via col-outer traversal: %.1f\n", sum_col_outer);
printf("identical sums, wildly different cache behavior: %s\n",
(sum_row_outer == sum_col_outer) ? "true" : "false");
printf("\n=== Column-major storage, SAME 8x64 shape: its own natural order vs crossing it ===\n");
ColMajorMatrix cm(8, 64);
for (int r = 0; r < cm.rows; r++)
for (int c = 0; c < cm.cols; c++)
cm.set(r, c, (float)(r * 64 + c));
// Reuse the identical line-touch machinery, but now the "natural" (fast)
// order for THIS storage is col-outer/row-inner, since col_major's
// offset formula makes consecutive ROWS at a fixed column contiguous.
LRUCacheSim cm_natural(4), cm_crossed(4);
for (int c = 0; c < cm.cols; c++)
for (int r = 0; r < cm.rows; r++)
cm_natural.access(((long long)c * cm.rows + r) / FLOATS_PER_LINE);
for (int r = 0; r < cm.rows; r++)
for (int c = 0; c < cm.cols; c++)
cm_crossed.access(((long long)c * cm.rows + r) / FLOATS_PER_LINE);
long long cm_total = (long long)cm.rows * cm.cols;
printf("col-major storage, col-outer traversal (matches storage): %lld hits, %lld misses / %lld (%.2f%% hit rate)\n",
cm_natural.hits, cm_natural.misses, cm_total, 100.0 * cm_natural.hits / cm_total);
printf("col-major storage, row-outer traversal (crosses storage): %lld hits, %lld misses / %lld (%.2f%% hit rate)\n",
cm_crossed.hits, cm_crossed.misses, cm_total, 100.0 * cm_crossed.hits / cm_total);
printf("matching traversal reaches the identical 93.75%% hit rate row-major's own matching case did: %s\n",
(cm_natural.hits == case1_row_hits) ? "true" : "false");
printf("\n=== Why the CROSSING numbers above are NOT a mirror image of each other ===\n");
printf("row-major 8x64, crossing stride = cols = 64 floats = %d lines (>= cache capacity 4): every access misses\n", 64 / FLOATS_PER_LINE);
printf("col-major 8x64, crossing stride = rows = 8 floats = half of one %d-float line: partial reuse survives\n", FLOATS_PER_LINE);
printf("row-major crossing: %lld hits, %lld misses (stride a whole multiple of line size)\n", case1_col_hits, case1_col_misses);
printf("col-major crossing: %lld hits, %lld misses (stride smaller than one line)\n", cm_crossed.hits, cm_crossed.misses);
printf("conclusion: the SEVERITY of crossing storage depends on the stride relative to line size,\n");
printf("not just on whether storage was crossed at all -- both are still worse than matching, but not equally worse.\n");
printf("\n=== A genuinely SYMMETRIC (square) shape DOES mirror exactly ===\n");
const int SQ = 32;
long long sq_row_hits, sq_row_misses, sq_col_hits, sq_col_misses;
simulate_row_outer(SQ, SQ, 4, sq_row_hits, sq_row_misses); // row-major storage
simulate_col_outer(SQ, SQ, 4, sq_col_hits, sq_col_misses);
LRUCacheSim sq_cm_natural(4), sq_cm_crossed(4);
for (int c = 0; c < SQ; c++)
for (int r = 0; r < SQ; r++)
sq_cm_natural.access(((long long)c * SQ + r) / FLOATS_PER_LINE);
for (int r = 0; r < SQ; r++)
for (int c = 0; c < SQ; c++)
sq_cm_crossed.access(((long long)c * SQ + r) / FLOATS_PER_LINE);
printf("row-major %dx%d, matching traversal: %lld hits, %lld misses\n", SQ, SQ, sq_row_hits, sq_row_misses);
printf("col-major %dx%d, matching traversal: %lld hits, %lld misses\n", SQ, SQ, sq_cm_natural.hits, sq_cm_natural.misses);
printf("row-major %dx%d, crossing traversal: %lld hits, %lld misses\n", SQ, SQ, sq_col_hits, sq_col_misses);
printf("col-major %dx%d, crossing traversal: %lld hits, %lld misses\n", SQ, SQ, sq_cm_crossed.hits, sq_cm_crossed.misses);
bool square_mirrors = (sq_row_hits == sq_cm_natural.hits) && (sq_col_hits == sq_cm_crossed.hits);
printf("for a square matrix, row-major and column-major produce IDENTICAL hit/miss counts both ways: %s\n",
square_mirrors ? "true" : "false");
return 0;
}
Worked Example K.1.1 — An 8x64 matrix, a 4-line cache, and a 4x difference measured exactly¶
g++ -std=c++17 -Wall -Wextra -O2 01_row_major_col_major_traversal.cpp -o 01_row_major_col_major_traversal
./01_row_major_col_major_traversal
Genuinely compiled and run:
=== Section K.1: row-major vs column-major TRAVERSAL of row-major-stored data ===
Case 1: 8x64 matrix, 4-line cache (16 floats/line -> 64 floats resident)
row-outer (matches row-major storage): 480 hits, 32 misses / 512 accesses (93.75% hit rate)
col-outer (crosses row-major storage): 0 hits, 512 misses / 512 accesses (0.00% hit rate)
Case 2 (larger, same shape ratio): 32x128 matrix, 4-line cache (16 floats/line -> 64 floats resident)
row-outer (matches row-major storage): 3840 hits, 256 misses / 4096 accesses (93.75% hit rate)
col-outer (crosses row-major storage): 0 hits, 4096 misses / 4096 accesses (0.00% hit rate)
=== Correctness check: traversal order never changes the answer, only the cost ===
sum via row-outer traversal: 130816.0
sum via col-outer traversal: 130816.0
identical sums, wildly different cache behavior: true
=== Column-major storage, SAME 8x64 shape: its own natural order vs crossing it ===
col-major storage, col-outer traversal (matches storage): 480 hits, 32 misses / 512 (93.75% hit rate)
col-major storage, row-outer traversal (crosses storage): 256 hits, 256 misses / 512 (50.00% hit rate)
matching traversal reaches the identical 93.75% hit rate row-major's own matching case did: true
=== Why the CROSSING numbers above are NOT a mirror image of each other ===
row-major 8x64, crossing stride = cols = 64 floats = 4 lines (>= cache capacity 4): every access misses
col-major 8x64, crossing stride = rows = 8 floats = half of one 16-float line: partial reuse survives
row-major crossing: 0 hits, 512 misses (stride a whole multiple of line size)
col-major crossing: 256 hits, 256 misses (stride smaller than one line)
conclusion: the SEVERITY of crossing storage depends on the stride relative to line size,
not just on whether storage was crossed at all -- both are still worse than matching, but not equally worse.
=== A genuinely SYMMETRIC (square) shape DOES mirror exactly ===
row-major 32x32, matching traversal: 960 hits, 64 misses
col-major 32x32, matching traversal: 960 hits, 64 misses
row-major 32x32, crossing traversal: 0 hits, 1024 misses
col-major 32x32, crossing traversal: 0 hits, 1024 misses
for a square matrix, row-major and column-major produce IDENTICAL hit/miss counts both ways: true
An 8x64, row-major-stored matrix with a 4-line (64-float) cache: row-outer traversal touches each row's four lines once each (32 compulsory misses total, one per line, ever) and then reuses every one of the other 480 accesses as a hit, a 93.75% hit rate — the theoretical best this cache size can do, since every line must miss at least once. Col-outer traversal on the exact same buffer strides by 64 floats (exactly 4 lines) between consecutive accesses within its inner loop: with only 4 lines of capacity and 8 rows to visit per column, every single one of the 8 lines a column touches is distinct, so the cache — which can hold only 4 — never has a chance to reuse any of them before the next column starts the identical pattern over. Zero hits, 512 misses, out of 512 accesses. The correctness check confirms both traversals compute the identical sum over the identical data; only the cache behavior differs.
The column-major mirror section does not simply repeat the row-major result with the axes swapped: for the same 8x64 shape, column-major storage's OWN crossing traversal (row-outer, which crosses column-major's natural column-outer order) strides by only 8 floats — half of one cache line — so partial reuse survives, landing at a 50% hit rate rather than row-major's 0%. The severity of crossing a layout's natural order depends on how large the resulting stride is relative to one line, not merely on whether the order was crossed at all. A genuinely square shape removes that asymmetry entirely: at 32x32, row-major and column-major produce bit-for-bit identical hit and miss counts in both the matching and the crossing direction, confirmed by the same simulator run on both layouts.
[COMMON TRAP] "Row-major storage" and "fast traversal" are not the same claim. A loop nest that iterates for (col) for (row)) over a row-major matrix — an easy typo away from the correct for (row) for (col), and something Chapter 3's row-major convention alone does nothing to prevent — computes exactly the same answer this section's correctness check confirmed, at exactly the catastrophic cache behavior this section's simulator measured. The storage convention only ever describes where data lives; the loop nest that was actually written is what decides whether that convention pays off.
K.2 Transpose on the GPU: Bank Conflicts Meet Boundary Checks [FOUNDATIONAL]¶
Intuition¶
Appendix C.3 built the classic shared-memory transpose and its bank-conflict fix, but only for a square matrix whose width is already a clean multiple of 32 — shared_transpose_naive and shared_transpose_padded there never check a single bound. Real matrices are rarely 32k x 32k. This section builds the kernel Appendix E.3's own [COMMON TRAP] named as missing: the same tile, the same padding fix, generalized to any rows x cols shape and fully guarded against reading or writing past either matrix's true edge.
Background¶
The two-phase structure — cooperatively load a tile into shared memory, synchronize, then write the transposed tile back out at swapped block coordinates — is Chapter 18.3's tiled convolution kernel's own shape, applied to a different operation. What changes for a non-square matrix is that the OUTPUT tile's block coordinates are not simply "the same block, transposed in place" the way Appendix C.3's square case could get away with; they are blockIdx.x and blockIdx.y swapped, each checked against the output's own, independently-sized bound.
File: 02_matrix_transpose_tiled_boundary.cu¶
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <set>
#include <vector>
// Appendix K.2 -- Appendix C.3 built the classic shared-memory tiled
// transpose and its bank-conflict fix (a [32][32] tile padded to [32][33]),
// but only for a SQUARE matrix whose width divides evenly by 32 -- neither
// `shared_transpose_naive` nor `shared_transpose_padded` there ever checks
// a bound, and both would silently read or write past the end of a matrix
// that isn't exactly 32k x 32k. Appendix E.3's own [COMMON TRAP] named
// this gap directly ("a real kernel handling a matrix size that doesn't
// divide evenly needs the same bounds-checking and padding logic") without
// building it. This file builds it: a real, non-square-capable,
// boundary-checked tiled transpose, genuinely compiled with nvcc for a
// real target architecture and, since launching it needs a device this
// sandbox does not have, checked instead by a host-side simulator that
// walks the identical grid/block/tile loop nest a real launch would --
// exactly Appendix E.3's own verification technique, applied to a kernel
// instead of a memory-traffic count.
#define TILE 32
// The naive version: no padding, so a column read of the shared tile
// still produces Appendix C.3's exact 32-way bank conflict -- but now
// ALSO correctly bounded for a rows x cols matrix that need not be
// square or a multiple of TILE. Genuinely compiled for sm_80; never
// launched in this no-GPU environment (Getting Started's own account of
// this book's build environment).
__global__ void transpose_naive_boundary_checked(float* out, const float* in, int rows, int cols) {
__shared__ float tile[TILE][TILE];
int in_row = blockIdx.y * TILE + threadIdx.y;
int in_col = blockIdx.x * TILE + threadIdx.x;
if (in_row < rows && in_col < cols) {
tile[threadIdx.y][threadIdx.x] = in[in_row * cols + in_col];
}
__syncthreads();
// out is cols x rows (the transposed shape): swap which block axis
// supplies the row vs column half of the output coordinate, exactly
// the way Appendix C.3's square version swaps x and y in one line --
// generalized here into two full, independently-bounded coordinates.
int out_row = blockIdx.x * TILE + threadIdx.y;
int out_col = blockIdx.y * TILE + threadIdx.x;
if (out_row < cols && out_col < rows) {
out[out_row * rows + out_col] = tile[threadIdx.x][threadIdx.y]; // column read -- the conflict
}
}
// The fix: Appendix C.3's identical one-float row padding, plus the same
// boundary checks as the naive version above. Padding and bounds-checking
// are orthogonal concerns -- one is about shared-memory bank arithmetic,
// the other about which threads own real data -- and this kernel is the
// genuine, compiled proof that combining them needs no new idea, just
// both fixes applied at once.
__global__ void transpose_padded_boundary_checked(float* out, const float* in, int rows, int cols) {
__shared__ float tile[TILE][TILE + 1];
int in_row = blockIdx.y * TILE + threadIdx.y;
int in_col = blockIdx.x * TILE + threadIdx.x;
if (in_row < rows && in_col < cols) {
tile[threadIdx.y][threadIdx.x] = in[in_row * cols + in_col];
}
__syncthreads();
int out_row = blockIdx.x * TILE + threadIdx.y;
int out_col = blockIdx.y * TILE + threadIdx.x;
if (out_row < cols && out_col < rows) {
out[out_row * rows + out_col] = tile[threadIdx.x][threadIdx.y];
}
}
// ---- Everything below is host C++, genuinely compiled AND genuinely run ----
// Appendix C.3's exact bank-index arithmetic, reused verbatim on this
// file's own row strides to confirm the fix still holds unchanged.
int compute_distinct_banks(int row_stride_words) {
std::set<int> banks;
for (int t = 0; t < 32; t++) {
banks.insert((t * row_stride_words) % 32);
}
return (int)banks.size();
}
// A host-side simulator of the exact grid/block/thread/two-phase loop nest
// both kernels above run on a device -- Appendix E.3's own verification
// technique (simulate the launch, don't just trust the formula), applied
// here to a transpose instead of a read count. tile_width is a RUNTIME
// parameter (like E.3's own TILE_WIDTH), independent of the kernels'
// compile-time TILE, specifically so a small, hand-checkable tile_width
// can exercise the boundary path without touching the real kernels' macro.
void simulate_tiled_transpose(const std::vector<float>& in, int rows, int cols,
int tile_width, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f); // out is cols x rows; same element count
int blocks_x = (cols + tile_width - 1) / tile_width; // ceiling division, Chapter 18.1
int blocks_y = (rows + tile_width - 1) / tile_width;
std::vector<float> tile((size_t)tile_width * tile_width);
for (int by = 0; by < blocks_y; by++) {
for (int bx = 0; bx < blocks_x; bx++) {
// Phase 1: cooperative load into the simulated shared tile.
for (int ty = 0; ty < tile_width; ty++) {
for (int tx = 0; tx < tile_width; tx++) {
int in_row = by * tile_width + ty;
int in_col = bx * tile_width + tx;
if (in_row < rows && in_col < cols) {
tile[ty * tile_width + tx] = in[(size_t)in_row * cols + in_col];
}
}
}
// __syncthreads(): on real hardware this barrier is why phase 2
// may not begin until every thread in the block has finished
// phase 1. A sequential host simulation is already ordered, so
// the barrier is implied by nothing more than "this loop already
// finished" -- the same reasoning Appendix E.3 relies on when it
// simulates a launch without a real device to synchronize on.
// Phase 2: write the transposed tile out, boundary-checked
// against the OUTPUT shape (cols x rows), swapped block axes.
for (int ty = 0; ty < tile_width; ty++) {
for (int tx = 0; tx < tile_width; tx++) {
int out_row = bx * tile_width + ty;
int out_col = by * tile_width + tx;
if (out_row < cols && out_col < rows) {
out[(size_t)out_row * rows + out_col] = tile[tx * tile_width + ty];
}
}
}
}
}
}
void scalar_transpose_reference(const std::vector<float>& in, int rows, int cols, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
out[(size_t)c * rows + r] = in[(size_t)r * cols + c];
}
bool run_case(const char* label, int rows, int cols, int tile_width) {
std::vector<float> in((size_t)rows * cols);
for (size_t i = 0; i < in.size(); i++) in[i] = (float)(i + 1);
std::vector<float> tiled_out, reference_out;
simulate_tiled_transpose(in, rows, cols, tile_width, tiled_out);
scalar_transpose_reference(in, rows, cols, reference_out);
bool match = (tiled_out == reference_out);
int blocks_x = (cols + tile_width - 1) / tile_width;
int blocks_y = (rows + tile_width - 1) / tile_width;
bool boundary_case = (rows % tile_width != 0) || (cols % tile_width != 0);
printf("%s: %dx%d matrix, tile_width=%d, grid=%dx%d blocks (boundary tiles: %s)\n",
label, rows, cols, tile_width, blocks_x, blocks_y, boundary_case ? "yes" : "no");
printf(" tiled+boundary-checked simulation matches scalar reference, element for element: %s\n",
match ? "true" : "false");
return match;
}
int main() {
printf("=== Section K.2: tiled transpose with boundary checks, genuinely simulated and verified ===\n\n");
bool all_ok = true;
all_ok &= run_case("Case 1 (rows and cols both divide tile_width evenly)", 8, 8, 4);
all_ok &= run_case("Case 2 (neither dimension divides tile_width evenly)", 7, 5, 4);
all_ok &= run_case("Case 3 (larger, still non-divisible, non-square)", 10, 13, 4);
all_ok &= run_case("Case 4 (a dimension smaller than one whole tile)", 3, 11, 4);
printf("\nall %d cases verified correct: %s\n", 4, all_ok ? "true" : "false");
printf("\n=== Bank-conflict arithmetic on THIS file's tile shapes, Appendix C.3's exact method ===\n");
int naive_banks = compute_distinct_banks(TILE); // tile[32][32]: row stride = 32 floats
int padded_banks = compute_distinct_banks(TILE + 1); // tile[32][33]: row stride = 33 floats
printf("transpose_naive_boundary_checked's tile[%d][%d], column access: %d distinct bank(s) (1 = full 32-way conflict)\n",
TILE, TILE, naive_banks);
printf("transpose_padded_boundary_checked's tile[%d][%d], column access: %d distinct bank(s) (32 = fully resolved)\n",
TILE, TILE + 1, padded_banks);
printf("adding boundary checks changes NOTHING about the bank math -- it depends only on threadIdx and\n");
printf("row stride, never on whether a given thread's data happened to be in-bounds: %s\n",
(naive_banks == 1 && padded_banks == 32) ? "confirmed, matches Appendix C.3 exactly" : "MISMATCH");
return all_ok ? 0 : 1;
}
Worked Example K.2.1 — Four boundary shapes, all genuinely simulated and verified¶
NVDIR=/usr/local/lib/python3.11/dist-packages/nvidia/cu13
nvcc -arch=sm_80 -I$NVDIR/include -L$NVDIR/lib -lcudart 02_matrix_transpose_tiled_boundary.cu -o 02_matrix_transpose_tiled_boundary
LD_LIBRARY_PATH=$NVDIR/lib ./02_matrix_transpose_tiled_boundary
Genuinely compiled and run:
=== Section K.2: tiled transpose with boundary checks, genuinely simulated and verified ===
Case 1 (rows and cols both divide tile_width evenly): 8x8 matrix, tile_width=4, grid=2x2 blocks (boundary tiles: no)
tiled+boundary-checked simulation matches scalar reference, element for element: true
Case 2 (neither dimension divides tile_width evenly): 7x5 matrix, tile_width=4, grid=2x2 blocks (boundary tiles: yes)
tiled+boundary-checked simulation matches scalar reference, element for element: true
Case 3 (larger, still non-divisible, non-square): 10x13 matrix, tile_width=4, grid=4x3 blocks (boundary tiles: yes)
tiled+boundary-checked simulation matches scalar reference, element for element: true
Case 4 (a dimension smaller than one whole tile): 3x11 matrix, tile_width=4, grid=3x1 blocks (boundary tiles: yes)
tiled+boundary-checked simulation matches scalar reference, element for element: true
all 4 cases verified correct: true
=== Bank-conflict arithmetic on THIS file's tile shapes, Appendix C.3's exact method ===
transpose_naive_boundary_checked's tile[32][32], column access: 1 distinct bank(s) (1 = full 32-way conflict)
transpose_padded_boundary_checked's tile[32][33], column access: 32 distinct bank(s) (32 = fully resolved)
adding boundary checks changes NOTHING about the bank math -- it depends only on threadIdx and
row stride, never on whether a given thread's data happened to be in-bounds: confirmed, matches Appendix C.3 exactly
transpose_naive_boundary_checked and transpose_padded_boundary_checked are both genuinely compiled for sm_80 above; consistent with every kernel launch elsewhere in this book, actually running either one needs a real device this sandbox does not have. What the compile step alone cannot prove — that the boundary-checked index math is actually correct at a matrix size that doesn't divide evenly by the tile width, in either dimension, or even a dimension smaller than one whole tile — is proven instead by simulate_tiled_transpose, a host-side function that walks the identical grid/block/tile-step loop nest a real launch would, genuinely executed here and checked element-for-element against a plain scalar transpose reference across four shapes, including one (3x11) where an entire tile row is partially empty in both directions at once.
The bank-conflict arithmetic closing this section reuses Appendix C.3's compute_distinct_banks function completely unchanged, applied to this file's own tile[32][32] and tile[32][33]: one distinct bank (a full 32-way conflict) for the naive tile, 32 distinct banks (fully resolved) for the padded one — the identical result Appendix C.3 found, confirming that adding a boundary check changes nothing about the bank math, because that math is a pure function of threadIdx and the tile's row stride, with no dependence at all on whether the thread doing the reading actually owned real data.
[COMMON TRAP] Appendix E.3's own worked solution already warned that "a real kernel handling a matrix size that doesn't divide evenly needs the same bounds-checking and padding logic" tiling elsewhere in this book already applies — this section is the concrete answer to that warning, not a new idea: the padding fix (a shared-memory bank-conflict concern) and the boundary check (a which-threads-own-real-data concern) are completely independent decisions that happen to both live inside the same kernel, and skipping either one alone — padding without bounds-checking, or bounds-checking without padding — leaves exactly the one problem the skipped fix was for.
K.3 Mapping Matrix Coordinates onto CUDA Grids [FOUNDATIONAL]¶
Intuition¶
Chapter 18.1's ceiling-division-plus-bounds-check pattern was built for a flat, one-dimensional problem. A matrix needs that exact decision made twice — once per axis — plus one decision Chapter 18.1 never had to make at all: which matrix axis becomes threadIdx.x and which becomes threadIdx.y. Every 2D kernel this book has launched since Chapter 18.3 has already picked x for columns and y for rows without pausing to justify it; this section justifies it, and proves the coverage claim exhaustively rather than by argument.
Background¶
LaunchConfig2D is Chapter 18.1's compute_launch_config applied once per axis; verify_exact_coverage genuinely runs every (blockIdx, threadIdx) pair a real launch would produce and counts how many times each logical matrix position gets touched — the same exhaustive-verification discipline Appendix H.3's for_each_index walk and Appendix E.3's simulate_tiled both already rely on in place of hardware this sandbox does not have.
File: 03_grid_mapping_launch_config_2d.cpp¶
#include <cstdio>
#include <vector>
// Appendix K.3 -- Chapter 18.1 built compute_launch_config for a flat,
// one-dimensional problem: ceiling-divide size by threads_per_block, then
// guard the leftover threads with if (idx < size). A matrix needs the
// identical decision made TWICE, once per axis, and it needs one more
// decision Chapter 18.1 never had to make at all: which matrix axis
// becomes threadIdx.x and which becomes threadIdx.y. This file makes both
// decisions explicit, extends Chapter 18.1's launch-config struct to two
// dimensions, and then PROVES -- by genuinely running every (blockIdx,
// threadIdx) pair a real launch would produce, not by arguing it in prose
// -- that the resulting grid covers every element of an M x N matrix
// exactly once, no matter how M and N divide against the chosen block
// shape.
struct LaunchConfig2D {
int blocks_x, blocks_y; // blocks_x covers cols, blocks_y covers rows
int block_dim_x, block_dim_y;
long long total_threads_x, total_threads_y;
long long wasted_x, wasted_y;
};
// Direct two-axis extension of Chapter 18.1's compute_launch_config --
// identical ceiling-division formula, applied once per axis. Convention
// (matching every 2D kernel this book has launched since Chapter 18.3 and
// Appendix C.3/E.3): x indexes COLUMNS, y indexes ROWS.
LaunchConfig2D compute_launch_config_2d(int rows, int cols, int block_dim_x, int block_dim_y) {
LaunchConfig2D cfg;
cfg.block_dim_x = block_dim_x;
cfg.block_dim_y = block_dim_y;
cfg.blocks_x = (cols + block_dim_x - 1) / block_dim_x;
cfg.blocks_y = (rows + block_dim_y - 1) / block_dim_y;
cfg.total_threads_x = (long long)cfg.blocks_x * block_dim_x;
cfg.total_threads_y = (long long)cfg.blocks_y * block_dim_y;
cfg.wasted_x = cfg.total_threads_x - cols;
cfg.wasted_y = cfg.total_threads_y - rows;
return cfg;
}
void print_config(const char* label, int rows, int cols, const LaunchConfig2D& cfg) {
printf("%s: %dx%d matrix, block=%dx%d\n", label, rows, cols, cfg.block_dim_x, cfg.block_dim_y);
printf(" grid = %d x %d blocks, total launched threads = %lld x %lld\n",
cfg.blocks_x, cfg.blocks_y, cfg.total_threads_x, cfg.total_threads_y);
long long total_launched = cfg.total_threads_x * cfg.total_threads_y;
long long total_needed = (long long)rows * cols;
printf(" wasted (out-of-bounds, guarded) threads: %lld along x, %lld along y, %lld total (%lld launched - %lld needed)\n",
cfg.wasted_x, cfg.wasted_y, total_launched - total_needed, total_launched, total_needed);
}
// Genuinely run every (blockIdx.x, blockIdx.y, threadIdx.x, threadIdx.y)
// combination a real <<<grid, block>>> launch would produce, apply
// Chapter 18.1's own bounds-check pattern, and count how many times each
// logical (row, col) position gets touched. Exhaustive, deterministic,
// and -- unlike a device launch -- fully executable in this no-GPU
// environment, exactly the way Appendix H.3's for_each_index walk and
// Appendix E.3's simulate_tiled both stand in for hardware this sandbox
// does not have.
bool verify_exact_coverage(int rows, int cols, const LaunchConfig2D& cfg) {
std::vector<int> touch_count((size_t)rows * cols, 0);
for (int by = 0; by < cfg.blocks_y; by++) {
for (int bx = 0; bx < cfg.blocks_x; bx++) {
for (int ty = 0; ty < cfg.block_dim_y; ty++) {
for (int tx = 0; tx < cfg.block_dim_x; tx++) {
int out_row = by * cfg.block_dim_y + ty; // blockIdx.y * blockDim.y + threadIdx.y
int out_col = bx * cfg.block_dim_x + tx; // blockIdx.x * blockDim.x + threadIdx.x
if (out_row < rows && out_col < cols) {
touch_count[(size_t)out_row * cols + out_col]++;
}
}
}
}
}
long long zero = 0, exactly_one = 0, more_than_one = 0;
for (int v : touch_count) {
if (v == 0) zero++;
else if (v == 1) exactly_one++;
else more_than_one++;
}
printf(" coverage check: %lld positions touched 0 times, %lld touched exactly once, %lld touched more than once\n",
zero, exactly_one, more_than_one);
return zero == 0 && more_than_one == 0 && exactly_one == (long long)rows * cols;
}
int main() {
printf("=== Section K.3: 2D ceiling-division launch configs, exhaustively verified ===\n\n");
struct { int rows, cols, bx, by; const char* label; } cases[] = {
{32, 32, 16, 16, "Case 1 (evenly divisible both axes)"},
{37, 61, 16, 16, "Case 2 (neither axis divides evenly)"},
{5, 200, 16, 16, "Case 3 (one axis much smaller than block)"},
{1000, 1000, 32, 8, "Case 4 (non-square block shape)"},
};
bool all_ok = true;
for (auto& c : cases) {
LaunchConfig2D cfg = compute_launch_config_2d(c.rows, c.cols, c.bx, c.by);
print_config(c.label, c.rows, c.cols, cfg);
all_ok &= verify_exact_coverage(c.rows, c.cols, cfg);
printf("\n");
}
printf("every case covers its matrix exactly once, wasting only the counted leftover threads: %s\n\n",
all_ok ? "true" : "false");
printf("=== Why x maps to columns: consecutive threadIdx.x must produce consecutive addresses ===\n");
int cols = 100, block_dim_x = 32, block_dim_y = 8;
int fixed_block_row = 2, fixed_block_col = 3, fixed_ty = 1;
printf("fixed blockIdx=(%d,%d), threadIdx.y=%d; varying threadIdx.x across one warp (0..7 shown):\n",
fixed_block_col, fixed_block_row, fixed_ty);
printf("CORRECT convention (out_col uses threadIdx.x, out_row uses threadIdx.y):\n");
for (int tx = 0; tx < 8; tx++) {
int out_row = fixed_block_row * block_dim_y + fixed_ty;
int out_col = fixed_block_col * block_dim_x + tx;
long long addr = (long long)out_row * cols + out_col;
printf(" threadIdx.x=%d -> (row=%d, col=%d) -> address %lld\n", tx, out_row, out_col, addr);
}
printf(" addresses for consecutive threadIdx.x are consecutive integers (stride 1): coalesced\n");
printf("\n**[COMMON TRAP]** swapping which axis threadIdx.x indexes still touches every element,\n");
printf("but destroys coalescing -- genuinely computed, not merely asserted:\n");
for (int tx = 0; tx < 8; tx++) {
// Swapped convention: out_row now uses threadIdx.x, out_col uses threadIdx.y.
int out_row = fixed_block_col * block_dim_x + tx;
int out_col = fixed_block_row * block_dim_y + fixed_ty;
long long addr = (long long)out_row * cols + out_col;
printf(" threadIdx.x=%d -> (row=%d, col=%d) -> address %lld\n", tx, out_row, out_col, addr);
}
printf(" addresses for consecutive threadIdx.x now differ by cols=%d, not 1: NOT coalesced\n", cols);
printf(" every element is still visited exactly once (correctness survives), but each of the 8 lanes\n");
printf(" above would issue its own separate memory transaction instead of sharing one\n");
return all_ok ? 0 : 1;
}
Worked Example K.3.1 — Exhaustive coverage across four shapes, and the coalescing reversal genuinely computed¶
g++ -std=c++17 -Wall -Wextra -O2 03_grid_mapping_launch_config_2d.cpp -o 03_grid_mapping_launch_config_2d
./03_grid_mapping_launch_config_2d
Genuinely compiled and run:
=== Section K.3: 2D ceiling-division launch configs, exhaustively verified ===
Case 1 (evenly divisible both axes): 32x32 matrix, block=16x16
grid = 2 x 2 blocks, total launched threads = 32 x 32
wasted (out-of-bounds, guarded) threads: 0 along x, 0 along y, 0 total (1024 launched - 1024 needed)
coverage check: 0 positions touched 0 times, 1024 touched exactly once, 0 touched more than once
Case 2 (neither axis divides evenly): 37x61 matrix, block=16x16
grid = 4 x 3 blocks, total launched threads = 64 x 48
wasted (out-of-bounds, guarded) threads: 3 along x, 11 along y, 815 total (3072 launched - 2257 needed)
coverage check: 0 positions touched 0 times, 2257 touched exactly once, 0 touched more than once
Case 3 (one axis much smaller than block): 5x200 matrix, block=16x16
grid = 13 x 1 blocks, total launched threads = 208 x 16
wasted (out-of-bounds, guarded) threads: 8 along x, 11 along y, 2328 total (3328 launched - 1000 needed)
coverage check: 0 positions touched 0 times, 1000 touched exactly once, 0 touched more than once
Case 4 (non-square block shape): 1000x1000 matrix, block=32x8
grid = 32 x 125 blocks, total launched threads = 1024 x 1000
wasted (out-of-bounds, guarded) threads: 24 along x, 0 along y, 24000 total (1024000 launched - 1000000 needed)
coverage check: 0 positions touched 0 times, 1000000 touched exactly once, 0 touched more than once
every case covers its matrix exactly once, wasting only the counted leftover threads: true
=== Why x maps to columns: consecutive threadIdx.x must produce consecutive addresses ===
fixed blockIdx=(3,2), threadIdx.y=1; varying threadIdx.x across one warp (0..7 shown):
CORRECT convention (out_col uses threadIdx.x, out_row uses threadIdx.y):
threadIdx.x=0 -> (row=17, col=96) -> address 1796
threadIdx.x=1 -> (row=17, col=97) -> address 1797
threadIdx.x=2 -> (row=17, col=98) -> address 1798
threadIdx.x=3 -> (row=17, col=99) -> address 1799
threadIdx.x=4 -> (row=17, col=100) -> address 1800
threadIdx.x=5 -> (row=17, col=101) -> address 1801
threadIdx.x=6 -> (row=17, col=102) -> address 1802
threadIdx.x=7 -> (row=17, col=103) -> address 1803
addresses for consecutive threadIdx.x are consecutive integers (stride 1): coalesced
**[COMMON TRAP]** swapping which axis threadIdx.x indexes still touches every element,
but destroys coalescing -- genuinely computed, not merely asserted:
threadIdx.x=0 -> (row=96, col=17) -> address 9617
threadIdx.x=1 -> (row=97, col=17) -> address 9717
threadIdx.x=2 -> (row=98, col=17) -> address 9817
threadIdx.x=3 -> (row=99, col=17) -> address 9917
threadIdx.x=4 -> (row=100, col=17) -> address 10017
threadIdx.x=5 -> (row=101, col=17) -> address 10117
threadIdx.x=6 -> (row=102, col=17) -> address 10217
threadIdx.x=7 -> (row=103, col=17) -> address 10317
addresses for consecutive threadIdx.x now differ by cols=100, not 1: NOT coalesced
every element is still visited exactly once (correctness survives), but each of the 8 lanes
above would issue its own separate memory transaction instead of sharing one
All four cases — a clean divide, neither axis dividing evenly, one axis far smaller than the block, and a deliberately non-square block shape — cover every logical matrix position exactly once, with the wasted-thread count matching total launched - total needed exactly, extending Chapter 18.1's own 1D accounting to two axes without any new idea. The coalescing worked example then computes real addresses, not a description of them: for a fixed blockIdx and threadIdx.y, consecutive threadIdx.x values under the correct x -> column convention produce consecutive flat addresses (stride 1), while the same warp under the swapped convention produces addresses cols apart — genuinely computed, not asserted, and still touching every element exactly once either way, since correctness and coalescing are independent properties of the same index formula.
K.4 Boundary Checks That Guard Only One Axis [FOUNDATIONAL]¶
Intuition¶
Chapter 13.1's [COMMON TRAP] showed a shape check that simply does not exist anywhere in the code. This section's bug is subtler: the check exists, and checks a real, correct condition — it just checks only one of the two axes a 2D kernel actually needs, an easy mistake to make by pattern-matching Chapter 18.1's 1D if (idx < size) onto two dimensions and forgetting the second half.
Background¶
An earlier draft of this exact section used element-wise addition and discovered something worth keeping rather than discarding: because out[row*cols+col] = a[row*cols+col] + b[row*cols+col] reads and writes through the identical flat-index formula for all three arrays, an overflowing (row, col) aliases to the exact same flat slot — and reads the exact same legitimate data — as the very next row's own first columns. That specific bug turns out to be silently harmless everywhere except the last row. The fill kernel below is chosen instead because its written value depends on the calling thread's own (row, col) identity rather than a re-aliased read, so a genuinely wrong answer is what the missing check produces, and the point generalizes correctly: whether a missing axis check produces silently wrong DATA depends on the specific operation, not on the missing check alone.
File: 04_boundary_check_2d_add_kernel.cu¶
#include <cstdio>
#include <cstdlib>
#include <vector>
// Appendix K.4 -- Chapter 13.1's [COMMON TRAP] showed a shape-check that
// simply doesn't exist (a.cols vs b.rows). This section shows the CUDA
// analogue: a boundary check that exists but is INCOMPLETE -- checking
// only one of a 2D kernel's two axes, a mistake easy to make by pattern-
// matching Chapter 18.1's 1D "if (idx < size)" onto two dimensions and
// forgetting one half. Both kernels below fill each matrix position with
// a distinctive, position-derived marker (row*1000+col) specifically so a
// corrupted position is obviously wrong rather than coincidentally right --
// an earlier draft of this exact section used element-wise addition
// instead, and discovered a genuinely different, subtler fact: because
// out[row*cols+col] = a[row*cols+col] + b[row*cols+col] reads and writes
// through the IDENTICAL flat-index formula for all three arrays, an
// overflowing (row, col) aliases to the exact same flat slot as the
// legitimate (row+1, 0) and reads the exact same legitimate data -- so
// that particular bug turns out to be silently HARMLESS for every row but
// the last. That is a real, worth-knowing fact about element-wise ops
// specifically (Section K.4's own Worked Example below proves it), not a
// property of boundary-check bugs in general -- this section's fill
// kernel is chosen precisely because its written value depends on the
// (possibly wrong) row/col identity itself, not on a re-aliased read, so
// a genuine wrong answer is what a missing check produces here.
//
// Both kernels are genuinely compiled with nvcc for a real architecture;
// neither is launched (this sandbox has no device), so both are checked
// the way K.2 and K.3 already checked their own kernels: a host-side
// simulator walking the identical grid/block/thread loop nest, run for
// real.
__global__ void matrix_fill_boundary_checked(float* out, int rows, int cols) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < rows && col < cols) {
out[row * cols + col] = (float)(row * 1000 + col);
}
}
// The bug: only the ROW half of Chapter 18.1's bounds-check pattern
// survived. For any cols that doesn't divide blockDim.x evenly, threads
// past the true last column are never turned away -- they still compute
// row*cols+col using their OWN (out-of-range) col, so the marker they
// write is wrong for whatever flat slot that overshoot actually lands on.
__global__ void matrix_fill_buggy_row_only_checked(float* out, int rows, int cols) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < rows) { // BUG: col is never checked against cols at all
out[row * cols + col] = (float)(row * 1000 + col);
}
}
// ---- Host C++ below: genuinely compiled AND genuinely run ----
void fill_reference(int rows, int cols, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
out[(size_t)r * cols + c] = (float)(r * 1000 + c);
}
// Simulates matrix_fill_boundary_checked's exact grid/block/thread loop nest.
void simulate_correct_fill(int rows, int cols, int block_dim_x, int block_dim_y, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
int blocks_x = (cols + block_dim_x - 1) / block_dim_x;
int blocks_y = (rows + block_dim_y - 1) / block_dim_y;
for (int by = 0; by < blocks_y; by++)
for (int bx = 0; bx < blocks_x; bx++)
for (int ty = 0; ty < block_dim_y; ty++)
for (int tx = 0; tx < block_dim_x; tx++) {
int row = by * block_dim_y + ty;
int col = bx * block_dim_x + tx;
if (row < rows && col < cols) out[(size_t)row * cols + col] = (float)(row * 1000 + col);
}
}
// Simulates matrix_fill_buggy_row_only_checked's exact (buggy) loop nest.
// `out` is sized (rows+1) x cols -- one whole extra PADDING row, filled
// with a known sentinel, rather than exactly rows x cols. This is Chapter
// 13.1's own sentinel technique (Worked Example 13.1.4): it turns what
// would otherwise be a genuine, undefined-behavior heap overflow on the
// very last real row into a real, safe, deterministic write onto memory
// this file legitimately owns, so the corruption is observed directly
// rather than merely asserted to be dangerous.
//
// process_blocks_in_reverse controls only the ORDER this host simulation
// visits blockIdx.y in -- CUDA itself makes zero guarantees about block
// execution order or concurrency (this book has relied on that fact since
// Chapter 18's own launch-configuration discussion), so a real GPU could
// legally run these blocks in either order, both orders interleaved, or
// fully in parallel. This parameter exists specifically to demonstrate
// that the bug's VISIBLE damage depends entirely on an order nothing
// promises you.
void simulate_buggy_fill(int rows, int cols, int block_dim_x, int block_dim_y,
bool process_blocks_in_reverse, std::vector<float>& out) {
int blocks_x = (cols + block_dim_x - 1) / block_dim_x;
int blocks_y = (rows + block_dim_y - 1) / block_dim_y;
for (int by_i = 0; by_i < blocks_y; by_i++) {
int by = process_blocks_in_reverse ? (blocks_y - 1 - by_i) : by_i;
for (int bx = 0; bx < blocks_x; bx++)
for (int ty = 0; ty < block_dim_y; ty++)
for (int tx = 0; tx < block_dim_x; tx++) {
int row = by * block_dim_y + ty;
int col = bx * block_dim_x + tx;
if (row < rows) { // the bug, reproduced exactly
out[(size_t)row * cols + col] = (float)(row * 1000 + col);
}
}
}
}
int main() {
printf("=== Section K.4: a boundary check that guards only ONE of two axes ===\n\n");
int rows = 4, cols = 10, block_dim_x = 4, block_dim_y = 1;
printf("matrix: %dx%d, block=%dx%d (cols=%d is NOT a multiple of block_dim_x=%d)\n\n",
rows, cols, block_dim_x, block_dim_y, cols, block_dim_x);
std::vector<float> reference;
fill_reference(rows, cols, reference);
printf("=== The CORRECT kernel, fully boundary-checked ===\n");
std::vector<float> correct_out;
simulate_correct_fill(rows, cols, block_dim_x, block_dim_y, correct_out);
bool correct_matches = (correct_out == reference);
printf("matches the row*1000+col reference, element for element: %s\n\n", correct_matches ? "true" : "false");
printf("=== Worked Example: why the SAME bug on ELEMENT-WISE ADDITION would be harmless ===\n");
printf("out[row*cols+col] = a[row*cols+col] + b[row*cols+col] reads a and b through the SAME\n");
printf("flat formula it writes out through. An overflowing (row=0, col=cols=%d) computes flat index\n", cols);
printf("%d, which is ALSO exactly (row=1, col=0)'s own flat index in every one of a, b, and out --\n", cols);
printf("so the value written is a(1,0)+b(1,0), the mathematically correct value for that slot,\n");
printf("computed one iteration early -- the SAME identity every interior overflow relies on.\n");
printf("This fill kernel's written value depends on the CALLING thread's own (row, col), not on a\n");
printf("re-aliased read, so it has no such accidental self-correction -- the corruption below is real.\n\n");
printf("=== The BUGGY fill kernel: row checked, column never checked ===\n");
const float SENTINEL = -999.0f;
int padded_rows = rows + 1; // one extra real row of storage, filled with a sentinel
int overflow_cols_per_row = (block_dim_x - (cols % block_dim_x)) % block_dim_x;
printf("overflowing columns per row (col >= cols but still written): %d\n\n", overflow_cols_per_row);
std::vector<float> buggy_forward((size_t)padded_rows * cols, SENTINEL);
std::vector<float> buggy_reverse((size_t)padded_rows * cols, SENTINEL);
simulate_buggy_fill(rows, cols, block_dim_x, block_dim_y, false, buggy_forward);
simulate_buggy_fill(rows, cols, block_dim_x, block_dim_y, true, buggy_reverse);
auto count_corrupted_rows = [&](const std::vector<float>& out) {
int corrupted = 0;
for (int r = 1; r < rows; r++) {
bool bad = (out[(size_t)r * cols + 0] != reference[(size_t)r * cols + 0]) ||
(out[(size_t)r * cols + 1] != reference[(size_t)r * cols + 1]);
if (bad) corrupted++;
}
return corrupted;
};
auto count_sentinel_overwritten = [&](const std::vector<float>& out) {
int n = 0;
for (int c = 0; c < cols; c++) if (out[(size_t)rows * cols + c] != SENTINEL) n++;
return n;
};
int fwd_corrupted_rows = count_corrupted_rows(buggy_forward);
int rev_corrupted_rows = count_corrupted_rows(buggy_reverse);
int fwd_sentinel = count_sentinel_overwritten(buggy_forward);
int rev_sentinel = count_sentinel_overwritten(buggy_reverse);
printf("processing blockIdx.y forward (0, 1, 2, ...): %d of %d interior rows visibly corrupted\n",
fwd_corrupted_rows, rows - 1);
printf("processing blockIdx.y in reverse (..., 2, 1, 0): %d of %d interior rows visibly corrupted\n",
rev_corrupted_rows, rows - 1);
printf("padding-row sentinel cells overwritten, forward order: %d\n", fwd_sentinel);
printf("padding-row sentinel cells overwritten, reverse order: %d\n\n", rev_sentinel);
printf("row-by-row detail, reverse order (the one that DOES show interior corruption):\n");
for (int r = 0; r < rows; r++) {
float got0 = buggy_reverse[(size_t)r * cols + 0];
float want0 = reference[(size_t)r * cols + 0];
printf(" row %d, col 0: wrote %.0f, correct is %.0f -- %s\n",
r, got0, want0, (got0 == want0) ? "correct" : "WRONG (overwritten by the row above's overflow)");
}
printf("\n**[COMMON TRAP]** in FORWARD order, each row's own correct write happens AFTER the row\n");
printf("above's overflow lands on it, so the correct value silently overwrites the bad one -- every\n");
printf("interior row looks fine by accident. In REVERSE order, each row's own correct write happens\n");
printf("BEFORE the row above's overflow arrives, so the bad value survives and every interior row is\n");
printf("visibly wrong. CUDA guarantees neither order (nor that blocks run one at a time at all) --\n");
printf("matrix_fill_buggy_row_only_checked compiles cleanly and produces a genuine, scheduling-\n");
printf("dependent data race, not a fixed wrong answer. The one piece of damage every order agrees on\n");
printf("is the final row's overflow, because nothing ever writes to that memory a second time to\n");
printf("cover it back up -- exactly the %d sentinel cells both orders show above.\n", overflow_cols_per_row);
bool orders_genuinely_differ = (fwd_corrupted_rows != rev_corrupted_rows);
bool forward_self_heals = (fwd_corrupted_rows == 0);
bool reverse_fully_corrupts = (rev_corrupted_rows == rows - 1);
bool last_row_overflow_always_visible = (fwd_sentinel == overflow_cols_per_row) && (rev_sentinel == overflow_cols_per_row);
bool all_ok = correct_matches && orders_genuinely_differ && forward_self_heals &&
reverse_fully_corrupts && last_row_overflow_always_visible;
printf("\nall counted quantities are internally consistent: %s\n", all_ok ? "true" : "false");
return all_ok ? 0 : 1;
}
Worked Example K.4.1 — A genuine, scheduling-dependent data race, run both ways¶
NVDIR=/usr/local/lib/python3.11/dist-packages/nvidia/cu13
nvcc -arch=sm_80 -I$NVDIR/include -L$NVDIR/lib -lcudart 04_boundary_check_2d_add_kernel.cu -o 04_boundary_check_2d_add_kernel
LD_LIBRARY_PATH=$NVDIR/lib ./04_boundary_check_2d_add_kernel
Genuinely compiled and run:
=== Section K.4: a boundary check that guards only ONE of two axes ===
matrix: 4x10, block=4x1 (cols=10 is NOT a multiple of block_dim_x=4)
=== The CORRECT kernel, fully boundary-checked ===
matches the row*1000+col reference, element for element: true
=== Worked Example: why the SAME bug on ELEMENT-WISE ADDITION would be harmless ===
out[row*cols+col] = a[row*cols+col] + b[row*cols+col] reads a and b through the SAME
flat formula it writes out through. An overflowing (row=0, col=cols=10) computes flat index
10, which is ALSO exactly (row=1, col=0)'s own flat index in every one of a, b, and out --
so the value written is a(1,0)+b(1,0), the mathematically correct value for that slot,
computed one iteration early -- the SAME identity every interior overflow relies on.
This fill kernel's written value depends on the CALLING thread's own (row, col), not on a
re-aliased read, so it has no such accidental self-correction -- the corruption below is real.
=== The BUGGY fill kernel: row checked, column never checked ===
overflowing columns per row (col >= cols but still written): 2
processing blockIdx.y forward (0, 1, 2, ...): 0 of 3 interior rows visibly corrupted
processing blockIdx.y in reverse (..., 2, 1, 0): 3 of 3 interior rows visibly corrupted
padding-row sentinel cells overwritten, forward order: 2
padding-row sentinel cells overwritten, reverse order: 2
row-by-row detail, reverse order (the one that DOES show interior corruption):
row 0, col 0: wrote 0, correct is 0 -- correct
row 1, col 0: wrote 10, correct is 1000 -- WRONG (overwritten by the row above's overflow)
row 2, col 0: wrote 1010, correct is 2000 -- WRONG (overwritten by the row above's overflow)
row 3, col 0: wrote 2010, correct is 3000 -- WRONG (overwritten by the row above's overflow)
**[COMMON TRAP]** in FORWARD order, each row's own correct write happens AFTER the row
above's overflow lands on it, so the correct value silently overwrites the bad one -- every
interior row looks fine by accident. In REVERSE order, each row's own correct write happens
BEFORE the row above's overflow arrives, so the bad value survives and every interior row is
visibly wrong. CUDA guarantees neither order (nor that blocks run one at a time at all) --
matrix_fill_buggy_row_only_checked compiles cleanly and produces a genuine, scheduling-
dependent data race, not a fixed wrong answer. The one piece of damage every order agrees on
is the final row's overflow, because nothing ever writes to that memory a second time to
cover it back up -- exactly the 2 sentinel cells both orders show above.
all counted quantities are internally consistent: true
CUDA guarantees nothing at all about the order or concurrency of block execution — this book has relied on that fact since Chapter 18's own launch-configuration discussion. This file's host simulation genuinely runs the identical buggy loop nest in two different, both entirely legal, orders. Processing blockIdx.y forward happens to let each interior row's own correct write land AFTER the row above's overflow already wrote something wrong there, silently overwriting the bad value with the right one — every interior row looks fine, by accident. Processing the exact same blocks in reverse flips that timing: each row's own correct write now happens BEFORE the row above's overflow arrives, so the wrong value survives and every interior row is visibly corrupted. The one piece of damage both orders agree on is the final row's overflow, landing on memory the correct kernel would never touch at all (locked here, following Chapter 13.1's own sentinel technique, as a padding row rather than a genuine out-of-bounds heap write) — because nothing in either order ever gets a second chance to write over it.
[COMMON TRAP] matrix_fill_buggy_row_only_checked compiles without a single warning and would launch without error on real hardware. Its bug is not a crash and not a fixed wrong answer — it is a genuine data race whose visible symptoms depend on a scheduling order this book's own Chapter 18 discussion already established CUDA never promises. A missing axis check that happens to look harmless in one test run (this section's own forward-order case shows zero visible interior corruption) is not evidence the check was unnecessary; it is evidence the test happened to run in the one order that hides the bug.
K.5 Row-Major Kernels Meeting Column-Major Data [FOUNDATIONAL]¶
Intuition¶
Chapter 7.2 already named the collision this section makes concrete: "cuBLAS, the CUDA math library every serious GEMM in this book eventually calls into, is column-major by convention." Every coalescing argument this appendix has made so far — including Section K.3's own x -> column convention — was never actually a rule about columns specifically; it was a rule about whichever axis is contiguous in memory. Column-major data simply has a different axis be the contiguous one, and the entire convention flips accordingly.
Background¶
RowMajorMatrix and ColMajorMatrix are Chapter 6 and Chapter 7.2's own two offset formulas, specialized to two dimensions exactly the way Chapter 13's Matrix::get already specialized row-major alone. Converting between them is nothing more than reading through one formula and writing through the other, for every logical (row, col) — and reinterpreting one layout's raw bytes as the other without that conversion is Chapter 13.3's naive-reshape trap, encountered here in the one place this book will actually hit it in practice.
File: 05_row_major_col_major_kernel_interop.cpp¶
#include <cstdio>
#include <cstdlib>
#include <vector>
// Appendix K.5 -- Chapter 7.2 named the one place this book's own row-major
// convention meets a genuinely column-major world head-on: "cuBLAS, the
// CUDA math library every serious GEMM in this book eventually calls into,
// is column-major by convention." K.3 established that x should index the
// FASTEST-VARYING (contiguous) axis for coalescing. For row-major data
// that axis is the column; for column-major data it's the ROW -- the exact
// same coalescing argument, pointed at the opposite axis, because the
// underlying rule was never "x means column," it was "x means whichever
// axis is actually contiguous in memory." This section makes that
// reversal concrete, builds a genuine row-major <-> column-major converter
// for handing data to (or receiving it from) a column-major API, and
// closes with the naive-reinterpretation trap Chapter 13.3 already proved
// is a silent transpose, now in the specific context this book will
// actually hit it: a cuBLAS call.
struct RowMajorMatrix {
std::vector<float> data;
int rows, cols;
RowMajorMatrix(int rows_, int cols_) : data((size_t)rows_ * cols_, 0.0f), rows(rows_), cols(cols_) {}
float get(int r, int c) const { return data[(size_t)r * cols + c]; } // row-major, Chapter 3
void set(int r, int c, float v) { data[(size_t)r * cols + c] = v; }
};
struct ColMajorMatrix {
std::vector<float> data;
int rows, cols;
ColMajorMatrix(int rows_, int cols_) : data((size_t)rows_ * cols_, 0.0f), rows(rows_), cols(cols_) {}
float get(int r, int c) const { return data[(size_t)c * rows + r]; } // column-major, Chapter 7.2
void set(int r, int c, float v) { data[(size_t)c * rows + r] = v; }
};
int main() {
printf("=== Section K.5: which axis is contiguous REVERSES between the two layouts ===\n\n");
int rows = 4, cols = 6;
printf("for a %dx%d matrix:\n", rows, cols);
printf(" row-major: stepping one COLUMN moves 1 flat slot; stepping one ROW moves %d slots\n", cols);
printf(" col-major: stepping one ROW moves 1 flat slot; stepping one COLUMN moves %d slots\n\n", rows);
printf("K.3's coalescing rule ('x should index whichever axis is contiguous') therefore assigns:\n");
printf(" row-major kernel: threadIdx.x -> column, threadIdx.y -> row (K.3's own convention)\n");
printf(" column-major kernel: threadIdx.x -> row, threadIdx.y -> column (the exact reversal)\n\n");
// Genuinely confirm the reversal with real address arithmetic, exactly
// the way K.3 confirmed its own row-major convention.
RowMajorMatrix rm(rows, cols);
ColMajorMatrix cm(rows, cols);
printf("row-major: consecutive columns at fixed row, addresses: ");
for (int c = 0; c < cols; c++) printf("%d ", 2 * cols + c); // row=2 fixed
printf(" (stride 1)\n");
printf("col-major: consecutive rows at fixed column, addresses: ");
for (int r = 0; r < rows; r++) printf("%d ", 2 * rows + r); // col=2 fixed
printf(" (stride 1)\n");
printf("row-major: consecutive rows at fixed column, addresses: ");
for (int r = 0; r < rows; r++) printf("%d ", r * cols + 2); // col=2 fixed
printf(" (stride %d)\n", cols);
printf("col-major: consecutive columns at fixed row, addresses: ");
for (int c = 0; c < cols; c++) printf("%d ", c * rows + 2); // row=2 fixed
printf(" (stride %d)\n\n", rows);
// Fill a row-major matrix with a distinctive per-position marker --
// this book's own default Tensor layout since Chapter 6.
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
rm.set(r, c, (float)(r * 100 + c));
printf("=== Converting row-major data (this book's Tensor default) for a column-major API ===\n");
// Two independently-written conversions, cross-checked against each
// other -- this book's own convention (e.g. Worked Example 13.1.1's
// scalar/SIMD cross-check) for trusting a result.
ColMajorMatrix converted_a(rows, cols);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
converted_a.set(r, c, rm.get(r, c)); // logical (r,c) copy, layout handled by set()/get()
ColMajorMatrix converted_b(rows, cols);
for (size_t flat_rm = 0; flat_rm < rm.data.size(); flat_rm++) {
int r = (int)(flat_rm / cols);
int c = (int)(flat_rm % cols);
converted_b.data[(size_t)c * rows + r] = rm.data[flat_rm]; // direct flat-index conversion
}
bool converters_agree = (converted_a.data == converted_b.data);
bool logical_values_preserved = true;
for (int r = 0; r < rows && logical_values_preserved; r++)
for (int c = 0; c < cols; c++)
if (converted_a.get(r, c) != rm.get(r, c)) { logical_values_preserved = false; break; }
printf("two independent row-major -> column-major converters agree, byte for byte: %s\n",
converters_agree ? "true" : "false");
printf("every logical (row, col) value is preserved across the layout change: %s\n\n",
logical_values_preserved ? "true" : "false");
printf("=== Round trip: column-major -> row-major -> column-major returns the original exactly ===\n");
RowMajorMatrix back_to_row(rows, cols);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
back_to_row.set(r, c, converted_a.get(r, c));
bool round_trip_matches = (back_to_row.data == rm.data);
printf("row-major -> column-major -> row-major reproduces the original flat buffer exactly: %s\n\n",
round_trip_matches ? "true" : "false");
printf("=== [COMMON TRAP] reinterpreting a column-major buffer AS row-major without converting ===\n");
printf("This is Chapter 13.3's own naive-reshape trap, hit again here for a different reason: a\n");
printf("cuBLAS call handed row-major bytes without conversion (or a caller that reads cuBLAS's real\n");
printf("column-major OUTPUT as if it were row-major) doesn't crash and doesn't raise an error -- it\n");
printf("silently reads the TRANSPOSE of the intended matrix, because reinterpreting a flat buffer's\n");
printf("layout is exactly a relabeling of which index is contiguous, and swapping which index is\n");
printf("contiguous is precisely what transpose means.\n\n");
RowMajorMatrix naive_reinterpretation(cols, rows); // note: dims swapped, since that's what the bug does
for (size_t i = 0; i < naive_reinterpretation.data.size(); i++)
naive_reinterpretation.data[i] = converted_a.data[i]; // same bytes, wrong layout tag
bool matches_transpose = true;
for (int r = 0; r < rows && matches_transpose; r++)
for (int c = 0; c < cols; c++)
if (naive_reinterpretation.get(c, r) != rm.get(r, c)) { matches_transpose = false; break; }
printf("naive_reinterpretation.get(c, r) genuinely equals rm.get(r, c) for every (r, c): %s\n",
matches_transpose ? "true" : "false");
printf("i.e. reading the column-major buffer's bytes as row-major recovers exactly A-transpose,\n");
printf("never A itself, silently and without any error -- the identical failure mode Chapter 13.3\n");
printf("already proved for a naive reshape of a transposed buffer, now shown for a layout mismatch.\n");
bool all_ok = converters_agree && logical_values_preserved && round_trip_matches && matches_transpose;
printf("\nall counted/checked quantities are internally consistent: %s\n", all_ok ? "true" : "false");
return all_ok ? 0 : 1;
}
Worked Example K.5.1 — Conversion, round trip, and the transpose that appears when you skip it¶
g++ -std=c++17 -Wall -Wextra -O2 05_row_major_col_major_kernel_interop.cpp -o 05_row_major_col_major_kernel_interop
./05_row_major_col_major_kernel_interop
Genuinely compiled and run:
=== Section K.5: which axis is contiguous REVERSES between the two layouts ===
for a 4x6 matrix:
row-major: stepping one COLUMN moves 1 flat slot; stepping one ROW moves 6 slots
col-major: stepping one ROW moves 1 flat slot; stepping one COLUMN moves 4 slots
K.3's coalescing rule ('x should index whichever axis is contiguous') therefore assigns:
row-major kernel: threadIdx.x -> column, threadIdx.y -> row (K.3's own convention)
column-major kernel: threadIdx.x -> row, threadIdx.y -> column (the exact reversal)
row-major: consecutive columns at fixed row, addresses: 12 13 14 15 16 17 (stride 1)
col-major: consecutive rows at fixed column, addresses: 8 9 10 11 (stride 1)
row-major: consecutive rows at fixed column, addresses: 2 8 14 20 (stride 6)
col-major: consecutive columns at fixed row, addresses: 2 6 10 14 18 22 (stride 4)
=== Converting row-major data (this book's Tensor default) for a column-major API ===
two independent row-major -> column-major converters agree, byte for byte: true
every logical (row, col) value is preserved across the layout change: true
=== Round trip: column-major -> row-major -> column-major returns the original exactly ===
row-major -> column-major -> row-major reproduces the original flat buffer exactly: true
=== [COMMON TRAP] reinterpreting a column-major buffer AS row-major without converting ===
This is Chapter 13.3's own naive-reshape trap, hit again here for a different reason: a
cuBLAS call handed row-major bytes without conversion (or a caller that reads cuBLAS's real
column-major OUTPUT as if it were row-major) doesn't crash and doesn't raise an error -- it
silently reads the TRANSPOSE of the intended matrix, because reinterpreting a flat buffer's
layout is exactly a relabeling of which index is contiguous, and swapping which index is
contiguous is precisely what transpose means.
naive_reinterpretation.get(c, r) genuinely equals rm.get(r, c) for every (r, c): true
i.e. reading the column-major buffer's bytes as row-major recovers exactly A-transpose,
never A itself, silently and without any error -- the identical failure mode Chapter 13.3
already proved for a naive reshape of a transposed buffer, now shown for a layout mismatch.
all counted/checked quantities are internally consistent: true
Two independently-written converters — one going through each struct's own get/set, one operating directly on flat indices — agree byte for byte, and every logical value survives a full row-major to column-major and back round trip exactly. The closing trap reuses Chapter 13.3's own already-proven fact rather than a new one: reinterpreting a flat buffer's bytes under the other layout's offset formula is exactly a relabeling of which index is contiguous, and swapping which index is contiguous is precisely what transpose means — genuinely confirmed here by checking naive_reinterpretation.get(c, r) == rm.get(r, c) for every position, not merely asserting it.
K.6 Spiral Traversal: A Strictly Sequential Order With an O(1) Closed Form [FOUNDATIONAL]¶
Intuition¶
Every traversal this appendix has built a CUDA kernel for so far — row-major, column-major, the tiled transpose's block order — has had a trivially linear index formula, row * cols + col or its column-major mirror. Spiral order looks like it cannot have one at all: knowing that a position is visited fifth depends on already knowing which four positions came before it, since the boundary of "the current ring" shrinks as the walk proceeds. That would make spiral order impossible to hand to an independent CUDA thread — one thread cannot ask "what is my rank" without first simulating every position visited before it. This section proves that appearance wrong: spiral order does have a closed form, derived from nothing more than a position's distance from the nearest edge, and once that form exists, every thread can compute its own rank with no dependency on any other thread at all.
Background¶
A position's ring number is k = min(row, col, rows-1-row, cols-1-col) — literally its distance from whichever edge is closest. Peeling away k full rings removes exactly rows*cols - (rows-2k)*(cols-2k) positions, a closed form with no loop, so that quantity alone is every ring-k position's shared starting rank. From there, a position's offset within its own ring depends only on which of the ring's four sides it sits on — a direct subtraction against that side's own fixed coordinate, never a re-walk of the ring itself.
File: 06_spiral_traversal.cu¶
#include <cstdio>
#include <vector>
#include <utility>
// Appendix K.6 -- every traversal order so far in this appendix (row-major,
// column-major, the tiled transpose's block order) has had a trivial,
// already-linear index formula: row*cols+col or its column-major mirror.
// Spiral order does not: which position a loop visits 5th depends on
// having already visited four others first, peeling the matrix's outer
// ring before the next ring in is even touched. That looks, at first,
// like a traversal with no closed form at all -- and therefore no way to
// hand one (row, col) to one CUDA thread and ask "what is your rank in
// this order" without that thread first simulating the whole peel. This
// section derives the closed form anyway, cross-checks it exhaustively
// against the straightforward sequential simulation, and then uses it to
// write a real CUDA kernel where every thread computes its own spiral
// rank independently -- no thread ever needs to know what any other
// thread's rank is.
// The textbook approach: shrink four boundaries inward, appending
// positions as each side of the current ring is walked. Genuinely
// sequential -- each appended position depends on the loop having
// already reached this exact point in this exact order.
std::vector<std::pair<int,int>> spiral_order_sequential(int rows, int cols) {
std::vector<std::pair<int,int>> order;
if (rows <= 0 || cols <= 0) return order;
int top = 0, bottom = rows - 1, left = 0, right = cols - 1;
while (top <= bottom && left <= right) {
for (int c = left; c <= right; c++) order.push_back({top, c});
top++;
for (int r = top; r <= bottom; r++) order.push_back({r, right});
right--;
if (top <= bottom) {
for (int c = right; c >= left; c--) order.push_back({bottom, c});
bottom--;
}
if (left <= right) {
for (int r = bottom; r >= top; r--) order.push_back({r, left});
left++;
}
}
return order;
}
// The closed form. Every position's RING is how far it sits from the
// nearest edge: k = min(row, col, rows-1-row, cols-1-col). Peeling k
// rings removes exactly rows*cols - (rows-2k)*(cols-2k) elements, so
// that quantity alone is every ring-k position's shared starting rank --
// computed in O(1), no simulation of the outer rings required. Within
// ring k, a position's offset depends only on which of the ring's four
// sides it sits on, each a direct subtraction against that side's own
// fixed coordinate.
long long spiral_rank_closed_form(int row, int col, int rows, int cols) {
int k = row;
if (col < k) k = col;
if (rows - 1 - row < k) k = rows - 1 - row;
if (cols - 1 - col < k) k = cols - 1 - col;
int top = k, bottom = rows - 1 - k, left = k, right = cols - 1 - k;
long long ring_start = (long long)rows * cols - (long long)(rows - 2 * k) * (cols - 2 * k);
int width = right - left + 1;
int height = bottom - top + 1;
if (height == 1) { // single remaining row: pure left-to-right
return ring_start + (col - left);
}
if (width == 1) { // single remaining column: pure top-to-bottom
return ring_start + (row - top);
}
if (row == top) { // top side
return ring_start + (col - left);
} else if (col == right) { // right side
return ring_start + (width - 1) + (row - top);
} else if (row == bottom) { // bottom side
return ring_start + (width - 1) + (height - 1) + (right - col);
} else { // left side
return ring_start + 2 * (width - 1) + (height - 1) + (bottom - row);
}
}
bool cross_check_shape(int rows, int cols) {
auto sequential = spiral_order_sequential(rows, cols);
std::vector<long long> closed_form_rank((size_t)rows * cols);
bool all_in_range = true;
for (int r = 0; r < rows; r++) {
for (int c = 0; c < cols; c++) {
long long rank = spiral_rank_closed_form(r, c, rows, cols);
if (rank < 0 || rank >= (long long)rows * cols) { all_in_range = false; continue; }
closed_form_rank[(size_t)r * cols + c] = rank;
}
}
bool matches = all_in_range;
if (all_in_range) {
for (size_t rank = 0; rank < sequential.size() && matches; rank++) {
auto [r, c] = sequential[rank];
if (closed_form_rank[(size_t)r * cols + c] != (long long)rank) matches = false;
}
}
printf(" %dx%d: closed form matches sequential simulation for all %d positions: %s\n",
rows, cols, rows * cols, matches ? "true" : "false");
return matches;
}
// ---- The CUDA kernel: every thread computes its OWN rank, independently ----
__global__ void spiral_scatter_kernel(float* out, const float* in, int rows, int cols) {
int row = blockIdx.y * blockDim.y + threadIdx.y; // K.3's own convention
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < rows && col < cols) {
int k = row;
if (col < k) k = col;
if (rows - 1 - row < k) k = rows - 1 - row;
if (cols - 1 - col < k) k = cols - 1 - col;
int top = k, bottom = rows - 1 - k, left = k, right = cols - 1 - k;
long long ring_start = (long long)rows * cols - (long long)(rows - 2 * k) * (cols - 2 * k);
int width = right - left + 1;
int height = bottom - top + 1;
long long rank;
if (height == 1) rank = ring_start + (col - left);
else if (width == 1) rank = ring_start + (row - top);
else if (row == top) rank = ring_start + (col - left);
else if (col == right) rank = ring_start + (width - 1) + (row - top);
else if (row == bottom) rank = ring_start + (width - 1) + (height - 1) + (right - col);
else rank = ring_start + 2 * (width - 1) + (height - 1) + (bottom - row);
out[rank] = in[row * cols + col];
}
}
// Host-side simulation of the exact grid/block/thread loop nest above --
// Section K.2 and K.3's own verification technique, applied here.
void simulate_spiral_scatter(const std::vector<float>& in, int rows, int cols,
int block_dim_x, int block_dim_y, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
int blocks_x = (cols + block_dim_x - 1) / block_dim_x;
int blocks_y = (rows + block_dim_y - 1) / block_dim_y;
for (int by = 0; by < blocks_y; by++)
for (int bx = 0; bx < blocks_x; bx++)
for (int ty = 0; ty < block_dim_y; ty++)
for (int tx = 0; tx < block_dim_x; tx++) {
int row = by * block_dim_y + ty;
int col = bx * block_dim_x + tx;
if (row < rows && col < cols) {
long long rank = spiral_rank_closed_form(row, col, rows, cols);
out[rank] = in[(size_t)row * cols + col];
}
}
}
int main() {
printf("=== Section K.6: spiral order's closed form, cross-checked exhaustively ===\n\n");
printf("Canonical 3x3 check (values 1..9), sequential order:\n");
auto canon = spiral_order_sequential(3, 3);
printf(" ");
for (auto& [r, c] : canon) printf("%d ", r * 3 + c + 1);
printf("\n (expected: 1 2 3 6 9 8 7 4 5)\n\n");
bool all_shapes_ok = true;
int shapes[][2] = {{1,1},{1,5},{5,1},{3,3},{4,4},{4,5},{5,4},{7,7},{6,9},{9,6}};
for (auto& s : shapes) all_shapes_ok &= cross_check_shape(s[0], s[1]);
printf("\nall %d shapes: closed form and sequential simulation agree everywhere: %s\n\n",
(int)(sizeof(shapes) / sizeof(shapes[0])), all_shapes_ok ? "true" : "false");
printf("=== The CUDA kernel: genuinely simulated across divisible and non-divisible shapes ===\n");
bool all_kernel_ok = true;
struct { int rows, cols, bx, by; } kernel_cases[] = {
{4, 4, 4, 4}, {7, 5, 4, 4}, {9, 6, 4, 4}, {6, 9, 4, 4},
};
for (auto& kc : kernel_cases) {
std::vector<float> in((size_t)kc.rows * kc.cols);
for (size_t i = 0; i < in.size(); i++) in[i] = (float)i;
std::vector<float> kernel_out;
simulate_spiral_scatter(in, kc.rows, kc.cols, kc.bx, kc.by, kernel_out);
auto sequential = spiral_order_sequential(kc.rows, kc.cols);
std::vector<float> reference((size_t)kc.rows * kc.cols);
for (size_t rank = 0; rank < sequential.size(); rank++) {
auto [r, c] = sequential[rank];
reference[rank] = in[(size_t)r * kc.cols + c];
}
bool match = (kernel_out == reference);
all_kernel_ok &= match;
printf(" %dx%d, block=%dx%d: kernel simulation matches sequential-derived reference: %s\n",
kc.rows, kc.cols, kc.bx, kc.by, match ? "true" : "false");
}
printf("\nall kernel cases verified correct: %s\n", all_kernel_ok ? "true" : "false");
bool all_ok = all_shapes_ok && all_kernel_ok;
return all_ok ? 0 : 1;
}
Worked Example K.6.1 — Ten shapes, exhaustively cross-checked, then a genuine CUDA scatter¶
NVDIR=/usr/local/lib/python3.11/dist-packages/nvidia/cu13
nvcc -arch=sm_80 -I$NVDIR/include -L$NVDIR/lib -lcudart 06_spiral_traversal.cu -o 06_spiral_traversal
LD_LIBRARY_PATH=$NVDIR/lib ./06_spiral_traversal
Genuinely compiled and run:
=== Section K.6: spiral order's closed form, cross-checked exhaustively ===
Canonical 3x3 check (values 1..9), sequential order:
1 2 3 6 9 8 7 4 5
(expected: 1 2 3 6 9 8 7 4 5)
1x1: closed form matches sequential simulation for all 1 positions: true
1x5: closed form matches sequential simulation for all 5 positions: true
5x1: closed form matches sequential simulation for all 5 positions: true
3x3: closed form matches sequential simulation for all 9 positions: true
4x4: closed form matches sequential simulation for all 16 positions: true
4x5: closed form matches sequential simulation for all 20 positions: true
5x4: closed form matches sequential simulation for all 20 positions: true
7x7: closed form matches sequential simulation for all 49 positions: true
6x9: closed form matches sequential simulation for all 54 positions: true
9x6: closed form matches sequential simulation for all 54 positions: true
all 10 shapes: closed form and sequential simulation agree everywhere: true
=== The CUDA kernel: genuinely simulated across divisible and non-divisible shapes ===
4x4, block=4x4: kernel simulation matches sequential-derived reference: true
7x5, block=4x4: kernel simulation matches sequential-derived reference: true
9x6, block=4x4: kernel simulation matches sequential-derived reference: true
6x9, block=4x4: kernel simulation matches sequential-derived reference: true
all kernel cases verified correct: true
The canonical 3x3 case (1 2 3 6 9 8 7 4 5) confirms the sequential simulation against the answer every reader of this problem already expects, before the closed form is trusted against anything. Ten shapes — including a single row, a single column, and both non-square orientations — show the closed form and the sequential simulation agreeing at every one of their positions, not merely on the final count. spiral_scatter_kernel is genuinely compiled for sm_80 above; actually launching it needs a device this sandbox does not have, so simulate_spiral_scatter — a host-side walk of the kernel's own exact grid/block/thread loop nest, Section K.2 and K.3's own verification technique — checks it instead, across both divisible and non-divisible block shapes, against a reference built from the independently-verified sequential order.
[COMMON TRAP] A traversal looking sequential in its most natural description — "peel the outer ring, then the next ring in" — does not mean it IS sequential in the sense that matters for parallelizing it. What actually blocks parallelization is a rank that cannot be computed without consulting other positions' ranks first; a rank that merely takes more arithmetic to derive, as spiral order's ring-based formula does, is still an entirely independent, per-thread computation. Mistaking "complicated closed form" for "no closed form" is why spiral order is so often implemented as a sequential CPU pre-pass feeding a trivial GPU copy, when the whole computation could run as one ordinary, boundary-checked kernel launch instead.
K.7 Diagonal Traversal: Independent Within a Diagonal, Dependent Across Diagonals [FOUNDATIONAL]¶
Intuition¶
Section K.6 found that spiral order, however intricate its closed form, is still a strict total order: every position has exactly one predecessor and one successor, and nothing about it lets two positions be computed independently of each other. Diagonal order is the opposite kind of traversal, and it is the one CUDA's actual wavefront and dynamic-programming kernels are built around: every position on the same anti-diagonal (row + col constant) is independent of every other position on that same diagonal, and depends only on positions on the diagonal immediately before it. That distinction — not formula complexity — is what decides whether a computation fits one ordinary kernel launch or needs one launch per diagonal.
Background¶
A position's diagonal index is simply row + col — no ring peeling, no distance-from-edge reasoning, because diagonal membership never depended on proximity to an edge to begin with. What DOES need edge reasoning is a diagonal's own valid row range for a non-square matrix: row_min = max(0, d-cols+1), row_max = min(d, rows-1), the exact same ceiling/clamping style of bound Chapter 18.1's launch configuration already uses, just applied per diagonal instead of per whole grid.
File: 07_diagonal_traversal.cu¶
#include <cstdio>
#include <vector>
#include <utility>
// Appendix K.7 -- Section K.6 found that spiral order, despite looking
// like it has no closed form, still produces a STRICT TOTAL ORDER: every
// position has exactly one predecessor and one successor, and nothing
// about spiral order lets two positions be computed independently of one
// another. Diagonal order is the opposite kind of traversal entirely, and
// it is the one CUDA's actual wavefront/DP kernels are built around: every
// position on the SAME anti-diagonal (row + col == d) is independent of
// every other position on that diagonal, and depends only on positions on
// the diagonal immediately before it. That difference -- not a difference
// of formula complexity -- is what decides whether a traversal can be one
// kernel launch (Chapter 18's ordinary elementwise kernels) or needs one
// launch per diagonal (this section's own kernel).
// The classic zigzag print order (alternating direction each diagonal) --
// genuinely sequential, kept here only as a correctness anchor for the
// closed-form diagonal/position formulas below, the same role Section
// K.6's sequential simulation played for its own closed form.
std::vector<std::pair<int,int>> diagonal_order_zigzag(int rows, int cols) {
std::vector<std::pair<int,int>> order;
if (rows <= 0 || cols <= 0) return order;
for (int d = 0; d <= rows + cols - 2; d++) {
int row_min = (d - cols + 1 > 0) ? (d - cols + 1) : 0;
int row_max = (d < rows - 1) ? d : (rows - 1);
if (d % 2 == 0) {
for (int r = row_max; r >= row_min; r--) order.push_back({r, d - r});
} else {
for (int r = row_min; r <= row_max; r++) order.push_back({r, d - r});
}
}
return order;
}
// Closed form: which diagonal a position belongs to, and its position
// WITHIN that diagonal counting from the diagonal's own row_min -- unlike
// spiral order, this needs no ring peeling at all, because a diagonal's
// membership depends only on row + col, not on distance from an edge.
int diagonal_of(int row, int col) { return row + col; }
int row_min_of_diagonal(int d, int rows, int cols) {
(void)rows;
return (d - cols + 1 > 0) ? (d - cols + 1) : 0;
}
int row_max_of_diagonal(int d, int rows, int cols) {
(void)cols;
return (d < rows - 1) ? d : (rows - 1);
}
// ---- The CUDA side: a genuine wavefront dependency, one launch per diagonal ----
//
// Lattice path counting -- paths(r, c) is the number of right/down-only
// paths from (0,0) to (r,c): paths(r,0) = paths(0,c) = 1, and otherwise
// paths(r,c) = paths(r-1,c) + paths(r,c-1). Both of a cell's dependencies
// sit on diagonal (row+col-1) -- one diagonal earlier -- and EVERY cell on
// diagonal (row+col) depends only on that earlier diagonal, never on
// another cell on its own diagonal. That is exactly why this kernel is
// launched once per diagonal instead of once for the whole grid: nothing
// inside one launch may safely assume another thread in the SAME launch
// has already finished, but Chapter 18's ordinary elementwise kernels
// never need that assumption to begin with.
__global__ void dp_wavefront_kernel(float* grid, int rows, int cols, int d) {
int row_min = (d - cols + 1 > 0) ? (d - cols + 1) : 0;
int row_max = (d < rows - 1) ? d : (rows - 1);
int t = blockIdx.x * blockDim.x + threadIdx.x;
int row = row_min + t;
if (row > row_max) return; // boundary check: past this diagonal's last valid row
int col = d - row;
if (row == 0 || col == 0) {
grid[row * cols + col] = 1.0f;
} else {
grid[row * cols + col] = grid[(row - 1) * cols + col] + grid[row * cols + (col - 1)];
}
}
// Host-side simulation of launching dp_wavefront_kernel once per diagonal,
// d = 0 .. rows+cols-2, each launch itself simulated as a full grid/block/
// thread loop nest -- Section K.2/K.3/K.6's own verification technique,
// applied here across a SEQUENCE of launches rather than just one.
void simulate_wavefront_dp(int rows, int cols, int block_dim_x, std::vector<float>& grid) {
grid.assign((size_t)rows * cols, 0.0f);
for (int d = 0; d <= rows + cols - 2; d++) {
int row_min = (d - cols + 1 > 0) ? (d - cols + 1) : 0;
int row_max = (d < rows - 1) ? d : (rows - 1);
int diagonal_len = row_max - row_min + 1;
int blocks_x = (diagonal_len + block_dim_x - 1) / block_dim_x;
for (int bx = 0; bx < blocks_x; bx++) {
for (int tx = 0; tx < block_dim_x; tx++) {
int t = bx * block_dim_x + tx;
int row = row_min + t;
if (row > row_max) continue;
int col = d - row;
if (row == 0 || col == 0) grid[(size_t)row * cols + col] = 1.0f;
else grid[(size_t)row * cols + col] = grid[(size_t)(row - 1) * cols + col] +
grid[(size_t)row * cols + (col - 1)];
}
}
}
}
void scalar_dp_reference(int rows, int cols, std::vector<float>& grid) {
grid.assign((size_t)rows * cols, 0.0f);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
grid[(size_t)r * cols + c] = (r == 0 || c == 0) ? 1.0f
: grid[(size_t)(r - 1) * cols + c] + grid[(size_t)r * cols + (c - 1)];
}
int main() {
printf("=== Section K.7: diagonal order, and the wavefront dependency it exposes ===\n\n");
printf("Canonical 3x3 check (values 1..9), zigzag diagonal order:\n");
auto canon = diagonal_order_zigzag(3, 3);
printf(" ");
for (auto& [r, c] : canon) printf("%d ", r * 3 + c + 1);
printf("\n (expected: 1 2 4 7 5 3 6 8 9)\n\n");
printf("=== Closed-form diagonal membership, cross-checked against the zigzag simulation ===\n");
bool membership_ok = true;
int shapes[][2] = {{3,3},{4,5},{5,4},{1,6},{6,1},{7,7}};
for (auto& s : shapes) {
int rows = s[0], cols = s[1];
auto order = diagonal_order_zigzag(rows, cols);
bool shape_ok = true;
int idx = 0;
for (int d = 0; d <= rows + cols - 2; d++) {
int row_min = row_min_of_diagonal(d, rows, cols);
int row_max = row_max_of_diagonal(d, rows, cols);
for (int k = 0; k < row_max - row_min + 1; k++) {
auto [r, c] = order[idx++];
if (diagonal_of(r, c) != d) shape_ok = false;
}
}
shape_ok &= (idx == rows * cols);
membership_ok &= shape_ok;
printf(" %dx%d: every position's zigzag-visit position falls on its own row+col diagonal: %s\n",
rows, cols, shape_ok ? "true" : "false");
}
printf("\nall shapes: closed-form diagonal membership matches the zigzag order exactly: %s\n\n",
membership_ok ? "true" : "false");
printf("=== The wavefront DP kernel: genuinely simulated launch-by-launch, per diagonal ===\n");
bool all_dp_ok = true;
struct { int rows, cols, bx; } dp_cases[] = { {5, 5, 4}, {6, 9, 4}, {9, 6, 4}, {1, 8, 4}, {8, 1, 4} };
for (auto& dc : dp_cases) {
std::vector<float> wavefront_grid, reference_grid;
simulate_wavefront_dp(dc.rows, dc.cols, dc.bx, wavefront_grid);
scalar_dp_reference(dc.rows, dc.cols, reference_grid);
bool match = (wavefront_grid == reference_grid);
all_dp_ok &= match;
printf(" %dx%d, block_dim_x=%d, %d launches (one per diagonal): matches scalar DP reference: %s\n",
dc.rows, dc.cols, dc.bx, dc.rows + dc.cols - 1, match ? "true" : "false");
}
printf("\nall DP cases verified correct: %s\n", all_dp_ok ? "true" : "false");
printf("\n**[COMMON TRAP]** dp_wavefront_kernel cannot be launched once for the whole grid the way\n");
printf("Chapter 18's elementwise kernels are -- a single launch gives CUDA zero obligation to have\n");
printf("any particular thread, let alone an entire earlier diagonal's worth of threads, finish\n");
printf("before another thread reads grid[row-1][col] or grid[row][col-1]. Each case above needed\n");
printf("exactly rows+cols-1 separate launches, one per diagonal, for exactly one reason: a kernel\n");
printf("launch boundary is the only guarantee, in plain CUDA C++ without cooperative-groups grid\n");
printf("synchronization, that every thread from one launch has completed before the next launch's\n");
printf("threads begin.\n");
bool all_ok = membership_ok && all_dp_ok;
return all_ok ? 0 : 1;
}
Worked Example K.7.1 — Lattice-path counting, genuinely simulated one launch per diagonal¶
NVDIR=/usr/local/lib/python3.11/dist-packages/nvidia/cu13
nvcc -arch=sm_80 -I$NVDIR/include -L$NVDIR/lib -lcudart 07_diagonal_traversal.cu -o 07_diagonal_traversal
LD_LIBRARY_PATH=$NVDIR/lib ./07_diagonal_traversal
Genuinely compiled and run:
=== Section K.7: diagonal order, and the wavefront dependency it exposes ===
Canonical 3x3 check (values 1..9), zigzag diagonal order:
1 2 4 7 5 3 6 8 9
(expected: 1 2 4 7 5 3 6 8 9)
=== Closed-form diagonal membership, cross-checked against the zigzag simulation ===
3x3: every position's zigzag-visit position falls on its own row+col diagonal: true
4x5: every position's zigzag-visit position falls on its own row+col diagonal: true
5x4: every position's zigzag-visit position falls on its own row+col diagonal: true
1x6: every position's zigzag-visit position falls on its own row+col diagonal: true
6x1: every position's zigzag-visit position falls on its own row+col diagonal: true
7x7: every position's zigzag-visit position falls on its own row+col diagonal: true
all shapes: closed-form diagonal membership matches the zigzag order exactly: true
=== The wavefront DP kernel: genuinely simulated launch-by-launch, per diagonal ===
5x5, block_dim_x=4, 9 launches (one per diagonal): matches scalar DP reference: true
6x9, block_dim_x=4, 14 launches (one per diagonal): matches scalar DP reference: true
9x6, block_dim_x=4, 14 launches (one per diagonal): matches scalar DP reference: true
1x8, block_dim_x=4, 8 launches (one per diagonal): matches scalar DP reference: true
8x1, block_dim_x=4, 8 launches (one per diagonal): matches scalar DP reference: true
all DP cases verified correct: true
**[COMMON TRAP]** dp_wavefront_kernel cannot be launched once for the whole grid the way
Chapter 18's elementwise kernels are -- a single launch gives CUDA zero obligation to have
any particular thread, let alone an entire earlier diagonal's worth of threads, finish
before another thread reads grid[row-1][col] or grid[row][col-1]. Each case above needed
exactly rows+cols-1 separate launches, one per diagonal, for exactly one reason: a kernel
launch boundary is the only guarantee, in plain CUDA C++ without cooperative-groups grid
synchronization, that every thread from one launch has completed before the next launch's
threads begin.
The canonical 3x3 zigzag check (1 2 4 7 5 3 6 8 9) anchors the sequential simulation against the well-known result before anything else in this section trusts it. The closed-form diagonal-membership check then confirms, across six shapes including single-row and single-column matrices, that every position the zigzag order visits really does sit on the row + col diagonal the closed form assigns it to. dp_wavefront_kernel computes lattice-path counts — paths(r,c) = paths(r-1,c) + paths(r,c-1), with both dependencies sitting on the immediately preceding diagonal — and is genuinely compiled for sm_80 above; simulate_wavefront_dp genuinely executes the entire launch SEQUENCE this kernel would need, one launch per diagonal, each launch itself a full grid/block/thread simulation, verified against a plain scalar DP fill across five shapes, including both a single row and a single column.
[COMMON TRAP] dp_wavefront_kernel cannot be launched once for the whole grid the way Chapter 18's ordinary elementwise kernels are, even though it uses the exact same "compute an index, bounds-check it" template every other kernel in this book follows. A single launch gives CUDA zero obligation that any particular thread — let alone an entire earlier diagonal's worth of threads — has finished before another thread reads grid[row-1][col] or grid[row][col-1]. Section K.7's own worked example needed rows+cols-1 separate launches, one per diagonal, for exactly one reason: a kernel launch boundary is the only guarantee, in plain CUDA C++ without cooperative-groups grid synchronization, that every thread from one launch has completed before the next launch's threads begin — the price a wavefront dependency charges that an embarrassingly parallel elementwise kernel never has to pay.
Complete Runnable Code¶
01_row_major_col_major_traversal.cpp¶
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <list>
#include <unordered_map>
// Appendix K.1 -- row-major storage (Chapter 3's convention, Chapter 13's
// Matrix::get) fixes WHERE every element lives. It says nothing at all
// about the order a loop nest is allowed to VISIT those elements in. This
// file keeps storage fixed and varies only the traversal order, then
// measures the consequence with a real, deterministic cache simulator --
// not a wall-clock timing, which Getting Started's own account of this
// book's build environment already notes will "jitter somewhat from run
// to run and machine to machine" (Chapter 19's own GFLOPS figures). A
// cache-line-touch count, computed from pure index arithmetic, reproduces
// bit-for-bit on every rerun, on any machine, exactly like Appendix E.3's
// genuinely counted (not timed) global-memory reads.
struct Matrix {
float* data;
int rows, cols;
Matrix(int rows_, int cols_) : rows(rows_), cols(cols_) {
data = (float*)malloc((size_t)rows * cols * sizeof(float));
}
~Matrix() { free(data); }
float get(int row, int col) const { return data[row * cols + col]; } // row-major, Chapter 3
void set(int row, int col, float value) { data[row * cols + col] = value; }
};
// The same 2x3 buffer, but stored COLUMN-major instead -- Chapter 7.2's
// col_major offset formula (offset = col*rows + row) specialized to two
// dimensions, the same way Matrix::get specialized Chapter 6's general
// row-major formula down to two.
struct ColMajorMatrix {
float* data;
int rows, cols;
ColMajorMatrix(int rows_, int cols_) : rows(rows_), cols(cols_) {
data = (float*)malloc((size_t)rows * cols * sizeof(float));
}
~ColMajorMatrix() { free(data); }
float get(int row, int col) const { return data[col * rows + row]; } // column-major, Chapter 7.2
void set(int row, int col, float value) { data[col * rows + row] = value; }
};
// A CPU cache line is a fixed-size chunk (this book's own glossary, Appendix
// H.7, states "commonly 64 bytes, or 8 doubles" -- for 4-byte float this is
// 16 elements per line) fetched as one unit on a miss. A tiny, fully
// deterministic LRU (least-recently-used) simulator, capacity given in
// LINES rather than bytes, reproduces the real replacement decision a
// hardware cache makes without needing to touch actual memory or a clock.
struct LRUCacheSim {
int capacity_lines;
std::list<long long> lru_order; // front = most recently used
std::unordered_map<long long, std::list<long long>::iterator> present;
long long hits = 0, misses = 0;
explicit LRUCacheSim(int capacity_lines_) : capacity_lines(capacity_lines_) {}
void access(long long line_id) {
auto it = present.find(line_id);
if (it != present.end()) {
hits++;
lru_order.erase(it->second);
lru_order.push_front(line_id);
present[line_id] = lru_order.begin();
return;
}
misses++;
if ((int)lru_order.size() >= capacity_lines) {
long long victim = lru_order.back();
lru_order.pop_back();
present.erase(victim);
}
lru_order.push_front(line_id);
present[line_id] = lru_order.begin();
}
};
static const int FLOATS_PER_LINE = 16; // 64-byte line / 4-byte float, Appendix H.7's glossary
long long line_of(int row, int col, int cols) {
long long flat = (long long)row * cols + col;
return flat / FLOATS_PER_LINE;
}
// Walk a row-major-stored matrix outer-row/inner-col ("natural" order for
// row-major data) and feed every touched line, in visit order, to a fresh
// cache simulator of the given capacity.
void simulate_row_outer(int rows, int cols, int capacity_lines, long long& hits, long long& misses) {
LRUCacheSim sim(capacity_lines);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
sim.access(line_of(r, c, cols));
hits = sim.hits;
misses = sim.misses;
}
// Walk the SAME row-major-stored matrix outer-col/inner-row instead --
// storage never changes, only which index the loop nest holds fixed.
void simulate_col_outer(int rows, int cols, int capacity_lines, long long& hits, long long& misses) {
LRUCacheSim sim(capacity_lines);
for (int c = 0; c < cols; c++)
for (int r = 0; r < rows; r++)
sim.access(line_of(r, c, cols));
hits = sim.hits;
misses = sim.misses;
}
void run_case(const char* label, int rows, int cols, int capacity_lines) {
long long row_hits, row_misses, col_hits, col_misses;
simulate_row_outer(rows, cols, capacity_lines, row_hits, row_misses);
simulate_col_outer(rows, cols, capacity_lines, col_hits, col_misses);
long long total = (long long)rows * cols;
printf("%s: %dx%d matrix, %d-line cache (%d floats/line -> %d floats resident)\n",
label, rows, cols, capacity_lines, FLOATS_PER_LINE, capacity_lines * FLOATS_PER_LINE);
printf(" row-outer (matches row-major storage): %lld hits, %lld misses / %lld accesses (%.2f%% hit rate)\n",
row_hits, row_misses, total, 100.0 * row_hits / total);
printf(" col-outer (crosses row-major storage): %lld hits, %lld misses / %lld accesses (%.2f%% hit rate)\n",
col_hits, col_misses, total, 100.0 * col_hits / total);
}
int main() {
printf("=== Section K.1: row-major vs column-major TRAVERSAL of row-major-stored data ===\n\n");
run_case("Case 1", 8, 64, 4);
long long case1_row_hits, case1_row_misses, case1_col_hits, case1_col_misses;
simulate_row_outer(8, 64, 4, case1_row_hits, case1_row_misses);
simulate_col_outer(8, 64, 4, case1_col_hits, case1_col_misses);
printf("\n");
run_case("Case 2 (larger, same shape ratio)", 32, 128, 4);
printf("\n=== Correctness check: traversal order never changes the answer, only the cost ===\n");
Matrix m(8, 64);
for (int r = 0; r < m.rows; r++)
for (int c = 0; c < m.cols; c++)
m.set(r, c, (float)(r * 64 + c));
double sum_row_outer = 0.0;
for (int r = 0; r < m.rows; r++)
for (int c = 0; c < m.cols; c++)
sum_row_outer += m.get(r, c);
double sum_col_outer = 0.0;
for (int c = 0; c < m.cols; c++)
for (int r = 0; r < m.rows; r++)
sum_col_outer += m.get(r, c);
printf("sum via row-outer traversal: %.1f\n", sum_row_outer);
printf("sum via col-outer traversal: %.1f\n", sum_col_outer);
printf("identical sums, wildly different cache behavior: %s\n",
(sum_row_outer == sum_col_outer) ? "true" : "false");
printf("\n=== Column-major storage, SAME 8x64 shape: its own natural order vs crossing it ===\n");
ColMajorMatrix cm(8, 64);
for (int r = 0; r < cm.rows; r++)
for (int c = 0; c < cm.cols; c++)
cm.set(r, c, (float)(r * 64 + c));
// Reuse the identical line-touch machinery, but now the "natural" (fast)
// order for THIS storage is col-outer/row-inner, since col_major's
// offset formula makes consecutive ROWS at a fixed column contiguous.
LRUCacheSim cm_natural(4), cm_crossed(4);
for (int c = 0; c < cm.cols; c++)
for (int r = 0; r < cm.rows; r++)
cm_natural.access(((long long)c * cm.rows + r) / FLOATS_PER_LINE);
for (int r = 0; r < cm.rows; r++)
for (int c = 0; c < cm.cols; c++)
cm_crossed.access(((long long)c * cm.rows + r) / FLOATS_PER_LINE);
long long cm_total = (long long)cm.rows * cm.cols;
printf("col-major storage, col-outer traversal (matches storage): %lld hits, %lld misses / %lld (%.2f%% hit rate)\n",
cm_natural.hits, cm_natural.misses, cm_total, 100.0 * cm_natural.hits / cm_total);
printf("col-major storage, row-outer traversal (crosses storage): %lld hits, %lld misses / %lld (%.2f%% hit rate)\n",
cm_crossed.hits, cm_crossed.misses, cm_total, 100.0 * cm_crossed.hits / cm_total);
printf("matching traversal reaches the identical 93.75%% hit rate row-major's own matching case did: %s\n",
(cm_natural.hits == case1_row_hits) ? "true" : "false");
printf("\n=== Why the CROSSING numbers above are NOT a mirror image of each other ===\n");
printf("row-major 8x64, crossing stride = cols = 64 floats = %d lines (>= cache capacity 4): every access misses\n", 64 / FLOATS_PER_LINE);
printf("col-major 8x64, crossing stride = rows = 8 floats = half of one %d-float line: partial reuse survives\n", FLOATS_PER_LINE);
printf("row-major crossing: %lld hits, %lld misses (stride a whole multiple of line size)\n", case1_col_hits, case1_col_misses);
printf("col-major crossing: %lld hits, %lld misses (stride smaller than one line)\n", cm_crossed.hits, cm_crossed.misses);
printf("conclusion: the SEVERITY of crossing storage depends on the stride relative to line size,\n");
printf("not just on whether storage was crossed at all -- both are still worse than matching, but not equally worse.\n");
printf("\n=== A genuinely SYMMETRIC (square) shape DOES mirror exactly ===\n");
const int SQ = 32;
long long sq_row_hits, sq_row_misses, sq_col_hits, sq_col_misses;
simulate_row_outer(SQ, SQ, 4, sq_row_hits, sq_row_misses); // row-major storage
simulate_col_outer(SQ, SQ, 4, sq_col_hits, sq_col_misses);
LRUCacheSim sq_cm_natural(4), sq_cm_crossed(4);
for (int c = 0; c < SQ; c++)
for (int r = 0; r < SQ; r++)
sq_cm_natural.access(((long long)c * SQ + r) / FLOATS_PER_LINE);
for (int r = 0; r < SQ; r++)
for (int c = 0; c < SQ; c++)
sq_cm_crossed.access(((long long)c * SQ + r) / FLOATS_PER_LINE);
printf("row-major %dx%d, matching traversal: %lld hits, %lld misses\n", SQ, SQ, sq_row_hits, sq_row_misses);
printf("col-major %dx%d, matching traversal: %lld hits, %lld misses\n", SQ, SQ, sq_cm_natural.hits, sq_cm_natural.misses);
printf("row-major %dx%d, crossing traversal: %lld hits, %lld misses\n", SQ, SQ, sq_col_hits, sq_col_misses);
printf("col-major %dx%d, crossing traversal: %lld hits, %lld misses\n", SQ, SQ, sq_cm_crossed.hits, sq_cm_crossed.misses);
bool square_mirrors = (sq_row_hits == sq_cm_natural.hits) && (sq_col_hits == sq_cm_crossed.hits);
printf("for a square matrix, row-major and column-major produce IDENTICAL hit/miss counts both ways: %s\n",
square_mirrors ? "true" : "false");
return 0;
}
02_matrix_transpose_tiled_boundary.cu¶
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <set>
#include <vector>
// Appendix K.2 -- Appendix C.3 built the classic shared-memory tiled
// transpose and its bank-conflict fix (a [32][32] tile padded to [32][33]),
// but only for a SQUARE matrix whose width divides evenly by 32 -- neither
// `shared_transpose_naive` nor `shared_transpose_padded` there ever checks
// a bound, and both would silently read or write past the end of a matrix
// that isn't exactly 32k x 32k. Appendix E.3's own [COMMON TRAP] named
// this gap directly ("a real kernel handling a matrix size that doesn't
// divide evenly needs the same bounds-checking and padding logic") without
// building it. This file builds it: a real, non-square-capable,
// boundary-checked tiled transpose, genuinely compiled with nvcc for a
// real target architecture and, since launching it needs a device this
// sandbox does not have, checked instead by a host-side simulator that
// walks the identical grid/block/tile loop nest a real launch would --
// exactly Appendix E.3's own verification technique, applied to a kernel
// instead of a memory-traffic count.
#define TILE 32
// The naive version: no padding, so a column read of the shared tile
// still produces Appendix C.3's exact 32-way bank conflict -- but now
// ALSO correctly bounded for a rows x cols matrix that need not be
// square or a multiple of TILE. Genuinely compiled for sm_80; never
// launched in this no-GPU environment (Getting Started's own account of
// this book's build environment).
__global__ void transpose_naive_boundary_checked(float* out, const float* in, int rows, int cols) {
__shared__ float tile[TILE][TILE];
int in_row = blockIdx.y * TILE + threadIdx.y;
int in_col = blockIdx.x * TILE + threadIdx.x;
if (in_row < rows && in_col < cols) {
tile[threadIdx.y][threadIdx.x] = in[in_row * cols + in_col];
}
__syncthreads();
// out is cols x rows (the transposed shape): swap which block axis
// supplies the row vs column half of the output coordinate, exactly
// the way Appendix C.3's square version swaps x and y in one line --
// generalized here into two full, independently-bounded coordinates.
int out_row = blockIdx.x * TILE + threadIdx.y;
int out_col = blockIdx.y * TILE + threadIdx.x;
if (out_row < cols && out_col < rows) {
out[out_row * rows + out_col] = tile[threadIdx.x][threadIdx.y]; // column read -- the conflict
}
}
// The fix: Appendix C.3's identical one-float row padding, plus the same
// boundary checks as the naive version above. Padding and bounds-checking
// are orthogonal concerns -- one is about shared-memory bank arithmetic,
// the other about which threads own real data -- and this kernel is the
// genuine, compiled proof that combining them needs no new idea, just
// both fixes applied at once.
__global__ void transpose_padded_boundary_checked(float* out, const float* in, int rows, int cols) {
__shared__ float tile[TILE][TILE + 1];
int in_row = blockIdx.y * TILE + threadIdx.y;
int in_col = blockIdx.x * TILE + threadIdx.x;
if (in_row < rows && in_col < cols) {
tile[threadIdx.y][threadIdx.x] = in[in_row * cols + in_col];
}
__syncthreads();
int out_row = blockIdx.x * TILE + threadIdx.y;
int out_col = blockIdx.y * TILE + threadIdx.x;
if (out_row < cols && out_col < rows) {
out[out_row * rows + out_col] = tile[threadIdx.x][threadIdx.y];
}
}
// ---- Everything below is host C++, genuinely compiled AND genuinely run ----
// Appendix C.3's exact bank-index arithmetic, reused verbatim on this
// file's own row strides to confirm the fix still holds unchanged.
int compute_distinct_banks(int row_stride_words) {
std::set<int> banks;
for (int t = 0; t < 32; t++) {
banks.insert((t * row_stride_words) % 32);
}
return (int)banks.size();
}
// A host-side simulator of the exact grid/block/thread/two-phase loop nest
// both kernels above run on a device -- Appendix E.3's own verification
// technique (simulate the launch, don't just trust the formula), applied
// here to a transpose instead of a read count. tile_width is a RUNTIME
// parameter (like E.3's own TILE_WIDTH), independent of the kernels'
// compile-time TILE, specifically so a small, hand-checkable tile_width
// can exercise the boundary path without touching the real kernels' macro.
void simulate_tiled_transpose(const std::vector<float>& in, int rows, int cols,
int tile_width, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f); // out is cols x rows; same element count
int blocks_x = (cols + tile_width - 1) / tile_width; // ceiling division, Chapter 18.1
int blocks_y = (rows + tile_width - 1) / tile_width;
std::vector<float> tile((size_t)tile_width * tile_width);
for (int by = 0; by < blocks_y; by++) {
for (int bx = 0; bx < blocks_x; bx++) {
// Phase 1: cooperative load into the simulated shared tile.
for (int ty = 0; ty < tile_width; ty++) {
for (int tx = 0; tx < tile_width; tx++) {
int in_row = by * tile_width + ty;
int in_col = bx * tile_width + tx;
if (in_row < rows && in_col < cols) {
tile[ty * tile_width + tx] = in[(size_t)in_row * cols + in_col];
}
}
}
// __syncthreads(): on real hardware this barrier is why phase 2
// may not begin until every thread in the block has finished
// phase 1. A sequential host simulation is already ordered, so
// the barrier is implied by nothing more than "this loop already
// finished" -- the same reasoning Appendix E.3 relies on when it
// simulates a launch without a real device to synchronize on.
// Phase 2: write the transposed tile out, boundary-checked
// against the OUTPUT shape (cols x rows), swapped block axes.
for (int ty = 0; ty < tile_width; ty++) {
for (int tx = 0; tx < tile_width; tx++) {
int out_row = bx * tile_width + ty;
int out_col = by * tile_width + tx;
if (out_row < cols && out_col < rows) {
out[(size_t)out_row * rows + out_col] = tile[tx * tile_width + ty];
}
}
}
}
}
}
void scalar_transpose_reference(const std::vector<float>& in, int rows, int cols, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
out[(size_t)c * rows + r] = in[(size_t)r * cols + c];
}
bool run_case(const char* label, int rows, int cols, int tile_width) {
std::vector<float> in((size_t)rows * cols);
for (size_t i = 0; i < in.size(); i++) in[i] = (float)(i + 1);
std::vector<float> tiled_out, reference_out;
simulate_tiled_transpose(in, rows, cols, tile_width, tiled_out);
scalar_transpose_reference(in, rows, cols, reference_out);
bool match = (tiled_out == reference_out);
int blocks_x = (cols + tile_width - 1) / tile_width;
int blocks_y = (rows + tile_width - 1) / tile_width;
bool boundary_case = (rows % tile_width != 0) || (cols % tile_width != 0);
printf("%s: %dx%d matrix, tile_width=%d, grid=%dx%d blocks (boundary tiles: %s)\n",
label, rows, cols, tile_width, blocks_x, blocks_y, boundary_case ? "yes" : "no");
printf(" tiled+boundary-checked simulation matches scalar reference, element for element: %s\n",
match ? "true" : "false");
return match;
}
int main() {
printf("=== Section K.2: tiled transpose with boundary checks, genuinely simulated and verified ===\n\n");
bool all_ok = true;
all_ok &= run_case("Case 1 (rows and cols both divide tile_width evenly)", 8, 8, 4);
all_ok &= run_case("Case 2 (neither dimension divides tile_width evenly)", 7, 5, 4);
all_ok &= run_case("Case 3 (larger, still non-divisible, non-square)", 10, 13, 4);
all_ok &= run_case("Case 4 (a dimension smaller than one whole tile)", 3, 11, 4);
printf("\nall %d cases verified correct: %s\n", 4, all_ok ? "true" : "false");
printf("\n=== Bank-conflict arithmetic on THIS file's tile shapes, Appendix C.3's exact method ===\n");
int naive_banks = compute_distinct_banks(TILE); // tile[32][32]: row stride = 32 floats
int padded_banks = compute_distinct_banks(TILE + 1); // tile[32][33]: row stride = 33 floats
printf("transpose_naive_boundary_checked's tile[%d][%d], column access: %d distinct bank(s) (1 = full 32-way conflict)\n",
TILE, TILE, naive_banks);
printf("transpose_padded_boundary_checked's tile[%d][%d], column access: %d distinct bank(s) (32 = fully resolved)\n",
TILE, TILE + 1, padded_banks);
printf("adding boundary checks changes NOTHING about the bank math -- it depends only on threadIdx and\n");
printf("row stride, never on whether a given thread's data happened to be in-bounds: %s\n",
(naive_banks == 1 && padded_banks == 32) ? "confirmed, matches Appendix C.3 exactly" : "MISMATCH");
return all_ok ? 0 : 1;
}
03_grid_mapping_launch_config_2d.cpp¶
#include <cstdio>
#include <vector>
// Appendix K.3 -- Chapter 18.1 built compute_launch_config for a flat,
// one-dimensional problem: ceiling-divide size by threads_per_block, then
// guard the leftover threads with if (idx < size). A matrix needs the
// identical decision made TWICE, once per axis, and it needs one more
// decision Chapter 18.1 never had to make at all: which matrix axis
// becomes threadIdx.x and which becomes threadIdx.y. This file makes both
// decisions explicit, extends Chapter 18.1's launch-config struct to two
// dimensions, and then PROVES -- by genuinely running every (blockIdx,
// threadIdx) pair a real launch would produce, not by arguing it in prose
// -- that the resulting grid covers every element of an M x N matrix
// exactly once, no matter how M and N divide against the chosen block
// shape.
struct LaunchConfig2D {
int blocks_x, blocks_y; // blocks_x covers cols, blocks_y covers rows
int block_dim_x, block_dim_y;
long long total_threads_x, total_threads_y;
long long wasted_x, wasted_y;
};
// Direct two-axis extension of Chapter 18.1's compute_launch_config --
// identical ceiling-division formula, applied once per axis. Convention
// (matching every 2D kernel this book has launched since Chapter 18.3 and
// Appendix C.3/E.3): x indexes COLUMNS, y indexes ROWS.
LaunchConfig2D compute_launch_config_2d(int rows, int cols, int block_dim_x, int block_dim_y) {
LaunchConfig2D cfg;
cfg.block_dim_x = block_dim_x;
cfg.block_dim_y = block_dim_y;
cfg.blocks_x = (cols + block_dim_x - 1) / block_dim_x;
cfg.blocks_y = (rows + block_dim_y - 1) / block_dim_y;
cfg.total_threads_x = (long long)cfg.blocks_x * block_dim_x;
cfg.total_threads_y = (long long)cfg.blocks_y * block_dim_y;
cfg.wasted_x = cfg.total_threads_x - cols;
cfg.wasted_y = cfg.total_threads_y - rows;
return cfg;
}
void print_config(const char* label, int rows, int cols, const LaunchConfig2D& cfg) {
printf("%s: %dx%d matrix, block=%dx%d\n", label, rows, cols, cfg.block_dim_x, cfg.block_dim_y);
printf(" grid = %d x %d blocks, total launched threads = %lld x %lld\n",
cfg.blocks_x, cfg.blocks_y, cfg.total_threads_x, cfg.total_threads_y);
long long total_launched = cfg.total_threads_x * cfg.total_threads_y;
long long total_needed = (long long)rows * cols;
printf(" wasted (out-of-bounds, guarded) threads: %lld along x, %lld along y, %lld total (%lld launched - %lld needed)\n",
cfg.wasted_x, cfg.wasted_y, total_launched - total_needed, total_launched, total_needed);
}
// Genuinely run every (blockIdx.x, blockIdx.y, threadIdx.x, threadIdx.y)
// combination a real <<<grid, block>>> launch would produce, apply
// Chapter 18.1's own bounds-check pattern, and count how many times each
// logical (row, col) position gets touched. Exhaustive, deterministic,
// and -- unlike a device launch -- fully executable in this no-GPU
// environment, exactly the way Appendix H.3's for_each_index walk and
// Appendix E.3's simulate_tiled both stand in for hardware this sandbox
// does not have.
bool verify_exact_coverage(int rows, int cols, const LaunchConfig2D& cfg) {
std::vector<int> touch_count((size_t)rows * cols, 0);
for (int by = 0; by < cfg.blocks_y; by++) {
for (int bx = 0; bx < cfg.blocks_x; bx++) {
for (int ty = 0; ty < cfg.block_dim_y; ty++) {
for (int tx = 0; tx < cfg.block_dim_x; tx++) {
int out_row = by * cfg.block_dim_y + ty; // blockIdx.y * blockDim.y + threadIdx.y
int out_col = bx * cfg.block_dim_x + tx; // blockIdx.x * blockDim.x + threadIdx.x
if (out_row < rows && out_col < cols) {
touch_count[(size_t)out_row * cols + out_col]++;
}
}
}
}
}
long long zero = 0, exactly_one = 0, more_than_one = 0;
for (int v : touch_count) {
if (v == 0) zero++;
else if (v == 1) exactly_one++;
else more_than_one++;
}
printf(" coverage check: %lld positions touched 0 times, %lld touched exactly once, %lld touched more than once\n",
zero, exactly_one, more_than_one);
return zero == 0 && more_than_one == 0 && exactly_one == (long long)rows * cols;
}
int main() {
printf("=== Section K.3: 2D ceiling-division launch configs, exhaustively verified ===\n\n");
struct { int rows, cols, bx, by; const char* label; } cases[] = {
{32, 32, 16, 16, "Case 1 (evenly divisible both axes)"},
{37, 61, 16, 16, "Case 2 (neither axis divides evenly)"},
{5, 200, 16, 16, "Case 3 (one axis much smaller than block)"},
{1000, 1000, 32, 8, "Case 4 (non-square block shape)"},
};
bool all_ok = true;
for (auto& c : cases) {
LaunchConfig2D cfg = compute_launch_config_2d(c.rows, c.cols, c.bx, c.by);
print_config(c.label, c.rows, c.cols, cfg);
all_ok &= verify_exact_coverage(c.rows, c.cols, cfg);
printf("\n");
}
printf("every case covers its matrix exactly once, wasting only the counted leftover threads: %s\n\n",
all_ok ? "true" : "false");
printf("=== Why x maps to columns: consecutive threadIdx.x must produce consecutive addresses ===\n");
int cols = 100, block_dim_x = 32, block_dim_y = 8;
int fixed_block_row = 2, fixed_block_col = 3, fixed_ty = 1;
printf("fixed blockIdx=(%d,%d), threadIdx.y=%d; varying threadIdx.x across one warp (0..7 shown):\n",
fixed_block_col, fixed_block_row, fixed_ty);
printf("CORRECT convention (out_col uses threadIdx.x, out_row uses threadIdx.y):\n");
for (int tx = 0; tx < 8; tx++) {
int out_row = fixed_block_row * block_dim_y + fixed_ty;
int out_col = fixed_block_col * block_dim_x + tx;
long long addr = (long long)out_row * cols + out_col;
printf(" threadIdx.x=%d -> (row=%d, col=%d) -> address %lld\n", tx, out_row, out_col, addr);
}
printf(" addresses for consecutive threadIdx.x are consecutive integers (stride 1): coalesced\n");
printf("\n**[COMMON TRAP]** swapping which axis threadIdx.x indexes still touches every element,\n");
printf("but destroys coalescing -- genuinely computed, not merely asserted:\n");
for (int tx = 0; tx < 8; tx++) {
// Swapped convention: out_row now uses threadIdx.x, out_col uses threadIdx.y.
int out_row = fixed_block_col * block_dim_x + tx;
int out_col = fixed_block_row * block_dim_y + fixed_ty;
long long addr = (long long)out_row * cols + out_col;
printf(" threadIdx.x=%d -> (row=%d, col=%d) -> address %lld\n", tx, out_row, out_col, addr);
}
printf(" addresses for consecutive threadIdx.x now differ by cols=%d, not 1: NOT coalesced\n", cols);
printf(" every element is still visited exactly once (correctness survives), but each of the 8 lanes\n");
printf(" above would issue its own separate memory transaction instead of sharing one\n");
return all_ok ? 0 : 1;
}
04_boundary_check_2d_add_kernel.cu¶
#include <cstdio>
#include <cstdlib>
#include <vector>
// Appendix K.4 -- Chapter 13.1's [COMMON TRAP] showed a shape-check that
// simply doesn't exist (a.cols vs b.rows). This section shows the CUDA
// analogue: a boundary check that exists but is INCOMPLETE -- checking
// only one of a 2D kernel's two axes, a mistake easy to make by pattern-
// matching Chapter 18.1's 1D "if (idx < size)" onto two dimensions and
// forgetting one half. Both kernels below fill each matrix position with
// a distinctive, position-derived marker (row*1000+col) specifically so a
// corrupted position is obviously wrong rather than coincidentally right --
// an earlier draft of this exact section used element-wise addition
// instead, and discovered a genuinely different, subtler fact: because
// out[row*cols+col] = a[row*cols+col] + b[row*cols+col] reads and writes
// through the IDENTICAL flat-index formula for all three arrays, an
// overflowing (row, col) aliases to the exact same flat slot as the
// legitimate (row+1, 0) and reads the exact same legitimate data -- so
// that particular bug turns out to be silently HARMLESS for every row but
// the last. That is a real, worth-knowing fact about element-wise ops
// specifically (Section K.4's own Worked Example below proves it), not a
// property of boundary-check bugs in general -- this section's fill
// kernel is chosen precisely because its written value depends on the
// (possibly wrong) row/col identity itself, not on a re-aliased read, so
// a genuine wrong answer is what a missing check produces here.
//
// Both kernels are genuinely compiled with nvcc for a real architecture;
// neither is launched (this sandbox has no device), so both are checked
// the way K.2 and K.3 already checked their own kernels: a host-side
// simulator walking the identical grid/block/thread loop nest, run for
// real.
__global__ void matrix_fill_boundary_checked(float* out, int rows, int cols) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < rows && col < cols) {
out[row * cols + col] = (float)(row * 1000 + col);
}
}
// The bug: only the ROW half of Chapter 18.1's bounds-check pattern
// survived. For any cols that doesn't divide blockDim.x evenly, threads
// past the true last column are never turned away -- they still compute
// row*cols+col using their OWN (out-of-range) col, so the marker they
// write is wrong for whatever flat slot that overshoot actually lands on.
__global__ void matrix_fill_buggy_row_only_checked(float* out, int rows, int cols) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < rows) { // BUG: col is never checked against cols at all
out[row * cols + col] = (float)(row * 1000 + col);
}
}
// ---- Host C++ below: genuinely compiled AND genuinely run ----
void fill_reference(int rows, int cols, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
out[(size_t)r * cols + c] = (float)(r * 1000 + c);
}
// Simulates matrix_fill_boundary_checked's exact grid/block/thread loop nest.
void simulate_correct_fill(int rows, int cols, int block_dim_x, int block_dim_y, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
int blocks_x = (cols + block_dim_x - 1) / block_dim_x;
int blocks_y = (rows + block_dim_y - 1) / block_dim_y;
for (int by = 0; by < blocks_y; by++)
for (int bx = 0; bx < blocks_x; bx++)
for (int ty = 0; ty < block_dim_y; ty++)
for (int tx = 0; tx < block_dim_x; tx++) {
int row = by * block_dim_y + ty;
int col = bx * block_dim_x + tx;
if (row < rows && col < cols) out[(size_t)row * cols + col] = (float)(row * 1000 + col);
}
}
// Simulates matrix_fill_buggy_row_only_checked's exact (buggy) loop nest.
// `out` is sized (rows+1) x cols -- one whole extra PADDING row, filled
// with a known sentinel, rather than exactly rows x cols. This is Chapter
// 13.1's own sentinel technique (Worked Example 13.1.4): it turns what
// would otherwise be a genuine, undefined-behavior heap overflow on the
// very last real row into a real, safe, deterministic write onto memory
// this file legitimately owns, so the corruption is observed directly
// rather than merely asserted to be dangerous.
//
// process_blocks_in_reverse controls only the ORDER this host simulation
// visits blockIdx.y in -- CUDA itself makes zero guarantees about block
// execution order or concurrency (this book has relied on that fact since
// Chapter 18's own launch-configuration discussion), so a real GPU could
// legally run these blocks in either order, both orders interleaved, or
// fully in parallel. This parameter exists specifically to demonstrate
// that the bug's VISIBLE damage depends entirely on an order nothing
// promises you.
void simulate_buggy_fill(int rows, int cols, int block_dim_x, int block_dim_y,
bool process_blocks_in_reverse, std::vector<float>& out) {
int blocks_x = (cols + block_dim_x - 1) / block_dim_x;
int blocks_y = (rows + block_dim_y - 1) / block_dim_y;
for (int by_i = 0; by_i < blocks_y; by_i++) {
int by = process_blocks_in_reverse ? (blocks_y - 1 - by_i) : by_i;
for (int bx = 0; bx < blocks_x; bx++)
for (int ty = 0; ty < block_dim_y; ty++)
for (int tx = 0; tx < block_dim_x; tx++) {
int row = by * block_dim_y + ty;
int col = bx * block_dim_x + tx;
if (row < rows) { // the bug, reproduced exactly
out[(size_t)row * cols + col] = (float)(row * 1000 + col);
}
}
}
}
int main() {
printf("=== Section K.4: a boundary check that guards only ONE of two axes ===\n\n");
int rows = 4, cols = 10, block_dim_x = 4, block_dim_y = 1;
printf("matrix: %dx%d, block=%dx%d (cols=%d is NOT a multiple of block_dim_x=%d)\n\n",
rows, cols, block_dim_x, block_dim_y, cols, block_dim_x);
std::vector<float> reference;
fill_reference(rows, cols, reference);
printf("=== The CORRECT kernel, fully boundary-checked ===\n");
std::vector<float> correct_out;
simulate_correct_fill(rows, cols, block_dim_x, block_dim_y, correct_out);
bool correct_matches = (correct_out == reference);
printf("matches the row*1000+col reference, element for element: %s\n\n", correct_matches ? "true" : "false");
printf("=== Worked Example: why the SAME bug on ELEMENT-WISE ADDITION would be harmless ===\n");
printf("out[row*cols+col] = a[row*cols+col] + b[row*cols+col] reads a and b through the SAME\n");
printf("flat formula it writes out through. An overflowing (row=0, col=cols=%d) computes flat index\n", cols);
printf("%d, which is ALSO exactly (row=1, col=0)'s own flat index in every one of a, b, and out --\n", cols);
printf("so the value written is a(1,0)+b(1,0), the mathematically correct value for that slot,\n");
printf("computed one iteration early -- the SAME identity every interior overflow relies on.\n");
printf("This fill kernel's written value depends on the CALLING thread's own (row, col), not on a\n");
printf("re-aliased read, so it has no such accidental self-correction -- the corruption below is real.\n\n");
printf("=== The BUGGY fill kernel: row checked, column never checked ===\n");
const float SENTINEL = -999.0f;
int padded_rows = rows + 1; // one extra real row of storage, filled with a sentinel
int overflow_cols_per_row = (block_dim_x - (cols % block_dim_x)) % block_dim_x;
printf("overflowing columns per row (col >= cols but still written): %d\n\n", overflow_cols_per_row);
std::vector<float> buggy_forward((size_t)padded_rows * cols, SENTINEL);
std::vector<float> buggy_reverse((size_t)padded_rows * cols, SENTINEL);
simulate_buggy_fill(rows, cols, block_dim_x, block_dim_y, false, buggy_forward);
simulate_buggy_fill(rows, cols, block_dim_x, block_dim_y, true, buggy_reverse);
auto count_corrupted_rows = [&](const std::vector<float>& out) {
int corrupted = 0;
for (int r = 1; r < rows; r++) {
bool bad = (out[(size_t)r * cols + 0] != reference[(size_t)r * cols + 0]) ||
(out[(size_t)r * cols + 1] != reference[(size_t)r * cols + 1]);
if (bad) corrupted++;
}
return corrupted;
};
auto count_sentinel_overwritten = [&](const std::vector<float>& out) {
int n = 0;
for (int c = 0; c < cols; c++) if (out[(size_t)rows * cols + c] != SENTINEL) n++;
return n;
};
int fwd_corrupted_rows = count_corrupted_rows(buggy_forward);
int rev_corrupted_rows = count_corrupted_rows(buggy_reverse);
int fwd_sentinel = count_sentinel_overwritten(buggy_forward);
int rev_sentinel = count_sentinel_overwritten(buggy_reverse);
printf("processing blockIdx.y forward (0, 1, 2, ...): %d of %d interior rows visibly corrupted\n",
fwd_corrupted_rows, rows - 1);
printf("processing blockIdx.y in reverse (..., 2, 1, 0): %d of %d interior rows visibly corrupted\n",
rev_corrupted_rows, rows - 1);
printf("padding-row sentinel cells overwritten, forward order: %d\n", fwd_sentinel);
printf("padding-row sentinel cells overwritten, reverse order: %d\n\n", rev_sentinel);
printf("row-by-row detail, reverse order (the one that DOES show interior corruption):\n");
for (int r = 0; r < rows; r++) {
float got0 = buggy_reverse[(size_t)r * cols + 0];
float want0 = reference[(size_t)r * cols + 0];
printf(" row %d, col 0: wrote %.0f, correct is %.0f -- %s\n",
r, got0, want0, (got0 == want0) ? "correct" : "WRONG (overwritten by the row above's overflow)");
}
printf("\n**[COMMON TRAP]** in FORWARD order, each row's own correct write happens AFTER the row\n");
printf("above's overflow lands on it, so the correct value silently overwrites the bad one -- every\n");
printf("interior row looks fine by accident. In REVERSE order, each row's own correct write happens\n");
printf("BEFORE the row above's overflow arrives, so the bad value survives and every interior row is\n");
printf("visibly wrong. CUDA guarantees neither order (nor that blocks run one at a time at all) --\n");
printf("matrix_fill_buggy_row_only_checked compiles cleanly and produces a genuine, scheduling-\n");
printf("dependent data race, not a fixed wrong answer. The one piece of damage every order agrees on\n");
printf("is the final row's overflow, because nothing ever writes to that memory a second time to\n");
printf("cover it back up -- exactly the %d sentinel cells both orders show above.\n", overflow_cols_per_row);
bool orders_genuinely_differ = (fwd_corrupted_rows != rev_corrupted_rows);
bool forward_self_heals = (fwd_corrupted_rows == 0);
bool reverse_fully_corrupts = (rev_corrupted_rows == rows - 1);
bool last_row_overflow_always_visible = (fwd_sentinel == overflow_cols_per_row) && (rev_sentinel == overflow_cols_per_row);
bool all_ok = correct_matches && orders_genuinely_differ && forward_self_heals &&
reverse_fully_corrupts && last_row_overflow_always_visible;
printf("\nall counted quantities are internally consistent: %s\n", all_ok ? "true" : "false");
return all_ok ? 0 : 1;
}
05_row_major_col_major_kernel_interop.cpp¶
#include <cstdio>
#include <cstdlib>
#include <vector>
// Appendix K.5 -- Chapter 7.2 named the one place this book's own row-major
// convention meets a genuinely column-major world head-on: "cuBLAS, the
// CUDA math library every serious GEMM in this book eventually calls into,
// is column-major by convention." K.3 established that x should index the
// FASTEST-VARYING (contiguous) axis for coalescing. For row-major data
// that axis is the column; for column-major data it's the ROW -- the exact
// same coalescing argument, pointed at the opposite axis, because the
// underlying rule was never "x means column," it was "x means whichever
// axis is actually contiguous in memory." This section makes that
// reversal concrete, builds a genuine row-major <-> column-major converter
// for handing data to (or receiving it from) a column-major API, and
// closes with the naive-reinterpretation trap Chapter 13.3 already proved
// is a silent transpose, now in the specific context this book will
// actually hit it: a cuBLAS call.
struct RowMajorMatrix {
std::vector<float> data;
int rows, cols;
RowMajorMatrix(int rows_, int cols_) : data((size_t)rows_ * cols_, 0.0f), rows(rows_), cols(cols_) {}
float get(int r, int c) const { return data[(size_t)r * cols + c]; } // row-major, Chapter 3
void set(int r, int c, float v) { data[(size_t)r * cols + c] = v; }
};
struct ColMajorMatrix {
std::vector<float> data;
int rows, cols;
ColMajorMatrix(int rows_, int cols_) : data((size_t)rows_ * cols_, 0.0f), rows(rows_), cols(cols_) {}
float get(int r, int c) const { return data[(size_t)c * rows + r]; } // column-major, Chapter 7.2
void set(int r, int c, float v) { data[(size_t)c * rows + r] = v; }
};
int main() {
printf("=== Section K.5: which axis is contiguous REVERSES between the two layouts ===\n\n");
int rows = 4, cols = 6;
printf("for a %dx%d matrix:\n", rows, cols);
printf(" row-major: stepping one COLUMN moves 1 flat slot; stepping one ROW moves %d slots\n", cols);
printf(" col-major: stepping one ROW moves 1 flat slot; stepping one COLUMN moves %d slots\n\n", rows);
printf("K.3's coalescing rule ('x should index whichever axis is contiguous') therefore assigns:\n");
printf(" row-major kernel: threadIdx.x -> column, threadIdx.y -> row (K.3's own convention)\n");
printf(" column-major kernel: threadIdx.x -> row, threadIdx.y -> column (the exact reversal)\n\n");
// Genuinely confirm the reversal with real address arithmetic, exactly
// the way K.3 confirmed its own row-major convention.
RowMajorMatrix rm(rows, cols);
ColMajorMatrix cm(rows, cols);
printf("row-major: consecutive columns at fixed row, addresses: ");
for (int c = 0; c < cols; c++) printf("%d ", 2 * cols + c); // row=2 fixed
printf(" (stride 1)\n");
printf("col-major: consecutive rows at fixed column, addresses: ");
for (int r = 0; r < rows; r++) printf("%d ", 2 * rows + r); // col=2 fixed
printf(" (stride 1)\n");
printf("row-major: consecutive rows at fixed column, addresses: ");
for (int r = 0; r < rows; r++) printf("%d ", r * cols + 2); // col=2 fixed
printf(" (stride %d)\n", cols);
printf("col-major: consecutive columns at fixed row, addresses: ");
for (int c = 0; c < cols; c++) printf("%d ", c * rows + 2); // row=2 fixed
printf(" (stride %d)\n\n", rows);
// Fill a row-major matrix with a distinctive per-position marker --
// this book's own default Tensor layout since Chapter 6.
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
rm.set(r, c, (float)(r * 100 + c));
printf("=== Converting row-major data (this book's Tensor default) for a column-major API ===\n");
// Two independently-written conversions, cross-checked against each
// other -- this book's own convention (e.g. Worked Example 13.1.1's
// scalar/SIMD cross-check) for trusting a result.
ColMajorMatrix converted_a(rows, cols);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
converted_a.set(r, c, rm.get(r, c)); // logical (r,c) copy, layout handled by set()/get()
ColMajorMatrix converted_b(rows, cols);
for (size_t flat_rm = 0; flat_rm < rm.data.size(); flat_rm++) {
int r = (int)(flat_rm / cols);
int c = (int)(flat_rm % cols);
converted_b.data[(size_t)c * rows + r] = rm.data[flat_rm]; // direct flat-index conversion
}
bool converters_agree = (converted_a.data == converted_b.data);
bool logical_values_preserved = true;
for (int r = 0; r < rows && logical_values_preserved; r++)
for (int c = 0; c < cols; c++)
if (converted_a.get(r, c) != rm.get(r, c)) { logical_values_preserved = false; break; }
printf("two independent row-major -> column-major converters agree, byte for byte: %s\n",
converters_agree ? "true" : "false");
printf("every logical (row, col) value is preserved across the layout change: %s\n\n",
logical_values_preserved ? "true" : "false");
printf("=== Round trip: column-major -> row-major -> column-major returns the original exactly ===\n");
RowMajorMatrix back_to_row(rows, cols);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
back_to_row.set(r, c, converted_a.get(r, c));
bool round_trip_matches = (back_to_row.data == rm.data);
printf("row-major -> column-major -> row-major reproduces the original flat buffer exactly: %s\n\n",
round_trip_matches ? "true" : "false");
printf("=== [COMMON TRAP] reinterpreting a column-major buffer AS row-major without converting ===\n");
printf("This is Chapter 13.3's own naive-reshape trap, hit again here for a different reason: a\n");
printf("cuBLAS call handed row-major bytes without conversion (or a caller that reads cuBLAS's real\n");
printf("column-major OUTPUT as if it were row-major) doesn't crash and doesn't raise an error -- it\n");
printf("silently reads the TRANSPOSE of the intended matrix, because reinterpreting a flat buffer's\n");
printf("layout is exactly a relabeling of which index is contiguous, and swapping which index is\n");
printf("contiguous is precisely what transpose means.\n\n");
RowMajorMatrix naive_reinterpretation(cols, rows); // note: dims swapped, since that's what the bug does
for (size_t i = 0; i < naive_reinterpretation.data.size(); i++)
naive_reinterpretation.data[i] = converted_a.data[i]; // same bytes, wrong layout tag
bool matches_transpose = true;
for (int r = 0; r < rows && matches_transpose; r++)
for (int c = 0; c < cols; c++)
if (naive_reinterpretation.get(c, r) != rm.get(r, c)) { matches_transpose = false; break; }
printf("naive_reinterpretation.get(c, r) genuinely equals rm.get(r, c) for every (r, c): %s\n",
matches_transpose ? "true" : "false");
printf("i.e. reading the column-major buffer's bytes as row-major recovers exactly A-transpose,\n");
printf("never A itself, silently and without any error -- the identical failure mode Chapter 13.3\n");
printf("already proved for a naive reshape of a transposed buffer, now shown for a layout mismatch.\n");
bool all_ok = converters_agree && logical_values_preserved && round_trip_matches && matches_transpose;
printf("\nall counted/checked quantities are internally consistent: %s\n", all_ok ? "true" : "false");
return all_ok ? 0 : 1;
}
06_spiral_traversal.cu¶
#include <cstdio>
#include <vector>
#include <utility>
// Appendix K.6 -- every traversal order so far in this appendix (row-major,
// column-major, the tiled transpose's block order) has had a trivial,
// already-linear index formula: row*cols+col or its column-major mirror.
// Spiral order does not: which position a loop visits 5th depends on
// having already visited four others first, peeling the matrix's outer
// ring before the next ring in is even touched. That looks, at first,
// like a traversal with no closed form at all -- and therefore no way to
// hand one (row, col) to one CUDA thread and ask "what is your rank in
// this order" without that thread first simulating the whole peel. This
// section derives the closed form anyway, cross-checks it exhaustively
// against the straightforward sequential simulation, and then uses it to
// write a real CUDA kernel where every thread computes its own spiral
// rank independently -- no thread ever needs to know what any other
// thread's rank is.
// The textbook approach: shrink four boundaries inward, appending
// positions as each side of the current ring is walked. Genuinely
// sequential -- each appended position depends on the loop having
// already reached this exact point in this exact order.
std::vector<std::pair<int,int>> spiral_order_sequential(int rows, int cols) {
std::vector<std::pair<int,int>> order;
if (rows <= 0 || cols <= 0) return order;
int top = 0, bottom = rows - 1, left = 0, right = cols - 1;
while (top <= bottom && left <= right) {
for (int c = left; c <= right; c++) order.push_back({top, c});
top++;
for (int r = top; r <= bottom; r++) order.push_back({r, right});
right--;
if (top <= bottom) {
for (int c = right; c >= left; c--) order.push_back({bottom, c});
bottom--;
}
if (left <= right) {
for (int r = bottom; r >= top; r--) order.push_back({r, left});
left++;
}
}
return order;
}
// The closed form. Every position's RING is how far it sits from the
// nearest edge: k = min(row, col, rows-1-row, cols-1-col). Peeling k
// rings removes exactly rows*cols - (rows-2k)*(cols-2k) elements, so
// that quantity alone is every ring-k position's shared starting rank --
// computed in O(1), no simulation of the outer rings required. Within
// ring k, a position's offset depends only on which of the ring's four
// sides it sits on, each a direct subtraction against that side's own
// fixed coordinate.
long long spiral_rank_closed_form(int row, int col, int rows, int cols) {
int k = row;
if (col < k) k = col;
if (rows - 1 - row < k) k = rows - 1 - row;
if (cols - 1 - col < k) k = cols - 1 - col;
int top = k, bottom = rows - 1 - k, left = k, right = cols - 1 - k;
long long ring_start = (long long)rows * cols - (long long)(rows - 2 * k) * (cols - 2 * k);
int width = right - left + 1;
int height = bottom - top + 1;
if (height == 1) { // single remaining row: pure left-to-right
return ring_start + (col - left);
}
if (width == 1) { // single remaining column: pure top-to-bottom
return ring_start + (row - top);
}
if (row == top) { // top side
return ring_start + (col - left);
} else if (col == right) { // right side
return ring_start + (width - 1) + (row - top);
} else if (row == bottom) { // bottom side
return ring_start + (width - 1) + (height - 1) + (right - col);
} else { // left side
return ring_start + 2 * (width - 1) + (height - 1) + (bottom - row);
}
}
bool cross_check_shape(int rows, int cols) {
auto sequential = spiral_order_sequential(rows, cols);
std::vector<long long> closed_form_rank((size_t)rows * cols);
bool all_in_range = true;
for (int r = 0; r < rows; r++) {
for (int c = 0; c < cols; c++) {
long long rank = spiral_rank_closed_form(r, c, rows, cols);
if (rank < 0 || rank >= (long long)rows * cols) { all_in_range = false; continue; }
closed_form_rank[(size_t)r * cols + c] = rank;
}
}
bool matches = all_in_range;
if (all_in_range) {
for (size_t rank = 0; rank < sequential.size() && matches; rank++) {
auto [r, c] = sequential[rank];
if (closed_form_rank[(size_t)r * cols + c] != (long long)rank) matches = false;
}
}
printf(" %dx%d: closed form matches sequential simulation for all %d positions: %s\n",
rows, cols, rows * cols, matches ? "true" : "false");
return matches;
}
// ---- The CUDA kernel: every thread computes its OWN rank, independently ----
__global__ void spiral_scatter_kernel(float* out, const float* in, int rows, int cols) {
int row = blockIdx.y * blockDim.y + threadIdx.y; // K.3's own convention
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < rows && col < cols) {
int k = row;
if (col < k) k = col;
if (rows - 1 - row < k) k = rows - 1 - row;
if (cols - 1 - col < k) k = cols - 1 - col;
int top = k, bottom = rows - 1 - k, left = k, right = cols - 1 - k;
long long ring_start = (long long)rows * cols - (long long)(rows - 2 * k) * (cols - 2 * k);
int width = right - left + 1;
int height = bottom - top + 1;
long long rank;
if (height == 1) rank = ring_start + (col - left);
else if (width == 1) rank = ring_start + (row - top);
else if (row == top) rank = ring_start + (col - left);
else if (col == right) rank = ring_start + (width - 1) + (row - top);
else if (row == bottom) rank = ring_start + (width - 1) + (height - 1) + (right - col);
else rank = ring_start + 2 * (width - 1) + (height - 1) + (bottom - row);
out[rank] = in[row * cols + col];
}
}
// Host-side simulation of the exact grid/block/thread loop nest above --
// Section K.2 and K.3's own verification technique, applied here.
void simulate_spiral_scatter(const std::vector<float>& in, int rows, int cols,
int block_dim_x, int block_dim_y, std::vector<float>& out) {
out.assign((size_t)rows * cols, 0.0f);
int blocks_x = (cols + block_dim_x - 1) / block_dim_x;
int blocks_y = (rows + block_dim_y - 1) / block_dim_y;
for (int by = 0; by < blocks_y; by++)
for (int bx = 0; bx < blocks_x; bx++)
for (int ty = 0; ty < block_dim_y; ty++)
for (int tx = 0; tx < block_dim_x; tx++) {
int row = by * block_dim_y + ty;
int col = bx * block_dim_x + tx;
if (row < rows && col < cols) {
long long rank = spiral_rank_closed_form(row, col, rows, cols);
out[rank] = in[(size_t)row * cols + col];
}
}
}
int main() {
printf("=== Section K.6: spiral order's closed form, cross-checked exhaustively ===\n\n");
printf("Canonical 3x3 check (values 1..9), sequential order:\n");
auto canon = spiral_order_sequential(3, 3);
printf(" ");
for (auto& [r, c] : canon) printf("%d ", r * 3 + c + 1);
printf("\n (expected: 1 2 3 6 9 8 7 4 5)\n\n");
bool all_shapes_ok = true;
int shapes[][2] = {{1,1},{1,5},{5,1},{3,3},{4,4},{4,5},{5,4},{7,7},{6,9},{9,6}};
for (auto& s : shapes) all_shapes_ok &= cross_check_shape(s[0], s[1]);
printf("\nall %d shapes: closed form and sequential simulation agree everywhere: %s\n\n",
(int)(sizeof(shapes) / sizeof(shapes[0])), all_shapes_ok ? "true" : "false");
printf("=== The CUDA kernel: genuinely simulated across divisible and non-divisible shapes ===\n");
bool all_kernel_ok = true;
struct { int rows, cols, bx, by; } kernel_cases[] = {
{4, 4, 4, 4}, {7, 5, 4, 4}, {9, 6, 4, 4}, {6, 9, 4, 4},
};
for (auto& kc : kernel_cases) {
std::vector<float> in((size_t)kc.rows * kc.cols);
for (size_t i = 0; i < in.size(); i++) in[i] = (float)i;
std::vector<float> kernel_out;
simulate_spiral_scatter(in, kc.rows, kc.cols, kc.bx, kc.by, kernel_out);
auto sequential = spiral_order_sequential(kc.rows, kc.cols);
std::vector<float> reference((size_t)kc.rows * kc.cols);
for (size_t rank = 0; rank < sequential.size(); rank++) {
auto [r, c] = sequential[rank];
reference[rank] = in[(size_t)r * kc.cols + c];
}
bool match = (kernel_out == reference);
all_kernel_ok &= match;
printf(" %dx%d, block=%dx%d: kernel simulation matches sequential-derived reference: %s\n",
kc.rows, kc.cols, kc.bx, kc.by, match ? "true" : "false");
}
printf("\nall kernel cases verified correct: %s\n", all_kernel_ok ? "true" : "false");
bool all_ok = all_shapes_ok && all_kernel_ok;
return all_ok ? 0 : 1;
}
07_diagonal_traversal.cu¶
#include <cstdio>
#include <vector>
#include <utility>
// Appendix K.7 -- Section K.6 found that spiral order, despite looking
// like it has no closed form, still produces a STRICT TOTAL ORDER: every
// position has exactly one predecessor and one successor, and nothing
// about spiral order lets two positions be computed independently of one
// another. Diagonal order is the opposite kind of traversal entirely, and
// it is the one CUDA's actual wavefront/DP kernels are built around: every
// position on the SAME anti-diagonal (row + col == d) is independent of
// every other position on that diagonal, and depends only on positions on
// the diagonal immediately before it. That difference -- not a difference
// of formula complexity -- is what decides whether a traversal can be one
// kernel launch (Chapter 18's ordinary elementwise kernels) or needs one
// launch per diagonal (this section's own kernel).
// The classic zigzag print order (alternating direction each diagonal) --
// genuinely sequential, kept here only as a correctness anchor for the
// closed-form diagonal/position formulas below, the same role Section
// K.6's sequential simulation played for its own closed form.
std::vector<std::pair<int,int>> diagonal_order_zigzag(int rows, int cols) {
std::vector<std::pair<int,int>> order;
if (rows <= 0 || cols <= 0) return order;
for (int d = 0; d <= rows + cols - 2; d++) {
int row_min = (d - cols + 1 > 0) ? (d - cols + 1) : 0;
int row_max = (d < rows - 1) ? d : (rows - 1);
if (d % 2 == 0) {
for (int r = row_max; r >= row_min; r--) order.push_back({r, d - r});
} else {
for (int r = row_min; r <= row_max; r++) order.push_back({r, d - r});
}
}
return order;
}
// Closed form: which diagonal a position belongs to, and its position
// WITHIN that diagonal counting from the diagonal's own row_min -- unlike
// spiral order, this needs no ring peeling at all, because a diagonal's
// membership depends only on row + col, not on distance from an edge.
int diagonal_of(int row, int col) { return row + col; }
int row_min_of_diagonal(int d, int rows, int cols) {
(void)rows;
return (d - cols + 1 > 0) ? (d - cols + 1) : 0;
}
int row_max_of_diagonal(int d, int rows, int cols) {
(void)cols;
return (d < rows - 1) ? d : (rows - 1);
}
// ---- The CUDA side: a genuine wavefront dependency, one launch per diagonal ----
//
// Lattice path counting -- paths(r, c) is the number of right/down-only
// paths from (0,0) to (r,c): paths(r,0) = paths(0,c) = 1, and otherwise
// paths(r,c) = paths(r-1,c) + paths(r,c-1). Both of a cell's dependencies
// sit on diagonal (row+col-1) -- one diagonal earlier -- and EVERY cell on
// diagonal (row+col) depends only on that earlier diagonal, never on
// another cell on its own diagonal. That is exactly why this kernel is
// launched once per diagonal instead of once for the whole grid: nothing
// inside one launch may safely assume another thread in the SAME launch
// has already finished, but Chapter 18's ordinary elementwise kernels
// never need that assumption to begin with.
__global__ void dp_wavefront_kernel(float* grid, int rows, int cols, int d) {
int row_min = (d - cols + 1 > 0) ? (d - cols + 1) : 0;
int row_max = (d < rows - 1) ? d : (rows - 1);
int t = blockIdx.x * blockDim.x + threadIdx.x;
int row = row_min + t;
if (row > row_max) return; // boundary check: past this diagonal's last valid row
int col = d - row;
if (row == 0 || col == 0) {
grid[row * cols + col] = 1.0f;
} else {
grid[row * cols + col] = grid[(row - 1) * cols + col] + grid[row * cols + (col - 1)];
}
}
// Host-side simulation of launching dp_wavefront_kernel once per diagonal,
// d = 0 .. rows+cols-2, each launch itself simulated as a full grid/block/
// thread loop nest -- Section K.2/K.3/K.6's own verification technique,
// applied here across a SEQUENCE of launches rather than just one.
void simulate_wavefront_dp(int rows, int cols, int block_dim_x, std::vector<float>& grid) {
grid.assign((size_t)rows * cols, 0.0f);
for (int d = 0; d <= rows + cols - 2; d++) {
int row_min = (d - cols + 1 > 0) ? (d - cols + 1) : 0;
int row_max = (d < rows - 1) ? d : (rows - 1);
int diagonal_len = row_max - row_min + 1;
int blocks_x = (diagonal_len + block_dim_x - 1) / block_dim_x;
for (int bx = 0; bx < blocks_x; bx++) {
for (int tx = 0; tx < block_dim_x; tx++) {
int t = bx * block_dim_x + tx;
int row = row_min + t;
if (row > row_max) continue;
int col = d - row;
if (row == 0 || col == 0) grid[(size_t)row * cols + col] = 1.0f;
else grid[(size_t)row * cols + col] = grid[(size_t)(row - 1) * cols + col] +
grid[(size_t)row * cols + (col - 1)];
}
}
}
}
void scalar_dp_reference(int rows, int cols, std::vector<float>& grid) {
grid.assign((size_t)rows * cols, 0.0f);
for (int r = 0; r < rows; r++)
for (int c = 0; c < cols; c++)
grid[(size_t)r * cols + c] = (r == 0 || c == 0) ? 1.0f
: grid[(size_t)(r - 1) * cols + c] + grid[(size_t)r * cols + (c - 1)];
}
int main() {
printf("=== Section K.7: diagonal order, and the wavefront dependency it exposes ===\n\n");
printf("Canonical 3x3 check (values 1..9), zigzag diagonal order:\n");
auto canon = diagonal_order_zigzag(3, 3);
printf(" ");
for (auto& [r, c] : canon) printf("%d ", r * 3 + c + 1);
printf("\n (expected: 1 2 4 7 5 3 6 8 9)\n\n");
printf("=== Closed-form diagonal membership, cross-checked against the zigzag simulation ===\n");
bool membership_ok = true;
int shapes[][2] = {{3,3},{4,5},{5,4},{1,6},{6,1},{7,7}};
for (auto& s : shapes) {
int rows = s[0], cols = s[1];
auto order = diagonal_order_zigzag(rows, cols);
bool shape_ok = true;
int idx = 0;
for (int d = 0; d <= rows + cols - 2; d++) {
int row_min = row_min_of_diagonal(d, rows, cols);
int row_max = row_max_of_diagonal(d, rows, cols);
for (int k = 0; k < row_max - row_min + 1; k++) {
auto [r, c] = order[idx++];
if (diagonal_of(r, c) != d) shape_ok = false;
}
}
shape_ok &= (idx == rows * cols);
membership_ok &= shape_ok;
printf(" %dx%d: every position's zigzag-visit position falls on its own row+col diagonal: %s\n",
rows, cols, shape_ok ? "true" : "false");
}
printf("\nall shapes: closed-form diagonal membership matches the zigzag order exactly: %s\n\n",
membership_ok ? "true" : "false");
printf("=== The wavefront DP kernel: genuinely simulated launch-by-launch, per diagonal ===\n");
bool all_dp_ok = true;
struct { int rows, cols, bx; } dp_cases[] = { {5, 5, 4}, {6, 9, 4}, {9, 6, 4}, {1, 8, 4}, {8, 1, 4} };
for (auto& dc : dp_cases) {
std::vector<float> wavefront_grid, reference_grid;
simulate_wavefront_dp(dc.rows, dc.cols, dc.bx, wavefront_grid);
scalar_dp_reference(dc.rows, dc.cols, reference_grid);
bool match = (wavefront_grid == reference_grid);
all_dp_ok &= match;
printf(" %dx%d, block_dim_x=%d, %d launches (one per diagonal): matches scalar DP reference: %s\n",
dc.rows, dc.cols, dc.bx, dc.rows + dc.cols - 1, match ? "true" : "false");
}
printf("\nall DP cases verified correct: %s\n", all_dp_ok ? "true" : "false");
printf("\n**[COMMON TRAP]** dp_wavefront_kernel cannot be launched once for the whole grid the way\n");
printf("Chapter 18's elementwise kernels are -- a single launch gives CUDA zero obligation to have\n");
printf("any particular thread, let alone an entire earlier diagonal's worth of threads, finish\n");
printf("before another thread reads grid[row-1][col] or grid[row][col-1]. Each case above needed\n");
printf("exactly rows+cols-1 separate launches, one per diagonal, for exactly one reason: a kernel\n");
printf("launch boundary is the only guarantee, in plain CUDA C++ without cooperative-groups grid\n");
printf("synchronization, that every thread from one launch has completed before the next launch's\n");
printf("threads begin.\n");
bool all_ok = membership_ok && all_dp_ok;
return all_ok ? 0 : 1;
}
Chapter Summary¶
Row-major storage (Chapter 3, Chapter 6) fixes where a matrix's values live; it says nothing about the order a loop nest is allowed to visit them in, and Section K.1's deterministic cache simulator measured the consequence exactly: a traversal matching storage order reached a 93.75% hit rate on an 8x64 matrix with a 4-line cache, while a traversal crossing that same storage reached exactly 0%, on the identical data, computing the identical sum. A genuinely square shape mirrors row-major and column-major results exactly; an asymmetric shape does not, because the severity of crossing a layout depends on the resulting stride relative to one cache line, not merely on whether the order was crossed. Section K.2 closed a gap this book's own Appendix E.3 named but never built: a shared-memory tiled transpose that is both bank-conflict-free (Appendix C.3's exact [32][32] to [32][33] padding fix, confirmed unchanged) and fully boundary-checked for any non-square, non-tile-divisible matrix shape, verified across four shapes by a host-side simulation of the real kernel's own grid/block/tile loop nest. Section K.3 extended Chapter 18.1's ceiling-division launch configuration to two axes, exhaustively confirmed that the resulting grid covers an M x N matrix exactly once for four differently-shaped cases, and computed real addresses proving why x conventionally indexes columns: consecutive threadIdx.x values produce consecutive, coalescable addresses only when x maps to the contiguous axis. Section K.4 showed a boundary check that guards only one of two axes is not a fixed bug with a fixed wrong answer but a genuine, scheduling-dependent data race — silently self-correcting in one execution order and fully visible in the reverse of that same order, with only the final row's overflow guaranteed to survive both — and separately proved that the identical missing check is completely harmless for element-wise addition specifically, because addition's flat-index aliasing happens to read and write the mathematically correct neighboring value regardless. Section K.5 closed the appendix's first arc where Chapter 7.2 first flagged it: cuBLAS's column-major convention needs every one of Section K.3's axis conventions reversed, a real converter moves data between the two layouts and back without loss, and skipping that conversion reproduces Chapter 13.3's own naive-reshape trap — a silent transpose — in the one place this book will actually run into it. Section K.6 proved that spiral order, despite looking inherently sequential, has a genuine O(1) closed form derived from a position's distance to the nearest edge, cross-checked exhaustively against a sequential simulation across ten shapes and then handed to a real CUDA kernel where every thread computes its own rank with no dependency on any other thread at all. Section K.7 closed the appendix with the opposite lesson: diagonal order's positions are independent WITHIN a diagonal but genuinely dependent ACROSS diagonals, demonstrated with a real lattice-path dynamic-programming kernel — paths(r,c) = paths(r-1,c) + paths(r,c-1) — that cannot be launched once for the whole grid the way this book's ordinary elementwise kernels are, and was instead genuinely simulated as the exact sequence of per-diagonal launches (rows+cols-1 of them) it would actually need, verified against a plain scalar DP fill across five shapes.
Self-Check Questions¶
- A row-major matrix is traversed with
colas the outer loop androwas the inner loop. Does this produce the correct sum over all elements? Does it produce good cache behavior? Justify both answers separately. - Section K.1 found that column-major storage's crossing traversal (50% hit rate on an 8x64 shape) is not as damaging as row-major storage's crossing traversal (0% hit rate) on the same shape. What single quantity, computed from the stride and the cache line size, explains the difference?
- Appendix C.3's
shared_transpose_naive/shared_transpose_paddedkernels take a singlewidthparameter. Section K.2's kernels take separaterowsandcolsparameters instead. Why does a singlewidthgenuinely suffice for Appendix C.3's case but not for Section K.2's? - In Section K.3's launch-configuration convention, which CUDA built-in variable indexes a matrix's columns, and why does making the opposite choice still produce a correct result?
- Section K.4 found that the SAME missing boundary check is harmless for element-wise addition but genuinely corrupting for the position-dependent fill kernel. What property of a kernel's written value determines which of those two outcomes occurs?
- Why does Section K.4's forward-order simulation show zero visibly corrupted interior rows, while the reverse-order simulation shows every interior row corrupted? Is either order something a real CUDA launch guarantees?
- A programmer converts a row-major
Tensorto column-major layout before calling a cuBLAS routine, but forgets the conversion on the way back for the routine's OUTPUT. What will they observe in the result, and why? - Section K.6 claims spiral order has a "closed form" even though computing a position's rank requires first computing its ring number
k. In what sense is this still a closed form, and how does it differ from a traversal that genuinely has none? naive_conv2d_kernel(Chapter 18.3) reads several neighboring INPUT positions per thread, anddp_wavefront_kernel(Section K.7) reads several neighboring GRID positions per thread. Both read more than one memory location. Why does Chapter 18.3's kernel work correctly as a single launch while Section K.7's kernel does not?
Where We Go Next¶
Part 7's own Chapter 22 prices options and other instruments over genuinely large simulated paths — batches of matrices large enough that Section K.1's traversal-order argument and Section K.3's launch-configuration accounting stop being an academic exercise and start being the difference between a pricing run that finishes and one that doesn't. Appendix J's pinned-versus-pageable memory distinction is the next piece that same real-hardware pricing pass will need: getting a large batch of these correctly-traversed, correctly-tiled matrices onto a device asynchronously in the first place.
Worked Solutions¶
1. Yes, the sum is correct regardless of loop order — Section K.1's own correctness check confirmed both orders produce the identical sum over the identical data, since addition is commutative and every element is still visited exactly once either way. No, it does not produce good cache behavior on row-major-stored data: the outer col loop strides across rows at a fixed column, and for a matrix whose row length exceeds the cache's line size, each such step lands on a different, likely-evicted cache line, exactly the 0%-hit-rate result Worked Example K.1.1 measured for row-major storage's crossing traversal.
2. The stride relative to the cache line size. Row-major storage's crossing traversal strides by cols = 64 floats — exactly 4 lines, at or beyond the simulator's own 4-line capacity, so every one of the 8 rows needed for one column occupies a distinct line and none can survive to be reused. Column-major storage's crossing traversal (on the same 8x64 shape) strides by only rows = 8 floats — half of one 16-float line — so two consecutive columns' accesses can still share a line before it gets evicted, producing partial reuse and the measured 50% hit rate instead of 0%.
3. Appendix C.3's kernels transpose a SQUARE matrix (an implicit assumption baked into using one width for both dimensions), so the output's row count and column count are both already width, and no separate tracking of "the other dimension" is ever needed. Section K.2's kernels handle a rows x cols matrix whose output is cols x rows — a genuinely different shape whenever rows != cols — so the output's bounds check needs its own, independently-tracked pair of dimensions; collapsing them back into one width would silently assume squareness Section K.2 was built specifically to not assume.
4. threadIdx.x, combined with blockIdx.x and blockDim.x, indexes a matrix's columns in this book's own established convention (every 2D kernel since Chapter 18.3). Making the opposite choice (x indexing rows instead) still produces a CORRECT result, because Worked Example K.3.1's exhaustive coverage check depends only on every (blockIdx, threadIdx) combination being enumerated and bounds-checked against the right matrix dimension somewhere — it does not depend on which named CUDA variable supplies which coordinate. What breaks under the opposite choice is coalescing, not correctness: Section K.3's own genuinely computed addresses show consecutive threadIdx.x values producing consecutive (stride-1) addresses only under the conventional assignment, and addresses cols apart under the swapped one.
5. Whether the value a thread writes is computed by RE-READING through the same flat-index formula the write itself uses, or is instead computed directly from that thread's own (possibly out-of-bounds) (row, col) identity. Element-wise addition's out[row*cols+col] = a[row*cols+col] + b[row*cols+col] re-reads a and b through the identical formula, so an overflowing index aliases to the same flat slot in a, b, AND out simultaneously and recomputes that slot's own correct value. The fill kernel's out[row*cols+col] = row*1000+col instead bakes the calling thread's own (wrong) row and col directly into the value, so an aliased slot receives a value that was never the correct one for that slot to begin with.
6. In forward order, blockIdx.y = 0 is processed first, so row 0's overflow lands on row 1 BEFORE row 1's own block runs and writes its own correct values there — the correct write happens later and silently overwrites the earlier bad one, for every interior row, in sequence. Reverse order processes the highest blockIdx.y first, so each row's own correct write happens BEFORE the row above it (processed later in this order) overflows onto it — nothing ever comes along afterward to overwrite the damage. Neither order is something a real CUDA launch guarantees: CUDA's programming model makes no promise about the order or concurrency of block execution, so both orders (and interleavings of them) are equally legal on real hardware.
7. They will observe the TRANSPOSE of the intended result matrix, not the intended matrix itself, and no error or crash of any kind — Section K.5's own [COMMON TRAP] proved this exactly: reading a column-major buffer's raw bytes as if they were row-major (or the reverse) is a relabeling of which index is contiguous, and swapping which index is contiguous is the literal definition of transpose. The bytes themselves are entirely correct; only the layout tag attached to them is wrong, which is precisely why this class of bug produces a plausible-looking, silently wrong matrix rather than any diagnosable failure.
8. It is a closed form in the sense that matters for parallelizing it: computing spiral_rank_closed_form(row, col, rows, cols) needs only row, col, rows, and cols — no other position's rank, and no simulation of any position visited earlier. That every thread can compute k and its ring offset independently, using nothing but its own coordinates, is exactly what makes the CUDA kernel possible without any inter-thread communication. A traversal with genuinely no closed form would instead require each rank to be computed FROM the previous rank (or from a data structure built by walking earlier positions first), the way spiral_order_sequential's own boundary-shrinking loop works — extra arithmetic inside a closed form is not the same obstacle as a genuine sequential dependency between positions.
9. naive_conv2d_kernel's neighboring reads are all against input — a buffer that is already fully populated before the kernel launches and that no thread in the kernel ever writes to, so there is no possibility that one thread reads a value another thread in the same launch hasn't produced yet. dp_wavefront_kernel's neighboring reads are against grid — the exact buffer OTHER THREADS in a genuinely single launch would be concurrently writing, specifically the threads assigned to the diagonal immediately before the reading thread's own diagonal. CUDA gives no guarantee within one launch that any other thread, let alone an entire earlier diagonal's worth of them, has finished before the current thread runs — which is precisely why Section K.7 needed a separate launch per diagonal, while Chapter 18.3 never needed more than one launch for its entire convolution.