1.8 GPU Architecture & The Memory Hierarchy¶

Every claim you will ever read about LLM training and inference — “this kernel is memory-bound,” “decode is bandwidth-limited,” “we need to keep the tensor cores fed,” “the KV cache spilled to a slower tier” — is downstream of the hardware described in this chapter. The GPU is not a magic box that does matrix multiplication fast. It is a very specific machine with a very specific shape: thousands of arithmetic units starved by a comparatively thin pipe to memory, organized into a strict hierarchy of caches and scratchpads, executing instructions in lockstep groups of 32. If you understand that shape, the entire performance-engineering half of this book becomes mechanical. If you do not, you will spend your career confused about why your “compute-heavy” workload runs at 5% of peak FLOPs.

This is the single most important systems chapter in the book. Almost everything in Part IV — Kernels, Efficiency & Quantization and Part VII — Inference & Serving is an exercise in respecting the constraints we lay out here. We will build up the GPU from threads to warps to streaming multiprocessors (SMs); walk down the memory hierarchy from registers to high-bandwidth memory (HBM); develop the central concept of arithmetic intensity that decides whether you are compute-bound or memory-bound; and ground all of it in the concrete specifications of the NVIDIA A100, H100, H200, and B200 datacenter GPUs. We assume you know Python and the linear algebra of Linear Algebra for Deep Learning; we explain everything GPU-specific.

Cross-references: Numerical Computing, Floating Point & Precision explains the FP16/BF16/FP8 number formats whose throughput we quote here; Parallel Computing & Collective Communication takes the NVLink and multi-GPU story to the cluster scale; The Roofline Model & Performance Engineering formalizes arithmetic intensity into the roofline; FlashAttention I: IO-Awareness & The Online Softmax is the canonical case study of trading FLOPs for memory traffic; and CUDA Programming Essentials for ML Engineers turns the execution model below into real kernels.

Why a GPU Looks the Way It Does¶

A CPU is optimized for latency: finish one thread’s instruction stream as fast as possible. It spends most of its transistor budget on machinery to make a single thread fast — large caches, branch predictors, out-of-order execution, deep speculation. A modern CPU core might have a few hundred kilobytes of private cache and devote enormous silicon to avoiding memory stalls for one thread.

A GPU is optimized for throughput: finish an enormous number of independent threads per second, and tolerate each individual thread being slow. It spends its transistor budget on arithmetic units and on the ability to keep thousands of threads in flight simultaneously. When one thread stalls waiting on memory, the hardware instantly switches to another ready thread. This is latency hiding through massive parallelism, and it is the single design decision from which everything else follows.

The consequence: a GPU only goes fast when you give it tens of thousands of independent units of work. A matrix multiply of two large matrices is perfect — billions of independent multiply-accumulates. Autoregressive decode of one token at a time for a single request is nearly the worst case — almost no parallelism, so the machine sits idle hiding latency it cannot hide. Half of inference engineering (Continuous Batching & Request Scheduling) is manufacturing parallelism the workload does not naturally have.

        CPU                                  GPU
  +---------------+              +-------------------------------+
  | big control   |              | tiny control, replicated many |
  | huge caches   |              | small caches, many ALUs        |
  | few fat cores |              | thousands of thin "cores"      |
  +---------------+              +-------------------------------+
  goal: 1 thread fast           goal: 10,000s of threads in flight
  hides latency by              hides latency by SWITCHING to
  predicting / caching          another ready warp every cycle

A GPU is also a co-processor. It hangs off the CPU (the “host”) over PCIe or, on Grace-Hopper-class systems, a fast chip-to-chip link. The host launches kernels (GPU functions) and shuffles data to and from GPU memory. A recurring performance bug is a workload that is secretly host-bound or PCIe-bound: the GPU is fast but it spends its life waiting for the CPU to feed it.

The Execution Model: Threads, Warps, Blocks, and SMs¶

NVIDIA’s programming model (CUDA — Compute Unified Device Architecture) exposes a hierarchy of parallel work. Understanding it is non-negotiable.

Threads and the SIMT model¶

The smallest unit is a thread. A thread runs the kernel body once, with its own registers and its own index, and typically computes one element (or a few) of the output. You launch many thousands of them.

Threads are grouped into warps of exactly 32 threads. This number is a hardware constant on every NVIDIA GPU to date and it dominates everything. The 32 threads in a warp execute in lockstep: they share one instruction fetch/decode unit and one program counter, and on each cycle all 32 lanes execute the same instruction on different data. NVIDIA calls this SIMT — Single Instruction, Multiple Threads. It is SIMD (Single Instruction, Multiple Data) with a thread-shaped programming abstraction layered on top.

Two consequences of the warp being the real unit of execution:

• Warp divergence. If threads in a warp take different sides of an if, the warp executes both paths serially, masking off the inactive lanes on each path. A branch that splits a warp 50/50 can halve throughput. Branchy, data-dependent code is GPU-hostile.
• Coalesced memory access. When a warp issues a load, the hardware tries to service all 32 lanes with as few memory transactions as possible. If lane \(i\) reads address base + i*4 (32 consecutive 4-byte words), the whole warp is satisfied by one or two 128-byte transactions — coalesced, near-peak bandwidth. If the 32 lanes read scattered addresses, you get up to 32 separate transactions and a 32× bandwidth penalty.

Blocks (CTAs) and the SM¶

Warps are grouped into a thread block, also called a cooperative thread array (CTA). A block is the unit of scheduling and the unit of cooperation: threads within a block can synchronize with a barrier (__syncthreads()) and share data through a fast on-chip scratchpad called shared memory (SMEM). Threads in different blocks cannot cheaply synchronize or share scratchpad; this isolation is what lets the GPU run blocks independently and scale.

