The GPU, for people who only ever call .cuda()

The hardware every other page assumes

Not a paper. The model underneath four other pages on this site, so they stop being four separate mysteries.

Download gpu_from_scratch.py — everything below runs from it. NumPy only, no GPU required.

I wrote the FlashAttention page and said SRAM is small and HBM is far. I wrote the KV cache page and said decode is memory-bound. I wrote the DeepSeek-V3 page and said the interconnect is the constraint. Three pages, three hand-waves at the same machine, and I never once explained the machine.

This is that page. There is no paper to summarise — but there is one 2009 paper that organises all of it, and I have built the whole thing as arithmetic you can run. The claim: you can predict which of your operations will be slow before you profile anything, from one ratio.

Read this before anything else on the page. Every hardware number here is a published specification — a vendor datasheet figure, rounded, approximate, and named to the part it belongs to. Nothing was measured. Peak numbers are ceilings nobody reaches; real kernels land well below them. This page tells you which wall you are standing against, not how fast your code will run.

A memory hierarchy is a supply chain

You already reason like this about a database. Some data is in the process, some in Redis, some on a disk in the same rack, some in another region. You never think "this query is slow because the CPU is slow". You think "how many round trips did I just buy?" A GPU is the same shape:

LevelPlain meaningSupply chain
RegistersScratch values a single thread is holding right nowIn your hands
Shared memory / L1A scratchpad the programmer manages, private to one SMThe shelf behind you
L2 cacheShared by the whole die, automatic, not yours to controlThe stockroom
HBMThe GPU's main memory. What .cuda() puts things inWarehouse across the yard
NVLinkThe link to another GPU in the same serverThe next building
InfiniBandThe link to a GPU in another serverAnother site entirely

Three NVIDIA words, defined once. A streaming multiprocessor (SM) is one of the ~130 independent processors on the die; you get the whole GPU only if your work fills all of them. A warp is 32 threads executing the same instruction in lockstep. Shared memory is a software-managed scratchpad on each SM — not a cache you hope hits, a buffer you explicitly load into. HBM is the stacked DRAM on the package: the 80 GB number in the marketing.

Don't take my word for it. Here's the actual output:

Accelerator: NVIDIA H100 SXM (published specification, approximate)
All figures PUBLISHED SPECIFICATIONS, not measurements.

level           capacity  bandwidth  latency* supply chain           FLOPs missed
---------------------------------------------------------------------------------
registers         256 KB  ~100 TB/s      1 ns on the bench                 990000
shared mem/L1     228 KB   ~20 TB/s     30 ns shelf behind you           29700000
L2 cache           50 MB   ~10 TB/s    200 ns stockroom                 198000000
HBM3               80 GB  3.35 TB/s    500 ns warehouse in the yard     495000000
NVLink peer            -   900 GB/s   2000 ns the next building        1980000000
InfiniBand             -    50 GB/s   5000 ns another site             4950000000

* latency is an order-of-magnitude figure, NOT a vendor spec.

The last column is the whole page in one number: peak FP16 throughput times the latency of that trip — the arithmetic you could have done while waiting. Half a billion floating point operations, for one round trip to the GPU's own main memory.

And the capacities run the other way: 228 KB of shared memory per SM against 80 GB of HBM. That asymmetry is not a detail — it is the entire reason FlashAttention exists, and why its tile size is what it is.

The chip is not slow at arithmetic. It is starving. Nearly every optimisation you will ever apply to a GPU is a way of taking fewer trips, or taking them further from the hot loop. Almost none of them are about doing less maths.

The one formula

Williams, Waterman and Patterson published this in 2009 — Roofline: An Insightful Visual Performance Model for Multicore Architectures, in Communications of the ACM. It predates CUDA's rise, predates tensor cores, predates every accelerator on this page, and still explains all of them, because it assumes only two things: you have a ceiling on arithmetic and a ceiling on bandwidth.

attainable FLOP/s = min( peak FLOP/s ,  bandwidth × I )    where I = FLOPs ÷ bytes moved

Read it as pipeline stages:

