GPU architecture tour

Track 07 - GPU programming · Notebook 01 · Runtime: ≈10 min on T4

Prerequisites (undergraduate CS): you’ve seen a computer-architecture course (CPU pipeline, cache hierarchy, memory bus) and you’ve written a tight loop in Python or C. No GPU experience required.

What you’ll know by the end: you can walk up to any NVIDIA GPU and say “this chip has X compute power, Y memory bandwidth, and these are the workloads where it’s bottlenecked on memory vs on math.” You’ll also have a mental model of why LLMs spend most of their time moving bytes around rather than computing.

---

Why start with hardware?

Every optimisation in the rest of this curriculum - KV caches, paged attention, quantisation, kernel fusion, speculative decoding - is really just a clever way of moving fewer bytes. If you don’t have a feel for why bytes are expensive, those techniques feel like magic. They’re not magic.

So: one notebook, one chip, before we touch language models.

A CS-major’s mental model of a GPU

You already know how a CPU works. A GPU is the same idea pushed to a different operating point:

Dimension	CPU (e.g. laptop)	GPU (e.g. T4)
Cores	4-16 wide, complex	thousands of tiny cores grouped into SMs (streaming multiprocessors)
What each core does	branch, speculate, out-of-order	run the same instruction on 32 threads at once (SIMT, like SIMD with more flexibility)
Main memory	DDR4/5, ~50 GB/s	HBM, ~300 GB/s to 3 TB/s
Cache hierarchy	L1 → L2 → L3 → DRAM	L1 (per-SM) → L2 (shared) → HBM
Cores are good at	unpredictable branchy code	massive, regular matrix-like work

Two numbers describe any GPU for our purposes:

1. Peak compute - how many FLOPs/second the math units can retire.

2. Peak memory bandwidth - how many bytes/second you can drag from HBM into the compute units.

That’s it. Everything else in this notebook - and in practice most GPU performance analysis - is about comparing your kernel’s measured throughput against those two ceilings.

The single most important insight

If your kernel reads more bytes from HBM than it does math on them, you will be bandwidth-bound. It won’t matter how fast the compute is.

This is why modern GPUs have so many TFLOPs: not because we need them for general-purpose code, but because the math is cheap provided you can feed it. The hard part - and the reason a 100 B LLM on an H100 isn’t 10 000× faster than on a phone - is getting the weights and activations into the compute units in time.

LLM inference is the textbook example. Generating one token requires reading every weight of the model from HBM. An 8 B model in fp16 is 16 GB of reads; at 3 TB/s that’s ≈5 ms of pure data movement per token. You cannot decode faster than that, no matter how much compute you throw at it.

We’re going to measure both ceilings on whatever GPU you’re running this notebook on, and then plot the roofline (Williams, Waterman, Patterson 2009) diagram that makes the tradeoff visible in one picture.

from llm_systems_cookbook.nb import bootstrap

import time

import numpy as np
import torch

s = bootstrap("07_gpu_01_gpu_architecture_tour")

Step 1 - Read the datasheet

For each GPU family we hard-code two numbers: advertised peak memory bandwidth (GB/s) and advertised peak FP16 throughput (TFLOPs). These come straight from the NVIDIA product brief PDFs. We’ll compare what we measure against these ceilings - which tells us how close the kernel is to the hardware limit.

An important subtlety: “peak” numbers are marketing ceilings. Real kernels hit 60-90 % of peak on a good day; if you measure 40 % it’s probably fine, if you measure 10 % it’s time to investigate.

# Approximate datasheet peaks per device family.
# FP16 throughput is tensor-core dense (no sparsity). BW is HBM.
# H100 row uses SXM5 numbers; H100 PCIe is ~2000 GB/s and ~756 TFLOPs.
DEVICE_PEAKS: dict[str, dict[str, float]] = {
    "T4":            {"peak_bw_gbps":  320.0, "peak_fp16_tflops":  65.0},
    "L4":            {"peak_bw_gbps":  300.0, "peak_fp16_tflops": 120.0},
    "A10G":          {"peak_bw_gbps":  600.0, "peak_fp16_tflops": 125.0},
    "A100":          {"peak_bw_gbps": 2039.0, "peak_fp16_tflops": 312.0},
    "H100":          {"peak_bw_gbps": 3350.0, "peak_fp16_tflops": 989.0},
    "V100":          {"peak_bw_gbps":  900.0, "peak_fp16_tflops": 125.0},
    "RTX 3090":      {"peak_bw_gbps":  936.0, "peak_fp16_tflops":  71.0},
    "RTX 4090":      {"peak_bw_gbps": 1008.0, "peak_fp16_tflops": 165.0},
}