Each block is assigned to exactly one streaming multiprocessor (SM) for its entire life. The SM is the fundamental compute building block of the GPU — think of it as a core, but a wide, throughput-oriented one. An H100 has 132 SMs; an A100 has 108. Inside one SM you find:

• A register file (e.g. 65,536 32-bit registers — 256 KB — per SM), partitioned among all resident threads.
• A pool of CUDA cores: scalar ALUs for FP32/INT32 arithmetic.
• Tensor Cores: dedicated matrix-multiply-accumulate units (more below).
• A block of on-chip SRAM split between shared memory and L1 cache (configurable, e.g. up to 228 KB total per SM on H100).
• Warp schedulers (typically 4 per SM) that pick a ready warp each cycle and issue its next instruction.

              GPU (e.g. H100: 132 SMs)
   +-----------------------------------------------+
   |  SM 0   SM 1   SM 2   ...            SM 131    |
   |   ||      ||                                   |
   |   ||  shared by all SMs:  L2 cache (~50 MB)    |
   |   ||              |                            |
   +-----------------------------------------------+
                       |
                  HBM (off-chip DRAM, 80 GB)

   Inside ONE SM:
   +-----------------------------------------------+
   |  4 warp schedulers                            |
   |  Register file (256 KB)                       |
   |  CUDA cores (FP32/INT)   Tensor Cores         |
   |  Shared memory / L1 (configurable SRAM)       |
   +-----------------------------------------------+

   Software hierarchy:    Grid  >  Block(CTA)  >  Warp(32)  >  Thread
   Hardware mapping:      GPU   >  SM          >  scheduler >  lane

The grid and the launch¶

A kernel launch creates a grid of blocks. You choose the grid dimensions and block dimensions. The GPU’s global scheduler hands blocks out to SMs as they free up. Because there are usually far more blocks than SMs, this gives automatic load balancing and forward scalability: the same kernel runs on a small GPU (few SMs) or a large one (many SMs) with no code change — bigger GPUs simply chew through more blocks concurrently.

# A from-scratch mental model of the GPU execution hierarchy in pure Python.
# This is NOT how you run on a GPU (you'd use CUDA/Triton) -- it is a faithful
# simulation of the indexing and lockstep semantics so the model is concrete.

WARP_SIZE = 32

def launch_kernel(grid_dim, block_dim, kernel, *args):
    """Emulate a CUDA grid launch: grid_dim blocks, each of block_dim threads."""
    for block_id in range(grid_dim):
        # Each block is independent -> could run on any SM, in any order.
        block_threads = list(range(block_dim))
        # Threads execute in warps of 32, in lockstep within a warp.
        for warp_start in range(0, block_dim, WARP_SIZE):
            warp = block_threads[warp_start: warp_start + WARP_SIZE]
            # All lanes in this warp run the SAME instruction stream together.
            for thread_id in warp:               # in HW: simultaneous, not serial
                global_tid = block_id * block_dim + thread_id
                kernel(global_tid, block_id, thread_id, *args)

def saxpy(global_tid, block_id, thread_id, a, x, y, out, n):
    """out[i] = a * x[i] + y[i], one element per thread (the classic 'SAXPY')."""
    i = global_tid
    if i < n:                                    # guard: grid may overshoot n
        out[i] = a * x[i] + y[i]

n = 100
x = [float(i) for i in range(n)]
y = [1.0] * n
out = [0.0] * n
block_dim = 32                                   # one warp per block here
grid_dim = (n + block_dim - 1) // block_dim      # ceil-div = 4 blocks for n=100
launch_kernel(grid_dim, block_dim, saxpy, 2.0, x, y, out, n)
print(out[:5])     # [1.0, 3.0, 5.0, 7.0, 9.0]  == 2*x + 1

This saxpy is the canonical memory-bound kernel: for every element it does 2 FLOPs (one multiply, one add) but reads 8 bytes (x[i], y[i]) and writes 4 bytes. We will quantify exactly why that is bad shortly.

CUDA Cores vs Tensor Cores¶

There are two fundamentally different arithmetic engines on a modern datacenter GPU, and conflating them is the most common beginner error in performance reasoning.

CUDA cores: scalar, general-purpose¶

A CUDA core is a scalar ALU. Each one does one fused multiply-add (FMA) per cycle on FP32 (or FP64 on the smaller FP64 units, or INT32). When you write c = a * b + d in a kernel on plain floats, that runs on CUDA cores. The FP32 throughput of a GPU is roughly (number of FP32 lanes) × (clock) × 2 (the FMA counts as 2 FLOPs). On an A100 this is on the order of 19.5 TFLOP/s of FP32. Respectable, but not where the headline numbers come from.

Tensor Cores: matrix-multiply-accumulate engines¶

A Tensor Core is a small systolic-array-like unit that computes an entire small matrix multiply-accumulate per operation: \(D = A \times B + C\) where \(A\), \(B\), \(C\), \(D\) are small tiles (e.g. \(16\times16\)). Instead of issuing \(16^3\) scalar FMAs, one Tensor Core instruction consumes whole tiles and produces a whole tile. This is why the BF16/FP16 Tensor Core throughput of an A100 (on the order of 312 TFLOP/s) is 16× its FP32 CUDA-core throughput. The H100 pushes this further with FP8 Tensor Cores; the Blackwell B200 adds FP4. Every modern training and inference workload lives or dies on Tensor Core utilization. The Tensor Core is a bolt-on matrix-multiply unit dropped into an otherwise general-purpose, programmable SIMT core; the opposite design point — a chip that is almost entirely systolic array, driven by a compiler instead of a warp/SMEM hierarchy — is Google’s TPU (and AWS Trainium), and that systolic-array-versus-SIMT contrast is developed in The Accelerator Landscape: TPUs, Trainium, AMD/ROCm & Gaudi.

