GPU Architecture and Programming

Volume I · The Machine

GPU Architecture & the SIMT Execution Model

Before you can size a cluster or debug a 40%-utilization training run, you have to know what a single GPU actually is: not a fast CPU, but a fundamentally different machine that trades latency for throughput. This chapter builds that machine from the transistor up.

Anyone building serious systems on GPUs is expected to move from “how many H100s does this need” all the way down to “why is this kernel only hitting 30% of peak FLOPS”, and back up to the business case, without losing the thread. That fluency starts here, with the execution model that every layer above inherits.

A modern server CPU and a datacenter GPU are both made of transistors on similar process nodes, but they spend those transistors on opposite goals. The CPU is a latency machine: it is built to finish a single stream of instructions as quickly as possible. The GPU is a throughput machine: it is built to finish an enormous number of identical operations per unit time, and it does not care how long any single one takes.

That difference cascades into the whole chip. A CPU core devotes most of its area to making one thread fast: deep out-of-order pipelines, branch predictors, speculative execution, and, above all, large caches (tens of megabytes of L3) to keep that thread fed without waiting on DRAM. A GPU inverts the ratio. It spends its area on arithmetic units, packs thousands of them onto the die, and accepts that any individual operation will stall waiting on memory. It hides that stall not with a cache but with sheer parallelism: when one group of threads stalls on a memory load, the hardware instantly switches to another group that is ready to run. This is the single most important idea in the chapter, and the dividing line between people who have read about GPUs and people who understand them.

Table 1.1 shows where the transistors go on each side of the divide:

NVIDIA calls its model SIMT, Single Instruction, Multiple Thread. It is a cousin of the SIMD (vector) model in CPUs, but with a crucial software-facing difference: you write code as if each thread is independent and scalar, and the hardware groups threads together and executes them in lockstep behind the scenes.

The unit of lockstep execution is the warp: exactly 32 threads on every NVIDIA GPU to date. All 32 threads in a warp share one program counter and issue the same instruction in the same cycle, each operating on its own data in its own registers. The warp is the true atom of scheduling, the hardware never schedules a single thread, only warps.

Threads are organised in a strict hierarchy, and every level maps to a physical resource. Internalise this, because the entire CUDA programming model in Chapter 3 is built on it:

The streaming multiprocessor (SM) is the GPU's fundamental compute building block, the rough analogue of a CPU core, except that one SM runs thousands of threads. An H100 has 132 of them; an A100 has 108. Scaling a GPU up across generations is largely a matter of adding SMs and feeding them more memory bandwidth.

The warp schedulers are where latency hiding happens. Each scheduler tracks a pool of resident warps. Every cycle it looks for a warp whose next instruction's operands are ready (not waiting on a pending memory load) and dispatches it. If a warp issues a load from HBM, hundreds of cycles of latency, the scheduler simply sets it aside and issues from a different ready warp. As long as there are enough resident warps to cover the latency, the arithmetic units never go idle. This is why “more threads” is the GPU answer to “hide more latency.”

Tensor Cores and the Transformer Engine

The plain CUDA cores do scalar FP32/INT math. The reason a GPU can train a transformer is the Tensor Core: a dedicated unit that performs a small matrix multiply-accumulate (e.g. a 4×4×4 or larger tile) in a single operation, at far higher throughput than issuing the equivalent scalar multiplies. Because virtually all the FLOPs in a neural network are matrix multiplications (the linear layers, the attention projections), Tensor Cores deliver the order-of-magnitude speedups that make large-model training feasible.

Tensor Cores have evolved each generation, and the headline is the precision they support: Volta introduced them at FP16; Ampere added TF32 and BF16; Hopper added FP8 with its Transformer Engine, a combination of FP8-capable Tensor Cores and runtime logic that automatically chooses FP8 vs 16-bit per layer and manages the scaling factors needed to keep FP8 numerically stable; Blackwell adds FP4 and FP6 and a second-generation Transformer Engine. We return to what these formats mean numerically in the precision section below.

