CUDA & GPU Programming

Your model's arithmetic looks trivial on paper: a few matrix multiplies and elementwise ops per layer.

Then you run it and the GPU utilization says 20%. Nobody is confused more than the person who just bought the chip.

So what is a GPU actually doing all day?

Here is the mental model you need, because almost every performance mystery in deep learning resolves through it:

Concretely, an H100 can execute around 3 thousand trillion simple operations per second — but its HBM memory hands over data at ~3 TB/s. Write an operation that touches 100 bytes per arithmetic op (typical for elementwise stacks, attention, decoding) and you starve the calculator 95% of the time.

To design around that starvation — to read a profiler flame graph, understand why fusing two ops doubles speed, or why LLM decode is bandwidth-bound — you need the GPU's thinking model: how it schedules threads, where its memories live, and how data moves between them. That is CUDA's worldview, and this whole topic is one guided walk through the pantry shelves.

A GPU is a throughput device. Where a CPU hides latency with big caches, branch prediction, and out-of-order execution (a handful of very smart cooks), a GPU hides latency with massive concurrency (thousands of obedient cooks).

The execution model, step by step:

• The chip contains 132–148 Streaming Multiprocessors (SMs) on H100 — each SM is a small kitchen with its own schedulers and memory.
• Work is issued in warps: 32 threads that execute the same instruction in lockstep (SIMT: Single Instruction, Multiple Threads). Like 32 cooks chopping 32 identical carrots simultaneously.
• Each SM juggles many warps at once. When one warp stalls waiting for data (a memory fetch takes hundreds of cycles), the hardware switches to another ready warp instantly — zero-cost context switch. Waiting is hidden rather than avoided.

That last trick is the whole philosophy:

Your CUDA program describes parallelism in a fixed vocabulary:

• Grid / Block / Thread: you carve the work into blocks of threads; the hardware packs blocks onto SMs and schedules warps inside them.
• Occupancy: ratio of active warps per SM versus the maximum. A kernel that hoards registers per thread leaves room for fewer warps → less latency hiding → the chip starves even at "100% utilization".
• Divergence: if threads inside one warp take different if branches, the hardware executes both paths serially, masking off inactive lanes. Branching across warps is cheap; branching within a warp is where throughput goes to die.

The most important number on a GPU is not clock speed — it is how many memory transactions a warp's load costs.

The hardware fetches memory in fixed-size chunks (32-byte sectors). A warp's 32 lanes each read 4 bytes (one float). Two worlds:

World 1 — contiguous reads (coalesced). Lane i reads address base + i·4. The 32 reads span 32 × 4 = 128 bytes = 4 sectors → 1–2 transactions. One cart picks up 32 adjacent boxes.

World 2 — strided reads. Lane i reads base + i·stride with stride 1 KB (say, column-major access of a row-stored matrix). Now the reads span 32 KB → up to 32 separate transactions. Thirty-two errands for thirty-two boxes.

Same math. Same bytes per lane. Up to 32× the memory traffic — and since memory traffic is the binding cost, up to 32× slower.

Now do one more: a "tiny" elementwise kernel over 1 MB, launched from Python.

• GPU time: 1 MB ÷ 3 TB/s ≈ 0.3 microseconds.
• Launch + orchestration overhead: ~5–10 microseconds.

You just built a Ferrari parked behind a bicycle. The compute was never the bottleneck; the paperwork was.

Both numbers point at the same life lesson, which sections 5–6 formalize: in GPU work, count bytes moved and kernels launched, not FLOPs performed.

Memory on a GPU is a stack of shelves, from "in your hand" to "across the warehouse." Slowest to fastest, with the H100 numbers you should be able to sketch from memory:

code
Copy
 1. HOST DRAM      over PCIe Gen5 x16: ~64 GB/s     ← never touch per-token
    │
 2. HBM3 "global"  80 GB, ~3.35 TB/s                ← THE number every
    │                                                  kernel minimizes
 3. L2 cache       ~50 MB, shared, automatic
    │
 4. Shared mem/L1  ~228 KB per SM, hand-managed     ← tiling arena
    │
 5. Registers      per-thread, fastest of all

The distances between rungs are brutal: host↔HBM is a 50× gap, HBM↔shared is another ~100× in effective per-access terms. Everything in GPU engineering is climbing this ladder — pulling data up once, using it many times, and never making a trip you can batch.

Where the roofline model fits: every kernel has an arithmetic intensity — FLOPs per byte touched. Plot it against the chip's FLOPs/bandwidth ratio and the kernel falls on one side of a bent line:

code
Copy
 perf ▲                    ╱ compute-bound (rare in DL)
      │                   ╱
      │  memory-bound ╱──╱   ← most attention/elementwise/decode lives HERE
      │______________╱
      └──────────────────────► FLOPs per byte moved

Modern deep learning is overwhelmingly the flat left side: memory-bound. Which is why the levers that matter are traffic reducers — coalescing, tiling, fusion — not "faster math."

Carry this picture through the rest: a banquet kitchen with ten thousand line cooks — and the pantry a long walk away (HBM).