The catch: Tensor Cores only accelerate matrix multiplication (and operations you can phrase as one), and only in reduced precision (BF16/FP16/FP8/FP4 inputs, usually with FP32 accumulation). They want tile-shaped, contiguous data with dimensions that are multiples of 8 or 16. A pointwise operation (add, GELU, layernorm) gets zero benefit from Tensor Cores — it runs on CUDA cores and is almost always memory-bound. This split — matmuls on Tensor Cores, everything else on CUDA cores and bandwidth-limited — is the structural reason kernel fusion (Kernel Fusion, torch.compile, CUDA Graphs & Compilers) matters so much.

Tile and wave quantization: why dimensions want to be round¶

Two “quantization” effects — nothing to do with numeric quantization — fall directly out of tile-shaped hardware, and they explain a surprising amount of real-world model design.

Tile quantization. A Tensor Core GEMM kernel carves the output matrix into fixed-size tiles (typically \(128\times128\), \(128\times256\), or similar, chosen by the library’s heuristics). A dimension that is not a multiple of the tile forces a final, partly-wasted tile: an output width of 257 needs \(\lceil 257/128 \rceil = 3\) tile-columns of work to produce 257 useful columns, so ~33% of the issued work is padding. Efficiency therefore drops in a sawtooth as dimensions cross tile boundaries. This is why hidden sizes, FFN widths, head dimensions, and vocabularies are padded to multiples of 64 or 128 in practice — Megatron-LM exposes exactly this as --make-vocab-size-divisible-by 128, and it is why you see vocabularies like 50,304 (a padded 50,257) rather than the tokenizer’s raw count. We apply the same rounding when choosing Stack-100M’s width and vocabulary (The Stack-100M Architecture: SOTA Components, Cited and Assembled).

Wave quantization. The grid is executed in waves: roughly (number of SMs) × (blocks resident per SM) blocks run concurrently, then the next wave starts. A grid of 133 blocks at one block per SM on a 132-SM H100 runs two waves, the second of which occupies a single SM — nearly half the wall-clock at under 1% utilization. Kernel time therefore rises in a staircase, and a problem size just past a wave boundary is the worst case. Practical consequence: prefer batch/tile shapes whose grid is an exact multiple of the SM count, and be suspicious when a small increase in problem size causes a large jump in time.

import torch

# Demonstrate the CUDA-core vs Tensor-core gap empirically (run on a GPU).
# FP32 matmul uses CUDA cores; BF16 matmul uses Tensor Cores.
def benchmark_matmul(dtype, n=8192, iters=50):
    if not torch.cuda.is_available():
        return None
    a = torch.randn(n, n, device="cuda", dtype=dtype)
    b = torch.randn(n, n, device="cuda", dtype=dtype)
    torch.cuda.synchronize()
    start = torch.cuda.Event(enable_timing=True)
    end = torch.cuda.Event(enable_timing=True)
    _ = a @ b                                  # warm up / autotune
    torch.cuda.synchronize()
    start.record()
    for _ in range(iters):
        c = a @ b
    end.record()
    torch.cuda.synchronize()
    ms = start.elapsed_time(end) / iters
    flops = 2 * n**3                           # 2*n^3 FLOPs for n x n matmul
    tflops = flops / (ms * 1e-3) / 1e12
    return ms, tflops

# On an A100 you would typically see something on the order of:
#   FP32 : ~19  TFLOP/s   (CUDA cores)
#   BF16 : ~280 TFLOP/s   (Tensor Cores) -- over 10x faster on the SAME hardware
for dt in (torch.float32, torch.bfloat16):
    res = benchmark_matmul(dt)
    if res:
        ms, tflops = res
        print(f"{str(dt):16s}  {ms:7.2f} ms  {tflops:7.1f} TFLOP/s")

The reason torch.set_float32_matmul_precision("high") and BF16 autocast exist is precisely to route FP32 logical math onto the Tensor Cores (via TF32 or BF16) and capture that order-of-magnitude speedup. See Mixed Precision, bf16 & FP8 Training.

The Memory Hierarchy: Registers, SMEM, L2, and HBM¶

Compute is cheap; moving data is expensive. The GPU memory hierarchy is a series of tiers, each roughly an order of magnitude larger and an order of magnitude slower than the one above it. The art of GPU kernel writing is keeping data in the fast tiers and minimizing trips to the slow ones.

The numbers are illustrative and vary by exact SKU and clock, but the ratios are the point: SMEM is roughly an order of magnitude faster than L2, which is faster than HBM, which dwarfs PCIe. Each step down is a cliff.

Registers¶

Each thread has private registers. These are the only operands the ALUs read directly, and accessing them is effectively free. The register file is a fixed resource per SM (e.g. 65,536 32-bit registers). If each thread uses 64 registers, an SM can host at most \(65536 / 64 = 1024\) threads. Use too many registers per thread and you reduce occupancy (next section). “Register spilling” — when the compiler runs out of registers and pushes variables to slow “local memory” (which actually lives in L1/L2/HBM) — is a silent performance killer.

Shared memory (SMEM): the programmer-managed scratchpad¶