Because a warp shares one program counter, a data-dependent branch where some of the 32 threads go one way and some go the other forces the hardware to execute both paths serially, masking off the inactive threads on each pass. This is warp divergence, and it is one of the most common reasons a GPU kernel underperforms. A branch on threadIdx.x % 2 can halve your throughput; a branch that is uniform across the warp costs nothing.

// Each thread computes one element of C = A + B.
// This is the canonical data-parallel pattern: no branches, no divergence.
__global__ void vecAdd(const float* A, const float* B, float* C, int n) {
    // Global thread index = which block + which lane within the block.
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n)                 // <-- uniform-ish guard; only the tail warp diverges
        C[i] = A[i] + B[i];    // one fused-ish op per thread, fully parallel
}

// GOOD: branch is uniform across the warp (all 32 threads take the same side),
// because warpId is constant for every lane in the warp.
__global__ void good_branch(const float* x, float* y, int n) {
    int i      = blockIdx.x * blockDim.x + threadIdx.x;
    int warpId = i / 32;                       // same value for all 32 lanes
    if (warpId % 2 == 0) y[i] = x[i] * 2.0f;   // entire warp goes one way -> no divergence
    else                 y[i] = x[i] * 3.0f;
}

// BAD: branch splits *within* the warp on the low bit of the lane id.
// The hardware executes BOTH sides serially with half the lanes masked: ~2x cost.
__global__ void bad_branch(const float* x, float* y, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (threadIdx.x % 2 == 0) y[i] = x[i] * 2.0f;   // lanes 0,2,4,...  active
    else                      y[i] = x[i] * 3.0f;   // lanes 1,3,5,...  active  (serialised!)
}

int main() {
    int n = 1 << 24;                 // 16M elements
    size_t bytes = n * sizeof(float);
    float *dA, *dB, *dC;
    cudaMalloc(&dA, bytes); cudaMalloc(&dB, bytes); cudaMalloc(&dC, bytes);
    // ... copy inputs H2D with cudaMemcpy ...

    int threads = 256;                              // block size (multiple of 32!)
    int blocks  = (n + threads - 1) / threads;      // enough blocks to cover n
    vecAdd<<<blocks, threads>>>(dA, dB, dC, n);     // triple-angle = grid,block launch
    cudaDeviceSynchronize();                        // wait for the GPU

    cudaFree(dA); cudaFree(dB); cudaFree(dC);
}

Occupancy is the ratio of warps actually resident on an SM to the hardware maximum (64 warps on recent architectures). It matters because it determines how much latency the schedulers can hide: more resident warps means more ready work to switch to when others stall. But occupancy is bounded by three finite per-SM resources, and a kernel hits the limit on whichever runs out first:

• 1
  Registers per thread. The SM has a fixed register file (65,536 32-bit registers on recent GPUs). If each thread needs 64 registers, the SM can host 65,536 / 64 ≈ 1,024 threads = 32 warps = 50% occupancy, regardless of anything else.
• 2
  Shared memory per block. If each block claims 48 KB of the 256 KB SRAM, at most ~5 blocks fit, capping resident threads.
• 3
  Block / warp slots. Hardware limits on blocks-per-SM and warps-per-SM.

The nuance that separates rote knowledge from operating experience: maximum occupancy is not always optimal. A kernel that uses more registers per thread to keep data in fast registers (doing more work per thread, a technique called register blocking) can beat a high-occupancy kernel that spills to memory. The failure mode to know by name is register spilling: when a kernel needs more registers than the SM can give at the desired occupancy, the compiler “spills” excess values to local memory (which lives in slow HBM, cached in L1), quietly destroying performance. NVIDIA's occupancy calculator and the --ptxas-options=-v compiler flag (which prints register and shared-memory usage) are the tools for tuning this.

// The CUDA runtime can tell you the best block size for max occupancy of a kernel.
int minGridSize, blockSize;
cudaOccupancyMaxPotentialBlockSize(&minGridSize, &blockSize,
                                   vecAdd, /*dynamicSMem=*/0, /*blockLimit=*/0);