def match_device(name: str) -> dict[str, float] | None:
    for key, peaks in DEVICE_PEAKS.items():
        if key.lower() in name.lower():
            return peaks
    return None


if IS_CUDA:
    dev_name = torch.cuda.get_device_name(0)
    peaks = match_device(dev_name) or DEVICE_PEAKS["A10G"]
    print(f"detected: {dev_name}")
    print(f"peaks:    {peaks}")
else:
    peaks = None
    print("no CUDA device; the roofline math still runs but the microbenchmarks are skipped.")

Step 2 - Measure bandwidth

The simplest possible way to measure “how fast can this chip move bytes”:

allocate two big buffers in HBM
for N iterations:
    copy buffer A → buffer B
bandwidth = (2 × bytes_per_copy × N) / wall_time

The factor of 2 is because each copy is one read and one write - both consume HBM bandwidth. A 500 MB copy loop should saturate the memory bus and give us a number close to the datasheet peak.

Why this works: there’s nothing for the compute units to do (we’re not transforming the data), so the kernel becomes a pure memory-bandwidth test. If the measured number is wildly below peak, the answer is almost always “you didn’t move enough bytes per launch to amortise kernel-launch overhead” - not that the chip is broken.

def measure_bandwidth_gbps(n_bytes: int = 512 * 1024 * 1024, n_iter: int = 20) -> float:
    '''Sustained HBM bandwidth from a contiguous float32 copy (in GB/s).'''
    if not torch.cuda.is_available():
        return 0.0
    n = n_bytes // 4
    a = torch.empty(n, device="cuda", dtype=torch.float32).uniform_()
    b = torch.empty_like(a)
    # Warm-up: first few launches pay for CUDA context init + code cache.
    for _ in range(3):
        b.copy_(a)
    torch.cuda.synchronize()
    t0 = time.perf_counter()
    for _ in range(n_iter):
        b.copy_(a)
    torch.cuda.synchronize()
    return 2 * n_bytes * n_iter / (time.perf_counter() - t0) / 1e9


bw_gbps = measure_bandwidth_gbps() if IS_CUDA else 0.0
if IS_CUDA:
    frac = bw_gbps / peaks["peak_bw_gbps"]
    print(f"measured bandwidth: {bw_gbps:.1f} GB/s  ({frac:.0%} of advertised peak)")
    s.benchmark(
        "bandwidth_at_least_40pct_peak",
        measured=bw_gbps,
        baseline=peaks["peak_bw_gbps"],
        must_beat=0.40,
        unit="GB/s",
    )
else:
    s.skip("bandwidth_at_least_40pct_peak", "no CUDA device available")

Step 3 - Measure compute

Now we do the opposite - arrange a workload that’s so arithmetically intense it can’t possibly be bandwidth-bound, and see how close we get to the compute ceiling.

A square fp16 matrix multiply is the go-to test. For A @ B where both are \(m \times m\):

• FLOPs: \(2m^3\) (every output element is a length-\(m\) dot product, and each dot product is \(m\) multiplies + \(m\) adds).

• Bytes: \(\approx 3 m^2 \cdot 2\) (two input matrices and one output, in fp16).

• Ratio (arithmetic intensity): \(2m^3 / (6 m^2) = m/3\).

For \(m = 4096\) that’s about 1365 FLOPs per byte - so much math per byte that only the compute units matter. Exactly the workload that sells GPUs.

def measure_matmul_tflops(m: int = 4096, n_iter: int = 30, dtype: torch.dtype = torch.float16) -> float:
    '''Sustained TFLOPs from square fp16 matmul of side ``m``.'''
    if not torch.cuda.is_available():
        return 0.0
    a = torch.randn((m, m), device="cuda", dtype=dtype)
    b = torch.randn((m, m), device="cuda", dtype=dtype)
    for _ in range(3):
        torch.matmul(a, b)
    torch.cuda.synchronize()
    t0 = time.perf_counter()
    for _ in range(n_iter):
        torch.matmul(a, b)
    torch.cuda.synchronize()
    return 2 * m**3 * n_iter / (time.perf_counter() - t0) / 1e12


tflops = measure_matmul_tflops() if IS_CUDA else 0.0
if IS_CUDA:
    frac = tflops / peaks["peak_fp16_tflops"]
    print(f"measured FP16 matmul: {tflops:.1f} TFLOPs  ({frac:.0%} of advertised peak)")
    s.benchmark(
        "matmul_tflops_at_least_40pct_peak",
        measured=tflops,
        baseline=peaks["peak_fp16_tflops"],
        must_beat=0.40,
        unit="TFLOPs",
    )