StepUnitsWhat it does
IFLOP/byteArithmetic intensity. A property of your operation, not the chip. How much work you get out of each byte you dragged in
bandwidth × IFLOP/sThe slanted roof. If you can only move B bytes a second and you do I flops per byte, that is your ceiling. No arithmetic unit can help
peak FLOP/sFLOP/sThe flat roof. The chip's published maximum. You cannot exceed it however many bytes you feed it
min(…)FLOP/sWhichever roof is lower is the one you hit. And which one it is tells you what to go and fix

That is it. Two straight lines and a min(). There is no curve to fit and nothing to measure. In code:

def roofline(intensity, peak_flops, bandwidth):
    memory_bound_ceiling = bandwidth * intensity
    attainable = min(peak_flops, memory_bound_ceiling)
    limiter = "compute" if memory_bound_ceiling >= peak_flops else "memory"
    return attainable, limiter

Six operations, and which wall each one hits

Here is the payoff. I counted FLOPs and bytes for six things you actually run, divided one by the other, and compared against two published H100 figures. Nothing profiled, nothing run on a GPU.

Byte counts assume the ideal case: each input read from HBM once, each output written once, everything else kept on chip. Real kernels move more, so these intensities are optimistic — which makes the bad news below worse, not better.

Ceiling  : 990 TFLOP/s dense FP16 tensor core (published peak)
Bandwidth: 3.35 TB/s HBM3 (published peak)
Ridge    : 295.5 FLOP/byte -- below this you are memory-bound

operation                         FLOP/byte     attainable   bound by
----------------------------------------------------------------------
element-wise add, n=67108864          0.167         0.6 TF/s     memory
mat-vec 4096x4096                     1.000         3.3 TF/s     memory
mat-mul 4096x4096x4096             1365.333       990.0 TF/s    compute ✓
attention, seq=512, d_head=64       256.000       857.6 TF/s     memory
attention, seq=8192, d_head=64     4096.000       990.0 TF/s    compute ✓
decode step, 7B params, batch=1       1.000         3.4 TF/s     memory

Four of six are memory-bound. The big square matrix multiply is the one operation here that genuinely uses the arithmetic units you paid for, and it is why a GPU looks fast in a benchmark. An element-wise add reaches 0.1% of the chip's arithmetic peak, because it touches every byte exactly once — which is why kernel fusion is not micro-optimisation for an add, it is the only lever that exists. Attention crosses the line as the sequence grows: seq=512 memory-bound, seq=8192 compute-bound, same kernel, same code. And the last row:

A single-token decode step sits at 1.00 FLOP/byte, roughly 296× below the ridge. Every weight is dragged out of HBM and used for two flops. The KV cache page argues decode is memory-bound from the software side, by counting what the cache costs. This is the same conclusion arriving from the hardware side, with no model in sight. Two independent routes, one answer — which is usually how you know a claim is structural rather than an artefact.

The ridge point, and why the field keeps getting more memory-bound

The two roofs meet somewhere. That crossover is the ridge point:

ridge = peak FLOP/s ÷ bandwidth   (FLOP per byte)

One number per chip per precision. Below it, memory-bound. Above it, compute-bound. It decides the fate of every kernel you will ever write, and you can compute it on a napkin from a datasheet.

Watch what it has been doing:

part                          HBM      FP16 peak ridge FLOP/byte
--------------------------------------------------------------
V100 SXM2 (2017)          0.90 TB/s        125 TF/s            139
A100 SXM  (2020)          2.04 TB/s        312 TF/s            153
H100 SXM  (2022)          3.35 TB/s        990 TF/s            296

Over those three parts:
  arithmetic throughput grew 7.9x
  memory bandwidth  grew 3.7x
  so the ridge point grew 2.1x

Arithmetic got cheap roughly twice as fast as bytes did. Which means the intensity you need in order to stay compute-bound goes up every generation. An operation that comfortably used a V100's arithmetic units can be memory-bound on an H100 without a single line of your code changing. You upgraded, and your kernel got relatively worse.

This is the thing I wish someone had told me five years ago. It explains why almost every interesting systems paper of the last five years is about moving fewer bytes rather than doing less arithmetic — FlashAttention, GQA, quantised caches, paged attention, all of it. The field did not decide to care about memory. The hardware decided for it.

Tensor cores: precision is a throughput knob

A tensor core is a unit that does a small fixed-size matrix multiply-accumulate as a single instruction, instead of one multiply-add at a time. It is much faster, and it is much faster the narrower the numbers get. Here is the same matrix multiply, four times, changing only the dtype:

Same operation: C = A @ B, 4096x4096x4096. Only the dtype changes.

precision              bytes/el        peak   ridge  FLOP/byte  bound by
------------------------------------------------------------------------
FP32 (no tensor core)         4     67 TF/s      20        683   compute
TF32 tensor core              4    495 TF/s     148        683   compute
FP16 tensor core              2    990 TF/s     296       1365   compute
FP8 tensor core               1   1979 TF/s     591       2731   compute ✓

Same maths, 137 GFLOPs either way.
Bytes moved: 201.3 MB at FP32, 50.3 MB at FP8 -- 4x less traffic.
Peak arithmetic: 67 TF/s at plain FP32 against 1979 TF/s at FP8
-- 30x more, from the published ceilings.

Look at what dropping precision does. It moves you on the roofline twice. The bytes halve, so your arithmetic intensity doubles and you slide right along the slanted roof. And the flat roof itself lifts, because the tensor cores do more per clock at narrower widths. Two independent wins from one change.

Quantisation gets sold almost entirely as a memory saving — the model fits on the card. On this hardware it is also an arithmetic saving of the same order, and that half rarely gets mentioned. The quantisation page is about what the precision drop costs you in accuracy and which numbers you are allowed to squash; this page only counts the throughput you get back.

The 30× is between two different things and I want to be honest about it. 67 TF/s is plain FP32 on the CUDA cores, no tensor core involved; 1979 TF/s is FP8 on the tensor cores. Most of that ratio is "you started using the tensor cores at all", not "you narrowed the dtype". The dtype-only comparison is the clean 2× per halving: 495 → 990 → 1979.

Occupancy: a small kernel wastes a large machine

The roofline says nothing about whether your work actually fills the chip. An H100 has 132 SMs. A thread block runs on exactly one SM. If your problem produces fewer blocks than there are SMs, the rest sit idle no matter which side of the ridge you are on.

Crudest possible model: equal-cost blocks dealt round-robin, kernel ends when the busiest SM does.

     batch     blocks    waves    utilisation       SMs idle
------------------------------------------------------------
         1          1        1           0.8%            131
         8          8        1           6.1%            124
        32         32        1          24.2%            100
        64         64        1          48.5%             68
       132        132        1         100.0%              0 ✓
       133        133        2          50.4%            131
       264        264        2         100.0%              0 ✓
       512        512        4          97.0%             16

At batch 1 you are using 0.76% of the machine and 131 of the 132 SMs have nothing to do. You did not buy a slow GPU. You bought a fast one and handed it one job.

The batch=133 row is the one that will actually bite you in production. One block past a full wave and utilisation halves, because the second wave is almost entirely empty but you wait for it anyway. Sizes that are not multiples of the SM count leave a quantised tail of idle time, and this is why tuning guides obsess over shapes that divide nicely.

This is the entire argument for batching. Batch-1 decode is memory-bound and leaves the machine idle — two separate failures with one fix. Batching raises intensity (the same weights serve more tokens) and fills the SMs. It costs you per-request latency, which is why continuous batching exists and why every serving stack has a knob labelled something like max_num_seqs. See speculative decoding for the other way out: more useful work per weight-load, without waiting for more requests.

Scaling out: NVLink is next door, InfiniBand is a lorry

One more level down the supply chain. Inside a server, GPUs talk over NVLink. Between servers, over InfiniBand. The published gap is about 18×, and that single ratio shapes every distributed training decision anyone makes.

A ring all-reduce — the gradient sync every data-parallel step — moves 2(N−1)/N × S bytes through each link, reduce-scatter then all-gather. That is bandwidth-optimal and is what the collectives library actually does. So:

Published link bandwidths (approximate, per GPU):
  NVLink (in a node)   : 900 GB/s
  InfiniBand NDR 400G  : 50 GB/s
  Ratio                : 18x

All-reduce of 7B BF16 gradients = 14.0 GB buffer.

configuration                      bytes/GPU      time
--------------------------------------------------------------
8 GPUs, 1 node, over NVLink          24.5 GB     27.2 ms ✓
  + 2 nodes, over InfiniBand         1.75 GB     62.2 ms  (56% on IB)
  + 4 nodes, over InfiniBand         2.62 GB     79.7 ms  (66% on IB)
  + 16 nodes, over InfiniBand        3.28 GB     92.8 ms  (71% on IB)