• FLOPs = knife skill. Enormous; the cooks can chop anything instantly.
• HBM bandwidth = the corridor to the pantry. The real constraint. A thousand trips of one carrot each loses to one trip of a cart.
• Registers = the cutting board. Yours alone, fastest to reach.
• Shared memory = the per-station prep table (~228 KB per SM): the team stages what it will reuse here, by hand.
• L2 = the walk-in fridge the whole kitchen shares.
• Host DRAM = the loading dock across the parking lot (PCIe). Going there mid-service is a disaster.
• Warp = a squad of 32 cooks doing identical chops. They move as one — including to the pantry: if their boxes are adjacent, one cart; if scattered, thirty-two carts (coalescing).
• Occupancy = number of squads on shift. Half-empty kitchen means standing around during every walk.
• Divergence = one squad told "choppers chop, stirrers stir." The squad must do the chopping round, then the stirring round, half idle each time.
• Kernel fusion = prepping an entire dish at the station instead of running each half-ingredient back to the pantry between steps.
• Streams = separate order rails so the expeditor can plate course 2 while course 1 walks to the pantry.
• MIG = walling the kitchen into seven independent small kitchens; MPS = one shared kitchen, tighter choreography between two head chefs.

Every technique below is one of: fewer walks, fuller carts, more squads on shift, or a shorter corridor.

With the pantry model loaded, the industry's favorite optimizations decode instantly.

Kernel fusion is bandwidth arithmetic. y = gelu(matmul(X, W)) unfused: matmul writes a giant intermediate to HBM; gelu reads it back; downstream ops read again. Each round-trip is pantry walks. Fused into one kernel, the intermediate dies in registers — one trip eliminated per fused op. (Recall the elementwise-chain example from section 3: GPU time was 0.3 µs; the launches were the story.)

FlashAttention is the archetype. Naive attention materializes the N×N score matrix in HBM, reads it for softmax, writes it back — O(N²) traffic for O(N²) math. FlashAttention tiles the problem through shared memory and registers, never materializing the full matrix: math unchanged, memory traffic cut from O(N²) to O(N). It is the canonical proof that an "algorithm" win in 2024 was mostly a memory-hierarchy win.

LLM decode is the same law at the serving level. Each step reads weights + the entire KV cache but does one token of arithmetic — bandwidth-bound (see the KV-cache topic). Every technique there (batching, GQA, quantized caches) is this topic's vocabulary: reduce bytes per useful FLOP.

Streams, overlap, and CUDA graphs. A CUDA stream is an ordered queue of work; independent streams run concurrently — the standard recipe is cudaMemcpyAsync on one stream, compute on another, pinned host memory so the DMA can actually overlap. Events synchronize across streams. And for launch-bound steady-state loops (serving), CUDA graphs capture and replay the entire launch sequence, deleting per-step Python/orchestration cost.

Sharing one GPU across models/tenants:

• MPS (Multi-Process Service): processes share SMs cooperatively at finer granularity than default time-slicing — better joint occupancy for several small workloads.
• MIG (Multi-Instance GPU, A100/H100): partitions the GPU into up to 7 hardware-isolated instances with their own memory slices and SM quotas — predictable, noisy-neighbor-free serving, at the cost of granularity.

Knowing the machine matters. Very few teams program it in raw __global__ CUDA today:

• OpenAI Triton (@triton.jit, Python syntax over blocks/tiles): you express parallelism over whole blocks; the compiler handles warp scheduling, coalescing hints, and shared-memory staging. Production kernels — flash-style attention variants, fused optimizers, PagedAttention flavors — are commonly Triton.
• torch.compile / Inductor auto-generates Triton for elementwise chains (the compiler writing your fused kernels, per the frameworks topic).
• cuBLAS / cuDNN / CUTLASS remain the gold standard for matmul-class ops; CUTLASS exposes template-level CUDA when you truly need tiling control.
• One reproducibility trap to keep: CUDA reductions are not bitwise-deterministic — atomicAdd in floating point adds in nondeterministic order, and float addition is not associative, so results wobble run-to-run. A recurring theme in regulated ML and in "why did my eval change by 0.0001?" incidents.

The sketch of a Triton kernel below is the whole pattern: one program per tile, coalesced loads within the tile, one kernel where PyTorch launched three.

import triton
import triton.language as tl

@triton.jit
def fused_scale_bias(x_ptr, y_ptr, n, BLOCK: tl.constexpr):
    pid = tl.program_id(0)
    offs = pid * BLOCK + tl.arange(0, BLOCK)
    mask = offs < n
    x = tl.load(x_ptr + offs, mask=mask)          # coalesced within block
    tl.store(y_ptr + offs, x * 1.25 + 0.1, mask=mask)  # one kernel, not three

What the knowledge is for, in day-to-day work:

1. Diagnose before optimizing. Profile with Nsight (systems then compute). The first question is always: is this kernel memory-bound or compute-bound (roofline)? Half of "GPU is slow" tickets end at the answer.
2. Count bytes, then trips, then squad utilization. Layout fixes first (channels-last, fused QKV weights, contiguous KV blocks are all coalescing choices), then fusion, then occupancy tuning — in that order of expected yield.
3. Attack the launch-bound cases at serving time: CUDA graphs for steady-state decode loops, batch/pinned-memory overlap for pipelines, MIG or MPS for multi-tenant boxes sized to the noisy-neighbor tolerance of each customer.
4. Write Triton, not CUDA, until a profiler proves otherwise. Raw CUDA is vendor-locked and fragile across architectures (SM80 vs SM90 tuning), and occupancy-driven tuning is a per-kernel black art with diminishing returns. Reserve CUDA/CUTLASS for provably hot paths — novel attention, custom fused ops — after Nsight says Inductor or cuBLAS is not at parity.
5. Design kernels with reproducibility in mind: atomicAdd-based reductions trade bit-exactness for speed; regulated or eval-sensitive pipelines need deterministic modes or tolerance bands.

The one-line worldview: the GPU will always be fast enough at math and too slow at errands. Every architecture, every API, and every optimization in this topic is a strategy for the errands.