// And how many blocks of a given size will be resident per SM:
int maxActiveBlocks;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&maxActiveBlocks, vecAdd, blockSize, 0);
// occupancy = (maxActiveBlocks * blockSize / 32) / maxWarpsPerSM

Half the art of large-model performance is using the lowest precision that still converges. Lower precision means fewer bytes moved (the binding constraint, see Chapter 2), more values per Tensor Core operation, and more model in the same HBM. But too few bits and training diverges. You need to be able to explain each format and when it is used.

*TF32 is stored in 32 bits but the Tensor Core only reads 19 of them (8 exponent, 10 mantissa). The two distinctions that matter most in practice:

FP16 vs BF16. They are the same size, but FP16 spends more bits on the mantissa (precision) and fewer on the exponent (range), while BF16 keeps FP32's full 8-bit exponent and sacrifices mantissa. In training, range matters more than precision: gradients span many orders of magnitude, and FP16's narrow exponent causes small gradients to underflow to zero, which is why FP16 training needs loss scaling (multiply the loss by a large constant before backprop to push gradients into the representable range, then unscale). BF16 has the same range as FP32, so it usually needs no loss scaling and has become the default. This is the cleanest “do you actually understand precision” test there is.

FP8's two flavours. Hopper's FP8 comes in two layouts because the forward and backward passes have different needs. E4M3 (4 exponent, 3 mantissa, range ±448) is used in the forward pass where precision helps; E5M2 (5 exponent, 2 mantissa, range ±57,344) is used for gradients in the backward pass where range matters. Because FP8's range is tiny, values are kept in range with per-tensor scaling factors tracked by the Transformer Engine (it keeps a short history of the maximum absolute value, the “amax,” and rescales). FP8 typically buys ~30–40% throughput and ~2× memory over BF16, but is prone to divergence when activations have outliers beyond the range, the practical mitigation is finer-grained (per-block or per-token) scaling.

This table is worth knowing cold, not every digit, but the shape: which generation introduced what, the memory capacity and bandwidth, and the relative compute. The progression tells a story: each generation adds memory bandwidth and a lower-precision format, because both directly attack the bottlenecks of ever-larger models. Figures are a 2026 snapshot, structures endure, numbers rot.

†Figures shown are the non-sparse (“dense”) rates; NVIDIA also quotes 2× “sparse” numbers that apply only with 2:4 structured sparsity. Treat all vendor TFLOPS as best-case. B200 memory is quoted as both 180 GB and 192 GB depending on the HBM3e stack configuration; prefer NVIDIA primary docs when precision matters.

Three machines worth understanding as systems, not chips