At two nodes the InfiniBand leg carries 7.1% of the bytes the NVLink leg does — and still takes the majority of the time. Inside a node, eight GPUs are effectively one machine. Across nodes they are posting parcels to each other.

That is why expert placement, all-to-all fan-out and overlapping communication with compute are genuine research problems rather than implementation details. DeepSeek-V3 is what one team did about it, on hardware where the NVLink number was cut in half again. This page is the constraint they were working against; that page is the response.

The code

Pure NumPy. No GPU required, no CUDA, nothing measured — every result is arithmetic against published datasheet figures, so it runs identically on your laptop. Every function under 12 lines.

Download the full script — about 60 lines of model plus eight runnable demos. Then:

python3 gpu_from_scratch.py
import numpy as np

# All figures are approximate, rounded, and taken from public vendor
# documentation (NVIDIA H100 Tensor Core GPU Architecture whitepaper and the
# H100 datasheet; NVLink and InfiniBand product pages). SXM form factor.
# They are PEAK numbers: the ceiling, not an achieved measurement.
H100 = {
    "sm_count":    132,              # streaming multiprocessors on the die
    "hbm_bytes":   80 * 1024**3,     # 80 GB of HBM3
    "hbm_bw":      3.35e12,          # ~3.35 TB/s memory bandwidth
    "l2_bytes":    50 * 1024**2,     # 50 MB L2, shared by all SMs
    "smem_bytes":  228 * 1024,       # up to 228 KB shared memory per SM
    "regs_bytes":  256 * 1024,       # 256 KB register file per SM
    # peak dense throughput, FLOP/s, no sparsity
    "fp32_novec":  67e12,            # plain FP32 CUDA cores, no tensor cores
    "tf32_tensor": 495e12,           # TF32 on tensor cores
    "fp16_tensor": 990e12,           # FP16/BF16 on tensor cores
    "fp8_tensor":  1979e12,          # FP8 on tensor cores
    "nvlink_bw":   900e9,            # ~900 GB/s NVLink per GPU, in a node
    "ib_bw":       50e9,             # 400 Gb/s InfiniBand NDR = ~50 GB/s
    "pcie_bw":     64e9,             # PCIe Gen5 x16, ~64 GB/s host<->device
}


# STAGE 1 -- arithmetic intensity. The whole page is this one ratio.
#     intensity = FLOPs performed / bytes moved from memory
# Units: FLOP per byte. It is a property of the OPERATION, not of the chip.
def arithmetic_intensity(flops, bytes_moved):
    return flops / bytes_moved


# STAGE 2 -- THE ROOFLINE. Williams, Waterman & Patterson, 2009.
#     attainable FLOP/s = min( peak FLOP/s,  bandwidth * intensity )
# Two straight lines. The slanted one is memory: if you can only move B bytes
# per second and you do I flops per byte, you can do at most B*I flops per
# second, no matter how good your arithmetic units are. The flat one is
# compute. The min of the two is the roof, and whichever term wins tells you
# what is limiting you.
def roofline(intensity, peak_flops, bandwidth):
    memory_bound_ceiling = bandwidth * intensity
    attainable = min(peak_flops, memory_bound_ceiling)
    limiter = "compute" if memory_bound_ceiling >= peak_flops else "memory"
    return attainable, limiter


# STAGE 3 -- the ridge point. Where the slanted roof meets the flat one.
#     ridge = peak FLOP/s / bandwidth      (units: FLOP per byte)
# Below this intensity you are memory-bound. Above it you are compute-bound.
def ridge_point(peak_flops, bandwidth):
    return peak_flops / bandwidth


# STAGE 5 -- occupancy, crudely. How much of the machine does the work fill?
# An SM is one of 132 independent processors on the die. A "warp" is a group
# of 32 threads that execute in lockstep; a kernel is chopped into thread
# blocks, and each block lands on exactly one SM. Fewer blocks than SMs and
# the leftover SMs idle. We model the crudest version: n_blocks of equal
# cost, dealt round-robin. The kernel finishes when the busiest SM finishes.
def utilisation(n_blocks, n_sms):
    if n_blocks == 0:
        return 0.0
    waves = int(np.ceil(n_blocks / n_sms))     # how many full passes needed
    capacity = waves * n_sms                    # block-slots we paid for
    return n_blocks / capacity                  # fraction actually used


