10. How GPUs run your code, and how to measure it

Why this part of the book exists

Your GPT works. Chapter 1 promised a fast engine, and “fast” means knowing what your hardware can do and making your code do it. The next six chapters teach GPU programming from the ground up: the execution model (this chapter), CUDA (11), matrix multiplication (12), reductions and floating point (13), Triton (14) and FlashAttention (15). Each technique gets used in the engine in Part IV.

Two philosophies of processor design

A CPU core is built to finish one stream of instructions as quickly as possible. Large caches, branch prediction and out-of-order execution all minimize the latency of each instruction. A desktop CPU has 8-32 of these sophisticated cores.

A GPU is built to finish an enormous amount of independent work, without caring how long any single piece takes. It has thousands of simple arithmetic units, small caches per unit, and no speculation. When one group of threads waits for memory, which takes hundreds of cycles, the hardware instantly switches to another group that’s ready. As long as there’s enough independent work, the arithmetic units stay busy. This is throughput-oriented design, and latency hiding through parallelism is the core idea. PMPP Chapter 1 frames the whole field this way.

Matrix multiplication across thousands of tokens and features is about as independent as work gets, which is why deep learning runs on GPUs.

Inside a GPU

An NVIDIA GPU is a collection of streaming multiprocessors (SMs): 132 on an H100, 128 on an RTX 4090. Each SM contains:

• Arithmetic units: FP32/INT32 cores, special-function units for exp and sin, and tensor cores, which perform small matrix multiplications (for example, 16×8×16 in BF16) in a single instruction.
• A large register file (256 KB per SM): the fastest storage, private to each thread.
• Shared memory / L1 cache (up to about 228 KB per SM on Hopper): fast, explicitly managed, shared among the threads of a block.
• Warp schedulers that pick, every cycle, which group of threads issues next.

All SMs share an L2 cache (tens of MB) and the main device memory (HBM or GDDR, tens of GB). Approximate figures for three machines used as examples in this book:

(Vendor datasheet figures; the dense BF16 numbers exclude structured sparsity.)

Threads execute in groups of 32 called warps. All 32 threads of a warp execute the same instruction at the same time on different data, a model NVIDIA calls SIMT (single instruction, multiple threads). Chapter 11 shows how you organize threads into blocks and grids. For now, the takeaway is the memory hierarchy: registers and shared memory are tiny but fast, device memory is huge but relatively slow, and making sure that data loaded from device memory is reused many times before being evicted is most of what “optimizing a kernel” means.

Your Python program submits work; the GPU does it later

When PyTorch runs y = x @ w on a CUDA tensor, the CPU doesn’t multiply anything. It enqueues a kernel on a stream (an ordered queue of GPU work) and returns immediately, usually within 5-20 µs. The GPU executes queued kernels in order, while the CPU races ahead submitting more.

This asynchrony is essential for performance: the GPU never waits for Python if the CPU keeps the queue full. But it has two consequences you must internalize:

1. Anything that needs a GPU value on the CPU waits for the queue to drain. print(y), y.item(), y.tolist(), if y > 0:, .cpu() and .numpy() all synchronize. In a decode loop, one .item() per token means the CPU stops, waits for the GPU, and only then submits the next step’s kernels, leaving the GPU idle while Python runs. Chapter 19 removes these.
2. Naive timing measures submission, not work.

Measuring time correctly

This looks reasonable and is wrong:

start = time.perf_counter()
y = x @ x                     # only *enqueued*
elapsed = time.perf_counter() - start   # ~10 µs, regardless of the matrix size

Two correct options:

• Synchronized wall time. Call torch.cuda.synchronize() before starting the clock (so earlier queued work isn’t counted) and after the operation (so this work is). This measures what a user experiences, including CPU overheads.
• CUDA events. Record an event before and after the work on the GPU’s own timeline and ask for the elapsed time between them. This excludes CPU-side gaps, so it measures kernel time.