else:
    s.skip("matmul_tflops_at_least_40pct_peak", "no CUDA device available")

Both ceilings side by side

Before we draw the full roofline, eyeball the two measurements against their datasheet peaks. 60-90% of peak is healthy; below 40% usually means the microbenchmark wasn’t big enough to amortise kernel-launch overhead, not that the chip is broken.

import matplotlib.pyplot as plt

if IS_CUDA:
    labels   = ["bandwidth (GB/s)", "compute (TFLOPs fp16)"]
    measured = [bw_gbps,             tflops]
    peak     = [peaks["peak_bw_gbps"], peaks["peak_fp16_tflops"]]

    fig, axes = plt.subplots(1, 2, figsize=(8, 3.2))
    for ax, name, m, p in zip(axes, labels, measured, peak, strict=True):
        ax.bar(["measured", "datasheet peak"], [m, p], color=["tab:blue", "tab:gray"])
        ax.set_title(f"{name}\n{m:.1f} / {p:.0f}  =  {m / p:.0%} of peak")
        ax.set_ylim(0, max(m, p) * 1.15)
    fig.suptitle(f"how close did we get? ({torch.cuda.get_device_name(0)})")
    fig.tight_layout()
    plt.show()
else:
    print("skipped - no CUDA device, so no measured numbers to plot against peak.")

Step 4 - The roofline diagram

We now have two ceilings: a bandwidth ceiling (\(\approx\) 320 GB/s on T4) and a compute ceiling (\(\approx\) 65 TFLOPs on T4). Every kernel on that chip lives under the lower of those two limits, depending on how much math it does per byte it reads.

The roofline plot draws both ceilings on one log-log chart:

• x-axis: arithmetic intensity (FLOPs per byte of data)

• y-axis: achievable throughput (TFLOPs)

• the “roof” is min(bandwidth × intensity, peak_compute)

The knee of the roof is at peak_compute / peak_bandwidth - call it the ridge intensity. Below the ridge you’re memory-bound; above it you’re compute-bound.

One-line glossary: prefill and decode

LLM inference splits into two phases that land on opposite sides of the ridge. They’ll come up constantly from here on, so pin them down now:

• Prefill - the model reads the whole input prompt once and computes hidden states (and the initial KV cache) for every token in parallel. Each weight read is amortised across many tokens of math, so prefill is usually compute-bound.

• Decode - the model then emits new tokens one at a time. Each step is a single query row attending over the cached keys/values, but it still has to read every weight from HBM. Tiny math-to-bytes ratio, so decode is almost always memory-bound.

On T4, ridge \(\approx\) 65 TFLOPs / 320 GB/s \(\approx\) 200 FLOPs/byte. Decode on a 7 B LLM is \(\approx\) 1-30 FLOPs/byte - deep in the memory-bound region, which is why decode is slow. Prefill on a long prompt is 700-2000 FLOPs/byte - squarely compute-bound.

import matplotlib.pyplot as plt


def matmul_intensity(m: int, k: int, n: int, dtype_bytes: int = 2) -> float:
    '''FLOPs per byte for ``A @ B`` with shapes (m,k) and (k,n).'''
    flops = 2 * m * k * n
    bytes_ = dtype_bytes * (m * k + k * n + m * n)
    return flops / bytes_


SHAPES = [
    (1,      4096, 4096),   # vector-matrix: tiny intensity, bandwidth-bound
    (16,     4096, 4096),
    (128,    4096, 4096),
    (1024,   4096, 4096),
    (4096,   4096, 4096),   # big square: huge intensity, compute-bound
]
intensities = [matmul_intensity(m, k, n) for (m, k, n) in SHAPES]
for shape, ai in zip(SHAPES, intensities, strict=True):
    print(f"  {shape[0]:>4}x{shape[1]}x{shape[2]}  AI = {ai:7.1f} FLOPs/byte")

peak_bw = (peaks or {"peak_bw_gbps": 320.0})["peak_bw_gbps"]
peak_compute = (peaks or {"peak_fp16_tflops": 65.0})["peak_fp16_tflops"]
ridge = peak_compute * 1e12 / (peak_bw * 1e9)
print(f"\nridge intensity = {ridge:.1f} FLOPs/byte")