# STAGE 6 -- ring all-reduce cost. The gradient sync every data-parallel step.
# It moves 2*(N-1)/N * S bytes through each link: reduce-scatter, then
# all-gather. Bandwidth-optimal, and what NCCL actually does.
def ring_allreduce_bytes(buffer_bytes, n_ranks):
    if n_ranks < 2:
        return 0.0
    return 2.0 * (n_ranks - 1) / n_ranks * buffer_bytes

Now break it

This is the part that actually makes it stick. Each row below is a belief I have personally held and been wrong about, and each one costs a number the script will print for you. Demo 7 runs all five back to back.

Break thisWatch for
Assume a faster chip fixes a memory-bound kernel: double "fp16_tensor" to 1980e12 and re-run demo 2 The decode row does not move. 3.35 TF/s before, 3.35 TF/s after — doubling the arithmetic buys exactly 0 extra TFLOP/s. Double "hbm_bw" instead and it goes to 6.70.
Serve at batch 1: call utilisation(1, 132) 0.76%. You are renting 132 processors and using one. 131 SMs idle, every step, for the whole request.
Cast to 16-bit but leave one op on the 32-bit path: compare op_matmul(4096, dtype_bytes=2) against dtype_bytes=4 with "fp32_novec" as the peak Intensity halves, 1365 → 683 FLOP/byte, and the ceiling drops 990 → 67 TF/s. A 14.8× loss of headroom from one stray operation. Both directions at once.
Move data to the GPU inside the loop: multiply the tensor bytes in demo 7 by iters instead of hoisting 33.6 MB becomes 33.6 GB over PCIe — 0.52 s of pure copying at 64 GB/s, against 0.01 s for the same bytes over HBM. 1000× more traffic than the version that hoists.
Shard across nodes when it would fit on one card: set nodes = 2 for a 14.0 GB model on an 80 GB part 1.75 GB per step over InfiniBand = 35.0 ms of interconnect that did not previously exist. Every step. Forever.
Then read the Roofline paper — it is six pages and mostly about CPUs, which is the point. Read the model and the ridge-point discussion; skip the CPU-specific ceilings and the kernel case studies, which are of their moment. After that, NVIDIA's own architecture whitepapers (the NVIDIA H100 Tensor Core GPU Architecture whitepaper is the one for the numbers on this page) are readable if you skip the marketing front matter and go to the specification tables. You'll find the paper is just the comments you've already read.

On the maths anxiety

There isn't any here, and that is worth saying plainly. This page has one division. Arithmetic intensity is FLOPs over bytes. The ridge point is peak over bandwidth. The roofline is a min() of two numbers. If you can estimate how many multiply-adds a matrix multiply does and how many bytes a tensor occupies, you have everything.

What you do need is the counting discipline, and that is a habit rather than a skill. Before you profile: how many bytes does this operation have to move, at minimum? How many flops does it do? Divide. Compare to the ridge point of the card you are on. You now know which of the two possible problems you have, and you have not started a profiler.

The interview question here is almost never "derive the roofline". It is "why is your batch-1 latency not improving on the newer GPU?" or "you doubled the FLOPs and got nothing, what happened?" The answer both times is a ratio you can compute in your head. Say the number, say which roof you are under, and say what you would change. That is the whole answer.

Where this page stops: it says nothing about warp scheduling, memory coalescing, bank conflicts, instruction-level parallelism, or the actual register-pressure arithmetic that decides real occupancy. Those matter when you are writing kernels. Almost none of them matter when you are choosing a batch size, a dtype, or a parallelism strategy — which is what most of us are actually doing.

Where to go next

This page gave you one instrument: a ratio and a ceiling. That is deliberately narrow, and it is upstream of a lot. Three of these go deeper than I have; one is the original.

Order I'd take them in: Horace He first, because it is the same model told better and adds overhead; then the NVIDIA guide to check my arithmetic against the vendor's; then the Roofline paper, which will feel obvious by then and should; then Simon Boehm's worklog only if you are actually going to write a kernel.