Shared memory is fast on-chip SRAM, private to a thread block, and — crucially — explicitly managed by you, the kernel author. It is not a cache that fills automatically; you stage data into it on purpose. The canonical use is tiling a matrix multiply: cooperatively load a tile of \(A\) and a tile of \(B\) from HBM into SMEM once, then have all threads in the block reuse those tiles many times from the fast scratchpad before fetching the next tile. This converts a flood of redundant HBM reads into a trickle. SMEM is physically banked (32 banks); if two lanes of a warp hit the same bank with different addresses you get a bank conflict and serialized access — the SMEM analogue of uncoalesced HBM access.

This staged reuse is the entire mechanism behind FlashAttention: keep the attention computation’s working set in SMEM and registers, never materialize the full \(S \times S\) attention matrix in HBM. See FlashAttention I: IO-Awareness & The Online Softmax.

L2 cache¶

The L2 cache is shared by all SMs and caches HBM transparently. You do not manage it directly, but you exploit it by structuring access patterns so that data fetched by one SM is reused by another while still resident. On Hopper, the L2 cache residency controls and the new thread block clusters / distributed shared memory let blocks on nearby SMs share SMEM, effectively a programmable extension of the SMEM tier.

HBM: the headline “GPU memory”¶

When someone says “an 80 GB H100,” they mean its HBM — High Bandwidth Memory, stacked DRAM sitting on the same package as the GPU die, connected by an extremely wide bus. HBM is where your model weights, activations, optimizer states, and KV cache live. Its capacity bounds how big a model you can hold; its bandwidth bounds how fast you can stream those weights through the compute units. For LLM decode, HBM bandwidth is almost always the binding constraint, because each generated token must read the entire model’s weights from HBM at least once and does very little arithmetic with them. This is the central fact of inference economics (The Anatomy of LLM Inference: Prefill, Decode & The KV Cache).

   fast & tiny  ----------------------------------------->  slow & huge
   +-----------+   +------------+   +--------+   +-----------------+
   | registers | < | SMEM / L1  | < |   L2   | < |   HBM (DRAM)    |
   | 256 KB/SM |   | ~228 KB/SM |   | ~50 MB |   |     80 GB       |
   |  ~1 cyc   |   |  ~30 cyc   |   |~200 cyc|   |   ~500 cyc      |
   +-----------+   +------------+   +--------+   +-----------------+
        ^                ^                            |
        |   you stage    |   you cache-friendly       |  you minimize
        |   explicitly   |   structure access         |  total bytes moved

Occupancy: Keeping the Machine Busy¶

The GPU hides memory latency by having many warps resident on each SM: when one warp stalls on a load, a warp scheduler instantly issues from another ready warp. Occupancy is the ratio of resident warps to the hardware maximum per SM. Higher occupancy generally means more latency-hiding headroom.

Occupancy is bounded by whichever per-SM resource runs out first:

• Registers. Resident threads × registers/thread ≤ register file size.
• Shared memory. Resident blocks × SMEM/block ≤ SMEM size per SM.
• Block/warp slots. A hardware cap on resident blocks and warps per SM.

A worked occupancy calculation makes the trade-offs concrete.

The crucial subtlety, established by Vasily Volkov’s work, is that maximum occupancy is not the goal — performance is. A kernel that uses lots of registers per thread to hold a big accumulator tile (low occupancy) can dramatically outperform a high-occupancy kernel, because it gets more reuse per byte loaded. FlashAttention and high-performance GEMM kernels deliberately run at modest occupancy. Occupancy buys you latency hiding; if you already hide latency another way (instruction-level parallelism, big register tiles), you do not need much of it. Treat occupancy as one lever, not the objective.