def measure_wall(function, samples=10, warmup=3, cuda=None):
    """Median wall-clock milliseconds of function(), synchronizing the GPU before starting and
    after finishing each sample, so queued work is neither excluded nor borrowed.  (Your engine: Chapter 10)"""
    cuda = torch.cuda.is_available() if cuda is None else cuda
    for _ in range(warmup):
        function()
    times = []
    for _ in range(samples):
        if cuda:
            torch.cuda.synchronize()
        start = time.perf_counter()
        function()
        if cuda:
            torch.cuda.synchronize()
        times.append((time.perf_counter() - start) * 1000)
    return {"median_ms": statistics.median(times), "min_ms": min(times), "max_ms": max(times), "samples": samples}


def measure_cuda(function, warmup=5, samples=20, repeats=1):
    """GPU-side milliseconds between two CUDA events recorded on the current stream."""
    for _ in range(warmup):
        function()
    torch.cuda.synchronize()
    times = []
    for _ in range(samples):
        start, end = torch.cuda.Event(enable_timing=True), torch.cuda.Event(enable_timing=True)
        start.record()
        for _ in range(repeats):
            function()
        end.record()
        end.synchronize()
        times.append(start.elapsed_time(end) / repeats)
    return {"median_ms": statistics.median(times), "min_ms": min(times), "max_ms": max(times), "samples": samples}

Three more rules make measurements trustworthy:

• Warm up. The first call pays one-time costs: CUDA context creation, library loading, kernel compilation (Triton, torch.compile) and memory-pool growth. Run a few iterations untimed.
• Report a distribution. Take the median of many samples and keep the min and max. GPU clocks, other processes and thermal throttling all add noise.
• Fix everything except the variable you’re testing. Shapes, dtype, device and input values (some kernels are faster on zeros!) must stay the same across the comparison.

Profiling: why is it slow?

A timer tells you how long; a profiler tells you why. PyTorch’s profiler records every operator and kernel, CPU side and GPU side:

python run.py profile          # writes runs/trace.json and prints the top operators