fig, ax = plt.subplots(figsize=(7, 4))
x = np.logspace(-1, 3, 200)
y_mem = peak_bw * 1e9 * x / 1e12
y_comp = np.full_like(x, peak_compute)
ax.loglog(x, np.minimum(y_mem, y_comp), label=f"roof (bw={peak_bw:.0f} GB/s, peak={peak_compute:.0f} TFLOPs)")
ax.axvline(ridge, color="tab:gray", linestyle=":", label=f"ridge @ {ridge:.0f}")
achieved = [peak_compute * min(1.0, ai / ridge) for ai in intensities]
ax.scatter(intensities, achieved, color="tab:red", zorder=5, label="matmul shapes")
ax.set_xlabel("arithmetic intensity (FLOPs / byte)")
ax.set_ylabel("throughput (TFLOPs)")
ax.set_title("roofline: bandwidth-bound left, compute-bound right")
ax.legend()
ax.grid(True, which="both", alpha=0.3)
plt.show()

What the checks below are asking

1. Small-batch matmul is memory-bound. When you multiply a single row by a 4096×4096 matrix, there’s barely any math per byte read. You’ll never hit peak compute, no matter how fast the chip is.

2. Large square matmul is compute-bound. At 4096×4096×4096 the ratio goes into the thousands and the chip finally gets to stretch its legs.

3. Intensity grows with batch size. This is the single most important lever for inference economics: batching amortises the weight read across many queries.

If all three hold, you’ve captured the core of GPU performance analysis in three predicates.

s.check(
    "small_batch_is_memory_bound",
    lambda: matmul_intensity(*SHAPES[0]) < ridge,
    msg=f"intensity={matmul_intensity(*SHAPES[0]):.1f}  ridge={ridge:.1f}",
)
s.check(
    "large_square_is_compute_bound",
    lambda: matmul_intensity(*SHAPES[-1]) > ridge,
    msg=f"intensity={matmul_intensity(*SHAPES[-1]):.1f}  ridge={ridge:.1f}",
)
s.check(
    "intensity_monotone_in_batch",
    lambda: all(intensities[i] <= intensities[i + 1] for i in range(len(intensities) - 1)),
    msg="arithmetic intensity should grow with batch at fixed k,n",
)

Step 5 - Who’s this chip?

Before moving on, note down the compute capability and SM count of whatever chip you’re on. These come up again in every later GPU notebook:

• Compute capability (e.g. 7.5 on T4, 8.6 on A10G, 9.0 on H100) gates feature availability. FlashAttention-2 needs ≥ 8.0; FP8 tensor cores need ≥ 8.9.

• SM count is the parallelism budget. A T4 has 40 SMs; an H100 has 132. Triton kernels will typically launch one program instance per SM.

if IS_CUDA:
    props = torch.cuda.get_device_properties(0)
    cc = torch.cuda.get_device_capability(0)
    print(f"compute capability: {cc[0]}.{cc[1]}")
    print(f"SM count:           {props.multi_processor_count}")
    print(f"total HBM:          {props.total_memory / 1024**3:.1f} GiB")
    s.check("compute_capability_detected", lambda: cc[0] >= 5, msg=f"cc={cc}")
    s.check("nontrivial_sm_count", lambda: props.multi_processor_count >= 20,
            msg=f"{props.multi_processor_count} SMs")
else:
    s.skip("compute_capability_detected", "no CUDA device available")
    s.skip("nontrivial_sm_count", "no CUDA device available")

Exercises

1. Change the dtype. Rerun measure_matmul_tflops with torch.float32 instead of torch.float16. On T4 you should see roughly 15× lower throughput because fp32 doesn’t use the tensor cores. What does this tell you about how much of the “T4 is 65 TFLOPs” number is really “T4 is 65 TFLOPs only if you use tensor cores”?

2. Find the ridge for your GPU. Plot peak_compute / peak_bandwidth for every GPU in the table. How does the ridge change from T4 to A100 to H100? What does a higher ridge mean for LLM inference economics?

3. Simulate LLM decode. Compute the arithmetic intensity of a decode step on Llama-3-8B at batch size 1 (weights: 16 GB fp16, compute: 2 × 8e9 FLOPs per token). Is it above or below your GPU’s ridge?

Further reading

• Williams, Waterman, Patterson 2009, Roofline: An Insightful Visual Performance Model for Multicore Architectures - the original paper, still readable in an evening.

• NVIDIA product briefs (one PDF per GPU) - authoritative source for the datasheet peaks.

• Horace He 2022, Making Deep Learning Go Brrrr From First Principles - the modern LLM-flavoured version of the roofline story.

What’s next

With this mental model in hand you’re ready for Triton basics (02_triton_101_softmax), where we write our first kernel and measure how close it gets to the roofline we just derived.

s.summary()
s.save()