Lab 1.2 — The Memory Hierarchy: Coalescing & Tiled Matrix Multiply
← Course 4 syllabus · Module 1 · Prev: « Lab 1.1 · Next: Lab 1.3 »
Goal
The GPU memory hierarchy, learned the same way Course 3 Module 2 learned the CPU’s cache ladder: by probing it with access patterns and reading the cliffs. Three ideas carry the lab: global-memory coalescing (a warp’s 32 loads becoming few wide transactions — or many narrow ones), shared memory as a programmer-managed cache with bank conflicts as its failure mode, and the classic naive → tiled matrix multiply ladder that turns those ideas into an order-of-magnitude speedup. The philosophical contrast with Course 3 is the point to hold onto: the CPU hid your bad access patterns with caches and speculation; the GPU hides latency with parallelism but makes bandwidth waste entirely your problem. Same experiment, different machine philosophy. Both tracks build it — CUDA C++ for the full ladder, Python for the Numba shared-memory version and the CuPy/cuBLAS ceiling you’re climbing toward.
Recommended reading
- Motta — the material on the CUDA memory hierarchy, global memory access, and shared-memory tiling (title-level reference; confirm against the copy in hand).
- CUDA C++ Programming Guide — the Memory Hierarchy section of the programming model, and the device-memory-access part of the Performance Guidelines chapter (coalescing rules, shared-memory banks).
- CUDA C++ Best Practices Guide — Coalesced Access to Global Memory and Shared Memory (including the shared-memory-in-matrix-multiply worked discussion — read it for the shape of the argument, then put it away and build yours).
- Numba CUDA documentation — Shared memory and thread synchronization (
cuda.shared.array,cuda.syncthreads()). - CuPy documentation — the linear-algebra routines (
cupy.matmul) and what they dispatch to (cuBLAS); worth knowing what you’re being compared against. - Nsight Compute documentation — the Memory Workload Analysis section reference, skimmed before the profiling step so its tables are recognizable.
Prerequisites
Project & environment setup
- New C++ targets in
cuda/:stride_bench,matmul(one executable, kernel selectable by flag is fine). SameCMAKE_CUDA_ARCHITECTURES 89block as Lab 1.1. - Compile with
--ptxas-options=-von for this lab — the per-kernel register and shared-memory report becomes input to Lab 1.5. - Python side: same venv; a
matmul_numba.pyand a short CuPy comparison script. Plots (bandwidth vs. stride) via matplotlib intobenchmarks/.
Where results go:
| Artifact | Path |
|---|---|
| Notes, predicted-vs-measured, bank-conflict postmortem | labs/lab-1-2/notes.md |
Nsight Compute reports (naive vs. tiled), ptxas -v output |
labs/lab-1-2/captures/ |
| Stride-sweep CSV + plot, matmul timing tables (C++, Numba, CuPy) | labs/lab-1-2/benchmarks/ |
Background
- Coalescing. A warp issues 32 addresses per load instruction; the hardware services them as some number of wide memory transactions. Consecutive 4-byte addresses → few transactions, full bus utilization; a stride walks the same data across many transactions, and effective bandwidth divides accordingly. This is Course 3’s stride experiment with the mechanism moved from the cache line to the warp transaction — predict the cliff positions before measuring.
- Shared memory. A small on-SM scratchpad, allocated per block, explicitly loaded and explicitly synchronized (
__syncthreads()). It is the GPU’s answer to “I will touch this data more than once”: stage a tile once from global memory, then reuse it at register-adjacent latency. - Bank conflicts. Shared memory is divided into 32 banks (4-byte words, striped). When threads of a warp hit distinct words in the same bank, accesses serialize. The classic symptom appears in column accesses of a square tile; the classic cure is padding the tile’s row length by one. This lab measures the conflict, not just names it.
- The matmul ladder. For \(C = AB\) with \(N \times N\) floats, the FLOP count is fixed at \(2N^3\). What the ladder changes is bytes: the naive kernel re-reads \(A\) and \(B\) from global memory for every product term — about \(2N^3 \cdot 4\) bytes, an arithmetic intensity near \(\tfrac{2N^3}{8N^3} = 0.25\) FLOP/byte. A shared-memory tile of width \(T\) makes each loaded element serve \(T\) multiply-adds, cutting global traffic to roughly \(2N^3 \cdot 4 / T\) bytes and raising intensity to about \(T/4\) FLOP/byte. Derive this in your notes before writing either kernel — it is the predicted-speedup row, and Lab 1.5 places both kernels on the roofline with it.
Tasks
CUDA C++
- Stride microbenchmark. A copy kernel reading with stride \(s \in \{1, 2, 4, \dots, 32\}\) (fixed total elements touched), event-timed, effective bandwidth computed per stride. Plot bandwidth vs. stride; annotate the plot with where the transaction model predicts each drop.
- Naive matmul. One thread per \(C_{ij}\), inner loop over \(k\), straight from global memory. Verify against a CPU reference at \(N = 256\) (tolerance argued from float rounding, Course 1 §3 — \(N\) same-magnitude accumulations), then event-time at \(N = 512\) and \(1024\).
- Tiled matmul. The shared-memory version: each block stages a \(T \times T\) tile of \(A\) and of \(B\), synchronizes, accumulates, advances. Structure is specified here; the kernel is yours. Choose \(T\) against your device table’s shared-memory-per-block figure and show the arithmetic. Same verification, same timings.
- Bank-conflict experiment. Take the tiled kernel’s transposed-access variant (or your \(B\)-tile column reads), measure with and without the
+1padding column, and capture the conflict counter difference in Nsight Compute — the point is seeing the counter move, not folklore about padding.
CUDA Python
- Numba tiled matmul. The same tiled structure with
cuda.shared.arrayandcuda.syncthreads(). Verify against NumPy; event-time it. Note where Numba made you hold tile dimensions as compile-time constants and what that implies. - The library ceiling. Time
cupy.matmulon the same sizes (events, warm-up excluded). This is cuBLAS — record it as the ceiling row of your table, and read off how far your best hand-written rung still is from it.
Deliverable & expected results
- Stride plot with annotated cliffs; matmul ladder table (naive C++, tiled C++, tiled Numba,
cupy.matmul) at \(N = 512, 1024\); verified correctness for every rung; bank-conflict before/after numbers.
| Quantity | Predicted | Measured |
|---|---|---|
| Bandwidth vs. stride shape | near-full bandwidth at stride 1, stepwise collapse as strides fragment warp transactions, floor by stride 32 — cliff positions from the transaction-size argument | … |
| Naive matmul arithmetic intensity | \(\approx 0.25\) FLOP/byte (\(2N^3\) FLOPs over \(\approx 2N^3\cdot4\) B) → firmly memory-bound | … |
| Tiled vs. naive speedup | large, direction certain — global traffic cut by \(\approx T\times\); the measured factor is what reconciliation explains | … |
| Bank-conflict padding effect | measurable improvement on the conflicted variant, conflict counters → ~0 with padding | … |
cupy.matmul vs. best hand kernel |
cuBLAS wins — it layers register tiling and deeper pipelining on the same ideas | … |
Profiling & performance
Nsight Compute, first serious use: ncu --section MemoryWorkloadAnalysis --section SpeedOfLight -o <report> ./matmul ... for the naive and tiled kernels. Read side by side: global load transactions per request (the coalescing verdict), achieved memory throughput vs. peak (naive should pin it; tiled should drop it while getting faster — say why that’s the healthy direction), and the shared-memory bank-conflict counters for task 4. Save both reports to captures/.
Analysis & reconciliation
Reconcile the stride plot against the transaction model and against Course 3’s CPU version of the same experiment — one paragraph on where the cliffs sit differently and which machine forgave which pattern. For matmul: from measured naive bandwidth and your \(T/4\) intensity derivation, compute the speedup the tile should have bought, compare with what it did, and name the leftovers (sync overhead, occupancy limits from ptxas -v — Lab 1.5’s territory). Close by placing your tiled kernel as a fraction of cupy.matmul and listing, from the Nsight report, the top candidate reason for the remaining gap.
Going further
- Rectangular tiles and a register-blocked inner loop (each thread computes a small sub-tile of \(C\)) — the next rung cuBLAS climbed.
- Re-run the stride sweep on pitched allocations (
cudaMallocPitch) and see what the alignment guarantee is worth for 2-D data. - Try the same stride sweep with the L2 hit counters visible in Nsight Compute — Ada’s L2 is generous (read the exact size off your device query), and it’s instructive to see where it quietly helped.