Open the trace in Perfetto (or chrome://tracing). You’ll see two timelines: the CPU thread submitting operators, and the GPU stream executing kernels. Look for:

• Gaps on the GPU timeline: the GPU is idle waiting for the CPU. The cause is Python overhead, synchronization, or too many tiny kernels. Fixes: fusion, CUDA graphs, batching (Chapters 14, 19, 24).
• Unexpected kernels: a copy_ or contiguous you didn’t mean to cause, or a dtype conversion.
• One kernel dominating: optimize that, and check with the roofline whether it’s already near its limit.

For kernel-level detail, NVIDIA’s Nsight Systems (nsys) shows the whole-program timeline, and Nsight Compute (ncu) shows one kernel’s achieved bandwidth, occupancy and stall reasons. GPU Mode L1 walks through all three tools.

Tip

Profilers add overhead. Use a profile to explain, and a separate un-profiled run to measure.

The roofline model

Before optimizing anything, ask what limits it. Every kernel does some arithmetic (FLOPs) and moves some bytes to and from device memory. Its arithmetic intensity is the ratio:

$$ I = \frac{\text{FLOPs}}{\text{bytes moved}} . $$

A device can do at most $F$ FLOP/s and move at most $W$ bytes/s. So the best achievable throughput is

$$ \text{attainable FLOP/s} = \min(F,; I \cdot W). $$

Plot that against $I$ and you get a “roofline”: a slanted line (memory-bound) that meets a flat line (compute-bound) at the ridge point $I^* = F/W$. For an H100, $I^* \approx 990{,}000/3{,}350 \approx 295$ FLOPs per byte. An operation with lower intensity can’t use the tensor cores fully, no matter how good the kernel is.

Now place real operations on it:

• Vector add in FP32 reads 8 bytes and writes 4 per element, and does 1 FLOP: $I = 1/12$. Hopelessly memory-bound. All you can do is approach the bandwidth limit.
• Matrix × vector (one token of decode through a [4096, 4096] BF16 weight) does $2 \cdot 4096^2$ FLOPs and reads the $2 \cdot 4096^2$-byte matrix: $I \approx 1$. Memory-bound by a factor of about 300 on an H100. This is the Chapter 1 result again, now with a name.
• Batched decode with $B$ tokens per weight read: $I \approx B$. That’s 1 at batch 1, 8 at batch 8, about 120 at batch 128, and about 410 at batch 512 (the input and output reads start to count). Batching is how decode climbs the roofline (Chapter 24).
• Prefill of a 4,096-token prompt through the same weight is a [4096, 4096] @ [4096, 4096] GEMM: $I \approx 1{,}365$, comfortably compute-bound.

This is the most important analytical tool in the book. Before writing any kernel, compute its intensity and its ceiling, so you know when to stop optimizing.

The milestone turns these estimates into functions you can call from any experiment:

def matmul_intensity(m, n, k, bytes_per_element=2):
    """FLOPs per byte of [M,K] @ [K,N] if every operand is read once and C written once.  (Your engine: Chapter 10)"""
    flops = 2 * m * n * k
    traffic = (m * k + k * n + m * n) * bytes_per_element
    return flops / traffic


def attainable_flops(intensity, peak_flops, bandwidth):
    """The roofline: you cannot exceed the compute peak, nor bandwidth x intensity.  (Your engine: Chapter 10)"""
    return min(peak_flops, intensity * bandwidth)


def decode_ceiling(weight_bytes, bandwidth, batch=1, kv_bytes_per_sequence=0):
    """Upper bound on decode tokens/s when every step must read all weights once (shared by the
    batch) plus each sequence's own KV cache.  (Your engine: Chapter 10)"""
    step_bytes = weight_bytes + batch * kv_bytes_per_sequence
    return batch * bandwidth / step_bytes

Amdahl’s law: optimize what matters

If an operation takes a fraction $f$ of the total time and you make it $s$ times faster, the whole program gets

$$ \text{speedup} = \frac{1}{(1-f) + f/s} $$

times faster. Make a 10% operation infinitely fast and the program speeds up by only 1.11x. That’s why you profile the whole engine first (Chapter 19), and then optimize the biggest bar, not the most interesting kernel.

Important

DGX Spark. CPU and GPU share 128 GB of LPDDR5x. Two consequences: the 273 GB/s bandwidth is shared too (a CPU process streaming memory slows your GPU decode), and Linux’s “free” memory figure is misleading. Watch “available” in free -h and use torch.cuda.mem_get_info(). Moving data between cpu and cuda tensors still copies, even though the bytes sit in the same physical memory.

Build it

pytest tests/test_ch10_performance.py
python run.py profile

Stretch exercises

1. ★ Time x @ x for a 4096² FP32 matrix three ways: unsynchronized wall clock, synchronized wall clock, and CUDA events. Explain each number. Where: experiments/ch10.py (create it), using engine.measure.measure_wall and measure_cuda.
2. ★★ Measure the achieved bandwidth of x + y for vector sizes from $2^{10}$ to $2^{28}$ elements. Plot GB/s against size. Where does launch overhead dominate, and what fraction of datasheet bandwidth do you reach at the top? Where: experiments/ch10.py (create it).
3. ★★ Measure the TFLOP/s of a @ b in BF16 for square sizes 256 to 8,192 and plot them on your device’s roofline. At what size do you reach 70% of peak? Where: experiments/ch10.py (create it), using engine.measure for timing and roofline arithmetic.
4. ★★★ Profile one decode step of your GPT-2 with torch.profiler. Count the kernels, and estimate the fraction of step time spent in launch gaps rather than kernels. Where: experiments/ch10.py (create it), wrapping a decode step in torch.profiler.profile.

Check your understanding

1. Why must you synchronize before starting a timer as well as after the operation?
2. Why should profiling and benchmarking be separate runs?
3. What is the arithmetic intensity of decode for a batch of 16 sequences, and is it memory-bound on an H100?
4. A kernel taking 30% of runtime becomes 3x faster. What’s the overall speedup?
5. On DGX Spark, why can a CPU-heavy job slow down GPU decoding?

Going deeper

• PMPP Chapter 1 (heterogeneous computing, latency vs throughput), Chapter 4 §§4.1-4.7 (GPU architecture, warps, scheduling, occupancy), and §22.5 (batching: latency versus throughput).
• GPU Mode L1 (Mark Saroufim, profiling and integrating kernels in PyTorch, with nsys/ncu examples) and L8 (the CUDA performance checklist).
• Williams, Waterman and Patterson, Roofline: An Insightful Visual Performance Model (2009).
• NVIDIA Hopper and Ada architecture whitepapers for exact SM resources.