def max_blocks_per_sm(threads_per_block, regs_per_thread, smem_per_block_bytes,
                      regfile=65536, smem_per_sm=101376,
                      max_threads=2048, max_blocks=32, warp=32):
    """Replicate NVIDIA's occupancy calculator logic: the min over all limits."""
    regs_per_block = threads_per_block * regs_per_thread
    by_regs   = regfile // regs_per_block if regs_per_block else max_blocks
    by_smem   = smem_per_sm // smem_per_block_bytes if smem_per_block_bytes else max_blocks
    by_threads = max_threads // threads_per_block
    blocks = min(by_regs, by_smem, by_threads, max_blocks)
    warps  = blocks * threads_per_block // warp
    occupancy = warps / (max_threads // warp)
    return blocks, warps, occupancy

for regs, smem in [(48, 24*1024), (32, 16*1024), (32, 8*1024)]:
    b, w, occ = max_blocks_per_sm(256, regs, smem)
    print(f"regs/thread={regs:2d} smem={smem//1024:2d}KB -> "
          f"{b} blocks, {w} warps, occupancy={occ:.0%}")

Arithmetic Intensity: The Master Concept¶

Here is the idea that unifies this entire chapter and most of Part IV. Every kernel does some amount of arithmetic and moves some number of bytes. The ratio is its arithmetic intensity (also called operational intensity):

The hardware has a fixed peak compute \(\pi\) (FLOP/s) and a fixed peak bandwidth \(\beta\) (byte/s). Their ratio defines the machine balance or ridge point:

The achievable performance of a kernel is bounded by:

This is the roofline (The Roofline Model & Performance Engineering). Read it as a simple dichotomy:

• If \(I < I_{\text{ridge}}\), the kernel is memory-bound: it would finish faster if HBM were faster, and adding more FLOPs (e.g. recomputation) is “free.” You are bandwidth-limited.
• If \(I > I_{\text{ridge}}\), the kernel is compute-bound: it is limited by the arithmetic units. You are FLOP-limited.

For an A100, \(\pi \approx 312\,\text{TFLOP/s}\) (BF16) and \(\beta \approx 2.0\,\text{TB/s}\), so \(I_{\text{ridge}} \approx 156\) FLOP/byte. For an H100 with \(\pi \approx 990\,\text{TFLOP/s}\) (BF16, dense) and \(\beta \approx 3.35\,\text{TB/s}\), \(I_{\text{ridge}} \approx 295\) FLOP/byte. The ridge point has been climbing every GPU generation — compute grows faster than bandwidth — which means more and more kernels fall on the memory-bound side over time. This is the secular trend in ML systems: we are increasingly bandwidth-starved, and techniques that trade FLOPs for bytes keep winning.

def arithmetic_intensity(flops, bytes_moved):
    return flops / bytes_moved

def roofline_perf(intensity, peak_flops, peak_bw):
    """Attainable FLOP/s under the roofline model."""
    return min(peak_flops, intensity * peak_bw)

# A100-class numbers (BF16 Tensor Core peak, HBM2e bandwidth).
PEAK_FLOPS = 312e12      # 312 TFLOP/s
PEAK_BW    = 2.0e12      # 2.0 TB/s
ridge = PEAK_FLOPS / PEAK_BW
print(f"A100 ridge point: {ridge:.0f} FLOP/byte")

# Case 1: 8192^3 BF16 matmul
N = 8192
mm_flops = 2 * N**3
mm_bytes = 3 * N**2 * 2          # 2 bytes/elem, read A,B write C
I_mm = arithmetic_intensity(mm_flops, mm_bytes)
print(f"matmul  I={I_mm:8.0f} FLOP/byte -> "
      f"{'COMPUTE' if I_mm > ridge else 'MEMORY'}-bound, "
      f"{roofline_perf(I_mm, PEAK_FLOPS, PEAK_BW)/1e12:.0f} TFLOP/s attainable")

# Case 2: 7B decode, batch=1
params = 7e9
dec_bytes = params * 2           # BF16 weights, read once
dec_flops = 2 * params           # ~2 FLOPs/param for the matrix-vector products
I_dec = arithmetic_intensity(dec_flops, dec_bytes)
print(f"decode  I={I_dec:8.1f} FLOP/byte -> "
      f"{'COMPUTE' if I_dec > ridge else 'MEMORY'}-bound, "
      f"per-token >= {dec_bytes/PEAK_BW*1e3:.1f} ms (bandwidth-bound)")

# Case 3: decode at batch B reuses the SAME weight bytes for B tokens
for B in (1, 8, 32, 128, 256):
    I_b = (2 * params * B) / (params * 2)   # bytes constant, flops scale with B
    bound = 'COMPUTE' if I_b > ridge else 'MEMORY'
    print(f"  batch={B:4d}: I={I_b:6.0f} FLOP/byte -> {bound}-bound")

The third loop above shows the magic of batching in one line: weight bytes read are constant in batch size, but useful FLOPs scale linearly with it, so arithmetic intensity scales linearly with batch size until you cross the ridge and finally light up the Tensor Cores.

The Datacenter GPUs: A100, H100, H200, B200, and NVLink¶

You will reason about these chips constantly. Here are the load-bearing specs and what each generation changed. Treat the figures as representative of the SXM datacenter variants; exact numbers vary by SKU, clock, and whether sparsity is assumed (we quote dense compute unless noted).

What each generation actually changed¶

• A100 (Ampere, 2020). Established the modern recipe: BF16/TF32 Tensor Cores, 2:4 structured sparsity, and large HBM. The reference point for almost all published scaling-law and kernel work. Ridge point ≈ 156 FLOP/byte.
• H100 (Hopper, 2022). Roughly tripled BF16 Tensor Core throughput and added FP8 Tensor Cores, doubling effective matmul throughput again for models that tolerate FP8 (Mixed Precision, bf16 & FP8 Training). Introduced the Transformer Engine (automatic per-tensor FP8 scaling), thread block clusters with distributed shared memory (blocks on neighboring SMs share SMEM), and the Tensor Memory Accelerator (TMA) for asynchronous bulk copies between HBM and SMEM that free the warps from address arithmetic. FlashAttention-3 (FlashAttention 2 & 3) is built specifically to exploit TMA, warp specialization, and FP8 on Hopper. But bandwidth only went from 2.0 to 3.35 TB/s while compute tripled — so the ridge point jumped to ≈295, and more workloads became memory-bound.
• H200 (Hopper refresh, 2024). Same compute die as H100 but with HBM3e: 141 GB at ~4.8 TB/s. The extra capacity and bandwidth are aimed squarely at LLM inference, where you are bandwidth- and KV-cache-capacity-bound. More memory means longer contexts and bigger batches per GPU; more bandwidth directly speeds up decode (which is bandwidth-bound, as we computed). Same FLOPs, but materially faster and roomier serving.
• B200 (Blackwell, 2024). A two-die package presented as one logical GPU, with much higher HBM3e capacity (~192 GB) and bandwidth (~8 TB/s), a second-generation Transformer Engine, and new low-precision formats including FP4 for inference. NVLink jumps to ~1.8 TB/s per GPU. The design philosophy is unchanged but the dials are turned up, and FP4 pushes effective inference throughput much higher for models that survive 4-bit (Quantization II: INT4/INT8/FP8, GGUF, bitsandbytes & QAT).
• Blackwell Ultra and beyond (2025–2026). The B200 generation was refreshed as Blackwell Ultra (B300/GB300) — same Blackwell architecture family and HBM3e generation, with more memory and FP4 throughput per GPU — and by 2026 it is the GPU behind NVIDIA’s flagship GB300 NVL72 racks. NVIDIA’s next architecture, Rubin, is slated to ship in the second half of 2026 with HBM4 memory and a roadmap claim of roughly 2.5× Blackwell’s FP4 throughput; as always, treat pre-launch vendor figures as directional until independently benchmarked. The trend from the last four generations — compute and low-precision throughput outrunning bandwidth — shows no sign of reversing.

NVLink and NVSwitch: the intra-node fabric¶

A single GPU cannot hold a frontier model, and even when it can, you want many GPUs cooperating on one forward/backward pass. The link between GPUs inside a server is NVLink — a dedicated high-bandwidth, low-latency interconnect far faster than PCIe. NVLink generations scale from ~600 GB/s per GPU on A100 to ~900 GB/s on H100/H200 to ~1.8 TB/s on B200 (bidirectional aggregate). An NVSwitch is a crossbar that connects all 8 GPUs in a node so any pair gets full NVLink bandwidth — an all-to-all topology rather than a ring.

Why this matters: tensor parallelism (Distributed Training II) splits each weight matrix across GPUs and must do an all-reduce on activations every layer, which is only viable over NVLink, not PCIe — hence tensor parallelism stays within an NVLink domain (typically 8 GPUs), while pipeline and data parallelism span the slower inter-node network (Parallel Computing & Collective Communication). The Grace-Hopper and Blackwell systems extend this with NVLink-C2C (a coherent CPU–GPU link) and NVLink Switch fabrics that let dozens of GPUs share memory at NVLink-class bandwidth, blurring the line between “one big GPU” and “a node.”

# Where does your 70B model physically fit? A capacity + bandwidth sanity check.
def model_footprint_gb(params_b, bytes_per_param):
    return params_b * 1e9 * bytes_per_param / 1e9   # -> GB

def kv_cache_gb(layers, kv_heads, head_dim, seq, batch, bytes_per_elem=2):
    # K and V, both, per layer, per (kv) head: factor of 2 for K&V.
    elems = 2 * layers * kv_heads * head_dim * seq * batch
    return elems * bytes_per_elem / 1e9

# Llama-70B-ish: 70B params, 80 layers, 8 KV heads (GQA), head_dim 128.
for dtype, bpp in [("bf16", 2), ("fp8", 1), ("int4", 0.5)]:
    w = model_footprint_gb(70, bpp)
    print(f"70B weights in {dtype:5s}: {w:6.1f} GB")

kv = kv_cache_gb(layers=80, kv_heads=8, head_dim=128,
                 seq=8192, batch=32, bytes_per_elem=2)
print(f"KV cache (8k ctx, batch 32, GQA-8, bf16): {kv:.1f} GB")
# 70B bf16 weights ~140 GB -> does NOT fit one 80GB H100; fits one 192GB B200,
# or shards across 2x H100 / 1x via int4. Bandwidth then sets decode speed.

This snippet captures the two questions you ask of every (model, GPU) pair: does it fit (HBM capacity, including the KV cache that grows with batch × context — see PagedAttention & KV-Cache Memory Management), and how fast does it decode (HBM bandwidth). H200’s jump to 141 GB and B200’s to 192 GB are aimed directly at the first question; their bandwidth bumps at the second.

Know your own GPU: querying and profiling the real device¶

Spec tables are for reasoning; the GPU in front of you is what you must actually feed. You will most likely train the book’s ~100M-parameter model on a single mid-range GPU (The Pretraining Run: A Complete Single-GPU Training Loop), so learn to interrogate it directly. PyTorch exposes the device properties CUDA reports, and a two-line copy benchmark gives you an achieved bandwidth number. Expect it to come in meaningfully under the spec-sheet peak — a plain copy is not a tuned streaming kernel, clocks throttle, and a shared GPU is slower still — so treat it as a realistic floor for what your own kernels can hope for, and as the \(\beta\) to use when you want a conservative roofline.

import torch

def gpu_report(dev=0):
    """Everything you need to compute your own ridge point, from the driver."""
    if not torch.cuda.is_available():
        print("no CUDA device visible"); return None
    p = torch.cuda.get_device_properties(dev)
    l2 = getattr(p, "L2_cache_size", None)          # present on recent PyTorch
    print(f"{p.name}  (compute capability {p.major}.{p.minor})")
    print(f"  SMs                     : {p.multi_processor_count}")
    print(f"  HBM capacity            : {p.total_memory / 1e9:.1f} GB")
    print(f"  max threads / SM        : {getattr(p, 'max_threads_per_multi_processor', '?')}")
    print(f"  registers / SM          : {getattr(p, 'regs_per_multiprocessor', '?')}")
    print(f"  SMEM / SM               : {getattr(p, 'shared_memory_per_multiprocessor', '?')} B")
    if l2: print(f"  L2 cache                : {l2 / 1e6:.0f} MB")
    return p

def measured_hbm_bandwidth(n_bytes=1 << 30, iters=20):
    """Empirical HBM bandwidth: a big device-to-device copy reads AND writes."""
    if not torch.cuda.is_available():
        return None
    x = torch.empty(n_bytes // 2, dtype=torch.float16, device="cuda")
    y = torch.empty_like(x)
    y.copy_(x)                                       # warm up
    torch.cuda.synchronize()
    start = torch.cuda.Event(enable_timing=True)
    end = torch.cuda.Event(enable_timing=True)
    start.record()
    for _ in range(iters):
        y.copy_(x)
    end.record()
    torch.cuda.synchronize()
    ms = start.elapsed_time(end) / iters
    return 2 * n_bytes / (ms * 1e-3) / 1e12          # TB/s (read x + write y)

if gpu_report():
    bw = measured_hbm_bandwidth()
    print(f"measured HBM bandwidth  : {bw:.2f} TB/s")
    # Pair this with the measured BF16 matmul TFLOP/s from earlier in the chapter:
    # ridge = measured_TFLOPs * 1e12 / (bw * 1e12)  -> YOUR chip's FLOP/byte.

Two more habits worth forming now. First, divide measured throughput by peak to get model FLOPs utilization (MFU) — the standard headline efficiency number for a training run, formalized in The Roofline Model & Performance Engineering; well-tuned dense BF16 transformer pretraining typically lands in the tens of percent, and anything in the low single digits means you are host-bound, memory-bound, or not on the Tensor Cores at all. Second, remember that nvidia-smi’s memory number reflects PyTorch’s caching allocator, not live tensors: use torch.cuda.memory_allocated() / max_memory_allocated() for the truth, and see Memory-Efficient Training: Checkpointing, Offloading & LoRA Math for allocator tuning and fragmentation.

The standard NVIDIA tooling, in ascending order of detail:

# 1. What hardware am I on, and is it busy? (ships with the driver)
nvidia-smi --query-gpu=name,memory.total,memory.used,utilization.gpu,clocks.sm \
           --format=csv
nvidia-smi topo -m            # NVLink vs PCIe topology between GPUs in the node

# 2. Whole-program timeline: are there gaps? is the GPU idle waiting on Python?
nsys profile -o trace python train.py        # NVIDIA Nsight Systems

# 3. One kernel, in forensic detail: occupancy, achieved bandwidth, Tensor Core
#    utilization, and a measured roofline point for that kernel.
ncu --set full -k regex:my_kernel python train.py   # NVIDIA Nsight Compute

Inside PyTorch, torch.profiler.profile(activities=[ProfilerActivity.CPU, ProfilerActivity.CUDA]) gives you a per-op CUDA-time table without leaving Python, and is usually the right first stop; reach for Nsight Systems when you suspect launch overhead or host-bound stalls, and Nsight Compute when a single kernel is the problem (CUDA Programming Essentials for ML Engineers).

Putting It Together: Reading a Workload Like the Hardware Does¶

When you profile or design a kernel, walk the hierarchy top-down and ask the questions the hardware forces on you:

1. Is this op a matmul? If yes, it can use Tensor Cores — aim for peak by tiling through SMEM and feeding contiguous, well-shaped BF16/FP8 tiles. If no (pointwise, normalization, softmax, reductions), it runs on CUDA cores and is almost certainly memory-bound; your only lever is to fuse it with neighbors so intermediate tensors never round-trip to HBM.
2. What is its arithmetic intensity? Estimate FLOPs and HBM bytes. Compare to the chip’s ridge point. This single ratio tells you whether to optimize for compute (better tiling, higher Tensor Core utilization) or for memory (fewer bytes moved: fusion, recomputation, lower precision, better reuse).
3. Is the working set staying in fast tiers? A memory-bound kernel that re-reads the same HBM data many times is leaving enormous performance on the table. Stage reused data into SMEM/registers (the FlashAttention insight).
4. Are the warps full and busy? Check for warp divergence, uncoalesced loads, and shared-memory bank conflicts. Check occupancy — but only to confirm you have enough latency-hiding, not to maximize it blindly.
5. Is the GPU even the bottleneck? Profile end-to-end. Host-side Python overhead, tiny serialized kernels (launch latency), and PCIe transfers routinely cap throughput well below the GPU’s capability. CUDA graphs and torch.compile exist to kill launch overhead (Kernel Fusion, torch.compile, CUDA Graphs & Compilers).

import torch

def classify_op(flops, hbm_bytes, peak_flops=990e12, peak_bw=3.35e12):
    """A 30-second back-of-envelope to decide where to spend optimization effort."""
    intensity = flops / hbm_bytes
    ridge = peak_flops / peak_bw
    attainable = min(peak_flops, intensity * peak_bw)
    bound = "compute" if intensity > ridge else "memory"
    return dict(intensity=round(intensity, 2),
                ridge=round(ridge, 1),
                bound=bound,
                attainable_TFLOPs=round(attainable / 1e12, 1),
                advice=("tile/fuse-for-TC; you can chase peak FLOPs"
                        if bound == "compute"
                        else "cut bytes: fuse, recompute, lower precision, reuse"))

# Example: a fused LayerNorm over a (batch*seq, hidden) tensor, bf16.
rows, hidden = 4096 * 2048, 8192       # tokens x hidden
elems = rows * hidden
# LayerNorm: ~a handful of FLOPs/elem; reads + writes the tensor once each.
ln = classify_op(flops=8 * elems, hbm_bytes=2 * elems * 2)   # read+write, 2B/elem
print("LayerNorm:", ln)   # -> memory-bound: the reason norms get fused into matmuls

Key Takeaways¶

Further Reading¶

• NVIDIA, CUDA C++ Programming Guide and the per-architecture whitepapers (Ampere GA100, Hopper GH100, Blackwell) — The authoritative source for the execution model, memory hierarchy, and the per-generation specs summarized in this chapter. Read the architecture whitepapers for the “what changed” story.
• Williams, Waterman & Patterson, Roofline: An Insightful Visual Performance Model for Multicore Architectures (2009) — The original roofline paper that formalizes arithmetic intensity, the ridge point, and the compute- vs memory-bound dichotomy.
• Vasily Volkov, Better Performance at Lower Occupancy (GTC 2010) — The classic demonstration that maximizing occupancy is not the same as maximizing performance; essential context for the occupancy section.
• Dao et al., FlashAttention: Fast and Memory-Efficient Exact Attention with IO-Awareness (2022) and Dao, FlashAttention-2 (2023) — The canonical case studies in trading FLOPs for HBM traffic by keeping the working set in SMEM/registers; the practical payoff of everything in this chapter.
• Horace He, Making Deep Learning Go Brrrr From First Principles (2022) — A widely-cited engineer’s walkthrough of compute-bound vs memory-bound vs overhead-bound regimes in PyTorch, with the same intensity-first mindset.
• NVIDIA CUTLASS and OpenAI Triton repositories — Production and pedagogical implementations of tiled, Tensor-Core GEMM and fused kernels; the best way to see the SMEM-tiling and warp-level mechanics described here turned into real code (see Writing GPU Kernels with Triton).

Exercises¶

1. A kernel contains the branch if (x[i] > 0) { y[i] = expensive_A(x[i]); } else { y[i] = expensive_B(x[i]); }, where expensive_A and expensive_B each take about the same time \(t\) and neither is trivial. For a warp in which 16 lanes take the if side and 16 take the else side, how long does the warp take relative to the best case where all 32 lanes take the same side? Explain the mechanism, and describe one way to restructure the data so the warp stops paying this cost.

2. Consider an SM with these limits: 65,536 32-bit registers, 100 KB (take as 100,000 bytes) of shared memory, a cap of 2,048 resident threads, and a cap of 32 resident blocks. You launch a kernel with 128 threads per block, using 40 registers per thread and 8 KB (8,192 bytes) of shared memory per block. How many blocks fit on the SM, and what is the resulting occupancy?

3. A dense 34B-parameter model is served on a single H100 (peak \(\approx 990\) TFLOP/s BF16 dense, HBM bandwidth \(\beta \approx 3.35\) TB/s), using FP8 weights (1 byte per parameter). For single-request (batch-1) autoregressive decode: (a) how many bytes must be read from HBM to generate one token, and what is the bandwidth-floor per-token latency? (b) What is the arithmetic intensity, and is decode compute- or memory-bound? © The H100 ridge point in FP8 is about \(1979 / 3.35 \approx 590\) FLOP/byte. Roughly what minimum batch size would push batch-scaled decode past that ridge?

4. Using the chapter’s A100 numbers (\(\pi \approx 312\) TFLOP/s BF16, \(\beta \approx 2.0\) TB/s, ridge \(\approx 156\) FLOP/byte), consider a 7B model in BF16. (a) At batch 1, the arithmetic intensity of decode is \(I = 1\) FLOP/byte. What is the smallest integer batch size at which decode becomes compute-bound? (b) At that batch size, is per-token latency still floored by the same 7 ms bandwidth number as batch 1? Explain what has (and has not) changed.

5. Implement a small helper, in the style of the chapter’s roofline code, that answers the two questions of batched decode at once: decode_analysis(params, bytes_per_param, batch, peak_flops, peak_bw) should return the arithmetic intensity, whether decode is compute- or memory-bound, and the per-step time under the roofline (the max of the bandwidth time and the compute time). Then add min_batch_for_compute_bound(bytes_per_param, peak_flops, peak_bw) returning the smallest integer batch that crosses the ridge. Run both for a 7B BF16 model on the A100 numbers from Exercise 4 and confirm the crossover matches your hand calculation.

import math

def decode_analysis(params, bytes_per_param, batch, peak_flops, peak_bw):
    bytes_moved = params * bytes_per_param        # weights read once per step
    flops       = 2 * params * batch              # ~2 FLOPs/param/token
    intensity   = flops / bytes_moved
    ridge       = peak_flops / peak_bw
    bw_time      = bytes_moved / peak_bw           # seconds to stream weights
    compute_time = flops / peak_flops              # seconds to do the math
    step_time    = max(bw_time, compute_time)      # roofline: the binding one
    return dict(
        intensity=round(intensity, 2),
        ridge=round(ridge, 1),
        bound=("compute" if intensity > ridge else "memory"),
        step_ms=round(step_time * 1e3, 3),
        per_token_ms=round(step_time * 1e3 / batch, 4),
    )

def min_batch_for_compute_bound(bytes_per_param, peak_flops, peak_bw):
    ridge = peak_flops / peak_bw
    # intensity(B) = (2 / bytes_per_param) * B ; want intensity > ridge
    b = ridge * bytes_per_param / 2.0
    return math.floor(b) + 1                       # first strictly-above integer

# 7B BF16 on A100 numbers from Exercise 4.
PARAMS   = 7e9
BPP      = 2                 # BF16
PEAK_FLOPS = 312e12          # 312 TFLOP/s
PEAK_BW    = 2.0e12          # 2.0 TB/s

bmin = min_batch_for_compute_bound(BPP, PEAK_FLOPS, PEAK_BW)
print("min batch for compute-bound:", bmin)        # -> 157

for B in (1, bmin - 1, bmin):
    print(B, decode_analysis(PARAMS, BPP, B, PEAK_FLOPS, PEAK_BW))

min batch for compute-bound: 157
1   {'intensity': 1.0,   'ridge': 156.0, 'bound': 'memory',  'step_ms': 7.0,    'per_token_ms': 7.0}
156 {'intensity': 156.0, 'ridge': 156.0, 'bound': 'memory',  'step_ms': 7.0,    'per_token_ms': 0.0449}
157 {'intensity': 157.0, 'ridge': 156.0, 'bound': 'compute', 'step_ms': 7.045,  'per_token_ms': 0.0449}