H200 is an H100 with more memory. Same GH100 compute die, same TFLOPS, the only change is 141 GB of HBM3e at 4.8 TB/s versus 80 GB at 3.35 TB/s. So H200 helps exactly when you are memory-capacity or memory-bandwidth bound (models that don't fit in 80 GB, long context, large-batch memory-bound inference) and does nothing for compute-bound work. Being able to say “H200 is a memory upgrade, not a compute upgrade” is a small thing that signals you read systems, not marketing.

GH200 Grace Hopper fuses one Grace ARM CPU (72 cores) with one Hopper GPU over NVLink-C2C at 900 GB/s, cache-coherent, roughly 7× the bandwidth of PCIe Gen5. This matters for workloads that spill to CPU memory or stream huge datasets, because the CPU's memory becomes almost a slow extension of GPU memory rather than a distant pool behind a thin PCIe straw.

GB200 NVL72 is the one to understand as an architectural shift. It is a liquid-cooled rack that wires 72 Blackwell GPUs and 36 Grace CPUs into a single NVLink domain via nine NVSwitch trays, giving ~130 TB/s of all-to-all bandwidth across the whole rack and presenting 13.5 TB of unified HBM3e. The significance: inside that rack, the usual two-tier hierarchy, fast NVLink within a server, slower InfiniBand between servers, collapses. Tensor parallelism, which turns out to be the most communication-hungry form of parallelism and normally must stay inside one 8-GPU server, can now span 72 GPUs over NVLink. That changes how you map a model onto the hardware, which is a Chapter 6 conversation.

A single H100 is often too much GPU for one small job, and sharing it well is a platform design decision. There are two mechanisms, and the distinction, hardware isolation vs cooperative sharing, is exactly what a multi-tenant platform conversation turns on.

The decision rule: MIG when you need isolation guarantees between untrusted or SLA-bound tenants; MPS when you need utilization among cooperative work you control. On every Kubernetes-based GPU platform, MIG profiles are exposed as distinct schedulable resources (e.g. nvidia.com/mig-1g.10gb), which we wire up in Chapter 8.

Why does a GPU need a quarter-million threads to stay busy? Why does FP8 stop at exactly 448? Both are short calculations. This section turns the qualitative claims above, the GPU hides latency with massive parallelism, occupancy is bounded by three resources, each format has a particular range, into the formulas behind them.

Little's Law: exactly how much parallelism hides the latency

The central claim of this chapter, “switch to another warp while this one waits on memory”, has a precise quantitative form. Little's Law, from queueing theory, relates the number of items simultaneously in a system to the rate they flow through and how long each stays:

L = \lambda\, W

Apply it to memory. To saturate HBM bandwidth, the number of bytes in flight must equal the bandwidth–latency product. This is the memory-system form of Little's Law and it dictates the minimum parallelism, the “memory-level parallelism”, a kernel needs. Derivation, for an H100:

• 1
  Take HBM bandwidth \lambda = 3.35\times10^{12} B/s and a representative global-memory latency of W \approx 400 cycles.
• 2
  Convert latency to seconds at a ~1.5 GHz memory clock: W \approx 400 / (1.5\times10^{9}) \approx 2.67\times10^{-7} s ≈ 267 ns. (The clock is approximate; the point is the order of magnitude.)
• 3
  Bytes in flight: L = \lambda W = 3.35\times10^{12}\times 2.67\times10^{-7} \approx 8.9\times10^{5} bytes, ≈890 KB must be “in the pipe” simultaneously to hide the latency.
• 4
  At 128 bytes per coalesced warp transaction, that is 8.9\times10^{5}/128 \approx 6{,}900 outstanding warp-loads, out of a hardware maximum of 132 \times 64 = 8{,}448 resident warps. You must keep the GPU ≈80% full of memory-issuing warps just to saturate HBM.

Amdahl and Gustafson: the limits of parallel speedup

Whenever you add GPUs (or SMs, or cores), two laws bound what you get. Let p be the fraction of work that is parallelizable and N the number of workers.

S(N) = \frac{1}{(1-p) + \dfrac{p}{N}} \qquad\Longrightarrow\qquad S_{\max} = \lim_{N\to\infty} S(N) = \frac{1}{1-p}

Worked example, the Amdahl ceiling for a training step. Suppose 8% of each training step is effectively serial (gradient all-reduce that doesn't overlap, optimizer step, data stalls), so p = 0.92. The maximum conceivable speedup from scaling out is:

S_{\max} = \frac{1}{1-0.92} = 12.5\times

At N = 64 GPUs you would get S(64) = 1/(0.08 + 0.92/64) \approx 10.6\times, a scaling efficiency of 10.6/64 \approx 17\%. The lesson: shrinking the serial fraction (overlap, bigger batches to amortize communication) matters far more than adding nodes once you are deep into the Amdahl regime. This is the mathematical reason communication and synchronization overheads, the “serial” parts of distributed training, are so damaging, a thread we pick up with collectives in Chapter 5 and parallelism strategies in Chapter 6.

Amdahl assumes a fixed problem. In practice we scale the problem with the machine, bigger models, bigger batches, which is the regime Gustafson's Law describes: S(N) = (1-p) + pN. Here speedup grows (nearly) linearly in N because the parallel work expands to fill the larger machine. This is why weak scaling (more GPUs for a proportionally bigger job) looks so much healthier than strong scaling (more GPUs for the same job), and why training clusters are sold on the former. The honest engineer quotes which scaling they mean.

Occupancy as a constrained optimization

The three resource limits from the occupancy section have a formal shape: resident warps per SM is the minimum over three integer-floored constraints, and occupancy is that divided by the hardware max W_{\max} = 64:

W_{\text{active}} = \min\!\left( \underbrace{\left\lfloor \frac{R_{\text{SM}}}{32\,r} \right\rfloor}_{\text{registers}},\; \underbrace{\left\lfloor \frac{M_{\text{SM}}}{m} \right\rfloor \cdot \frac{b}{32}}_{\text{shared memory}},\; \underbrace{W_{\max}}_{\text{warp slots}} \right), \qquad \text{Occ} = \frac{W_{\text{active}}}{W_{\max}}

Worked example, the register cliff. A kernel uses r = 64 registers per thread on an SM with R_{\text{SM}} = 65{,}536 registers:

W_{\text{active}} \le \left\lfloor \frac{65{,}536}{32 \times 64} \right\rfloor = 32\ \text{warps} \;\Rightarrow\; \text{Occ} = \frac{32}{64} = 50\%

Push register use to r = 32 and the ceiling doubles to 64 warps (100%). Spill to r = 128 (often what happens when the compiler runs out) and it halves to 16 warps (25%). Occupancy is a step function of register pressure, small changes cross integer thresholds and produce sudden “cliffs,” which is why -maxrregcount and the occupancy calculator exist.

Worked example, which constraint binds? A kernel uses 40 registers/thread and 16 KB shared memory per 256-thread block, on an SM with 65,536 registers, 228 KB shared memory, and a 64-warp cap. Registers: \lfloor 65{,}536/(32\times40) \rfloor = 51 warps. Shared memory: \lfloor 228/16 \rfloor = 14 blocks × 8 warps/block = 112 warps (not binding). Warp slots: 64. The minimum is \min(51, 112, 64) = 51 warps, the register limit binds, giving occupancy 51/64 \approx 80\%. Dropping to 32 registers/thread would lift the register ceiling to 64 and saturate the warp-slot cap (100%).

IEEE-754: the value of a floating-point number

Every precision format in this chapter is an instance of one formula. A floating-point number with sign bit s, a k-bit exponent field with value E, and a p-bit mantissa field with value M represents, for normal numbers:

x = (-1)^{s}\,\Bigl(1 + \tfrac{M}{2^{p}}\Bigr)\, 2^{\,E - \text{bias}}, \qquad \text{bias} = 2^{\,k-1}-1

From this single formula, every format's range and precision follow by plugging in k and p:

Notice BF16 and FP32 share an exponent (k = 8, bias 127) and therefore the same dynamic range, the formal version of “BF16 keeps FP32's range,” which is why BF16 avoids loss scaling.

Deriving FP8 E4M3's maximum of 448. E4M3 deviates slightly from strict IEEE to reclaim encodings (it has no infinities), so its top exponent code is usable for finite numbers, and only the all-ones mantissa at the top exponent is reserved for NaN:

• 1
  Exponent bits k = 4 \Rightarrow \text{bias} = 2^{3}-1 = 7. Mantissa bits p = 3.
• 2
  The largest exponent code is 1111_2 = 15; E4M3 reserves only S.1111.111 for NaN, so the largest finite mantissa code at the top exponent is 110_2 = 6. (Unlike IEEE FP16, E4M3 uses the top exponent for finite values except that one NaN slot.)
• 3
  Plug into the value formula: x_{\max} = \bigl(1 + \tfrac{6}{8}\bigr)\, 2^{\,15-7} = 1.75 \times 2^{8}, since M/2^{p} = 6/8 = 0.75 and E - \text{bias} = 8.
• 4
  x_{\max} = 1.75 \times 256 = \mathbf{448}, exactly the ±448 quoted above. By the same steps E5M2 gives 1.75 \times 2^{15} = 57{,}344.

Machine epsilon, ULP, and rounding error

The precision of a format is captured by machine epsilon \varepsilon = 2^{-p}: the gap between 1.0 and the next representable number. The unit roundoff u = \varepsilon/2 (for round-to-nearest) bounds the relative error of representing any real number:

\operatorname{fl}(x) = x\,(1+\delta), \qquad |\delta| \le u = \tfrac{1}{2}\, 2^{-p}

A ULP (unit in the last place) is the spacing between adjacent representable numbers at a given magnitude; near a value x it is about 2^{\lfloor \log_2 |x| \rfloor - p}. The key engineering fact: the error is relative, so low-precision formats are fine for values of any magnitude, until you subtract nearly-equal numbers.

Summing n numbers naïvely accumulates rounding error that grows with n. The worst-case bound for naïve sequential summation versus compensated (Kahan) summation:

\underbrace{\bigl|\hat{S}-S\bigr| \le (n-1)\,u \sum_i |x_i|}_{\text{naïve, grows with } n} \qquad\qquad \underbrace{\bigl|\hat{S}_{\text{Kahan}}-S\bigr| \le \bigl(2u + O(nu^{2})\bigr)\sum_i |x_i|}_{\text{Kahan, independent of } n}

Worked example, why a long FP16 reduction drifts. Reduce n = 10^{6} values in FP16 (u = \tfrac{1}{2} 2^{-10} \approx 4.9\times10^{-4}). The naïve relative error bound scales as (n-1)u \approx 10^{6} \times 4.9\times10^{-4} \approx 490, i.e. total loss of accuracy; the result is meaningless. The same reduction accumulating in FP32 (u \approx 6\times10^{-8}) gives (n-1)u \approx 0.06, a few percent worst case, typically far better in practice. This is the formal reason NCCL and cuBLAS accumulate in FP32, and why “compute in BF16, accumulate in FP32” is the universal mixed-precision rule, the qualitative claim of the precision section, made quantitative.

Q. Why does adding more threads improve GPU performance, when on a CPU oversubscribing threads usually hurts?

A. Because the GPU hides memory latency by switching between resident warps rather than by caching. More threads (higher occupancy) give the warp schedulers a deeper pool of ready work to issue while other warps wait on HBM, keeping the arithmetic units busy. A CPU hides latency with big caches and out-of-order execution per core, so extra software threads mostly add contention. The follow-up worth volunteering: occupancy is bounded by registers, shared memory and warp slots, and maximum occupancy isn't always optimal because register-blocked kernels can win.

Q. A model needs 100 GB of weights. Walk from the GPU spec to a recommendation.

A. 100 GB exceeds a single H100's 80 GB, so either move to H200/B200 with larger HBM, or, more likely, shard the model across GPUs, which is a parallelism decision, not just a hardware one (Chapter 6). If it's a capacity problem and the compute fits, H200's 141 GB is the cheap fix; if it's also bandwidth-bound, H200's 4.8 TB/s helps; if it's compute-bound at frontier scale, Blackwell's FP8/FP4 throughput is the lever. Before committing, establish the parameter count, sequence length and batch size, because those decide whether this is a memory, bandwidth, or compute conversation. The discipline is naming which bottleneck you are actually solving.

Q. What's the difference between FP16 and BF16, and why did the industry move to BF16?

A. Same 16 bits, different split: FP16 is 5 exponent / 10 mantissa, BF16 is 8 exponent / 7 mantissa. BF16 keeps FP32's exponent range, so gradients spanning many orders of magnitude don't underflow, which means BF16 usually trains without loss scaling while FP16 needs it. The industry favours BF16 because range matters more than mantissa precision for gradient-based training, and it's simpler and more robust.

MIG gives true spatial isolation, one tenant cannot starve another's bandwidth or see its memory, which is what untrusted, SLA-bound multi-tenancy requires. MPS raises utilization for cooperative work you control, but offers no hard memory or bandwidth isolation. Isolation vs utilization is the axis.

The next chapter takes the single most important performance idea in all of GPU computing, that almost everything is bounded by moving bytes, not by doing math, and turns it into the mental model (the roofline, Chapter 2) you will use to debug real clusters.