There is a particular kind of frustration in reading GPU kernel code. Every line is simple — an index computation, a load, a barrier — and yet the whole is opaque, because the reasons live in hardware behaviour that nobody writes down next to the code. Why 128×128 tiles? Why is the shared-memory array indexed with an XOR? Why does the fast attention kernel have two warpgroups doing different things?

This series answers those questions in order, from the bottom. It is built the way the MMA swizzle-layout post that inspired it is built: numbered sections, a diagram for every mechanism, and a refusal to hand-wave past the confusing part.

Every diagram in the series is drawn from one small palette, and colour always means the same thing. There are only two hues. Blue is the execution side of the machine — lanes, warps, warpgroups, blocks. Amber is the storage side — registers, shared memory, L2, HBM. Within each, darker means further out: a darker blue is a wider scope, a darker amber is further from the ALU and slower to reach.

That leaves pink for the one thing that is neither, the tensor core, and grey for anything idle, masked or wasted. Green and red are held back for verdicts — conflict-free versus serialised, above the ridge versus below it — and never used to name a component.

• interactive in a caption means the figure has controls. Drag them; the readout under the chart updates. Nothing is a video — every figure is live SVG.
• Code appears in tabs. CUDA is the ground truth, Triton is what most people write, and PyTorch is the reference semantics you are trying to match. Not every block has all three.
• Key numbers are sourced inline. Where NVIDIA quotes a figure "with sparsity", the dense number is shown for an apples-to-apples roofline. The higher sparse peak is real when operands satisfy NVIDIA's supported 2:4 structured-sparsity pattern; an ordinary dense LLM kernel does not receive it automatically.
• Solid fills are measurements — a bar's length or a cell's value. Pale washes with a coloured outline are containers — a box in a floorplan or a dataflow. If a shape is pale, its size means nothing.
• The dark button in the top bar switches themes; the diagrams recolour with it.

Almost every design decision in the series is downstream of one of these. If you only take one thing from Parts 1–3, take this table.

* NVIDIA does not publish B200's SM count; 148 is a third-party measurement. NVIDIA's Blackwell Tuning Guide publishes the 126 MB L2 capacity. Tensor figures here are dense; NVIDIA's 2:4 sparse figures are twice the dense peak and require compatible structured operands. Ridge point = dense peak ÷ HBM bandwidth: the arithmetic intensity above which a kernel becomes compute-bound.

You need to be able to read C and Python. You do not need prior CUDA. If you have written a Triton kernel and want to know what it compiles into, or you have read the FlashAttention paper and want to know what "warp specialization" physically means, this is aimed at you.

The series leans on primary documentation throughout — the CUDA Programming Guide, PTX ISA, CUTLASS/CuTe source, and NVIDIA's architecture tuning guides — plus the following, which are worth reading in full:

• Yifan Yang, Tensor Core MMA Swizzle Layout — the post that prompted this series.
• Simon Boehm, How to Optimize a CUDA Matmul Kernel — the canonical optimization ladder.
• Horace He, Making Deep Learning Go Brrrr From First Principles — the compute / bandwidth / overhead taxonomy.
• Vasily Volkov, Better Performance at Lower Occupancy (GTC 2010) — still the correct framing of occupancy.
• Dao et al., FlashAttention-2, FlashAttention-3, and FlashAttention-4.
• Colfax Research's CUTLASS tutorials — the best writing on Hopper-era kernel structure.

@misc{mapika2026gpukernels,
  author = {Marosi, Mark},
  title = {An Illustrated Guide to GPU Kernels for ML},
  year = {2026},
  url = {https://gpu-kernels.pages.dev/}
}

Almost every performance mistake in a GPU kernel comes from one of two wrong beliefs: that a thread is like a CPU thread, or that "more parallelism" is automatically better. Both are wrong in specific, fixable ways. This part replaces them.

We are going to build the model bottom-up, and every claim will be tied to a number you can look up. By the end you should be able to answer, for any kernel you write: how many warps are resident, what are they waiting on, and what is the hardware doing while they wait?

A single CUDA thread is slow. On an H100 the clock is around 1.98 GHz — slower than a laptop CPU — and there is no deep out-of-order machinery, no aggressive branch prediction, no big private cache to hide behind. A dependent chain of fused multiply-adds retires at roughly one every 4 cycles, same as it did a decade ago.

What the GPU has instead is thousands of independent instruction streams that cost nothing to switch between. When one stalls on memory, the hardware issues from another in the very next cycle. The CUDA Programming Guide is blunt about why this is cheap:

A CPU hides latency with caches and speculation, and it does so per instruction stream. A GPU hides latency with other instruction streams. That single sentence explains most of the design: the enormous register file, the small caches, the wide memory bus, and the fact that occupancy is a thing you have to think about at all.

The consequence for you: your job as a kernel author is almost never to make one thread faster. It is to make sure that at every cycle, some warp on every SM has work it can issue.

Keep the actual latencies in view, because every design decision later in the series is a negotiation with this table. These are measured cycle counts from published microbenchmarks, not marketing numbers:

Source: Luo et al., Benchmarking and Dissecting the NVIDIA Hopper GPU Architecture, Table IV. The Hopper column is an H800 PCIe measurement; the H100/H200 time-domain extrapolation used below is labelled separately rather than presented as a measurement of those cards.

An HBM access costs roughly 120× a register access. If a warp issues a load and then immediately needs the result, that warp is dead for ~479 cycles. Whether the SM is dead too depends entirely on how many other warps are resident — which is what the rest of this part is about.

CUDA gives you five levels of grouping. They are not equally real: some are hardware, some are pure software convention, and confusing the two is a classic source of wrong intuitions.

1.2.1 The warp is the only unit that matters for correctness intuition

Threads do not execute independently. A warp fetches one instruction and applies it across 32 lanes. That is why the warp partitioning rule is worth memorising exactly:

Warp 0 always contains threads 0–31 of the block, in order. This means a 2D block of (32, 8) gives you warps that are rows — consecutive x, fixed y. A block of (8, 32) gives you warps that straddle 4 rows of y each. That choice silently determines whether your global loads coalesce, which is the single biggest factor in Part 2.

1.2.2 Block size, and the warps you waste

Because warps are cut at 32, a block size that is not a multiple of 32 permanently wastes lanes. A block of 100 threads occupies 4 warps — 128 lanes — and 28 of them are inactive for the entire lifetime of the block. The hardware will not reclaim them.

An H100 has 132 Streaming Multiprocessors. Every one of them, on every architecture from Maxwell to Blackwell, is built the same way at the top level: four sub-partitions, each with its own warp scheduler, its own slice of the register file, and its own execution units — including memory load/store units. The four sub-partitions converge on shared SM-wide resources such as the instruction cache and the L1/shared-memory array.

That last clause matters: a warp is pinned to one sub-partition. Its registers come from that sub-partition's 64 KB slice, and its instructions issue from that sub-partition's scheduler at a maximum of one per cycle. On Hopper, four schedulers × 1 instruction gives a peak of 4 warp instructions per clock per SM.

1.3.1 Why the FP32 lane count is a 2× story

An A100 sub-partition has 16 FP32 lanes; a warp is 32 threads. So a single FP32 multiply occupies the datapath for two cycles. Hopper doubled it to 32 lanes, so the same instruction retires in one. This is directly readable from the Programming Guide's per-clock throughput table:

Two rows deserve a bookmark. INT32 throughput is half of FP32 on Hopper. Address arithmetic, loop counters and predicate math all run on that path, and in a heavily unrolled kernel they can become the bottleneck — which is exactly why Hopper added TMA to move address generation into fixed-function hardware (Part 6).

And the special-function unit does 16 results per clock per SM. On an H100 that works out to about 3.9 Tops/s of exp against 989 TFLOP/s of tensor-core matmul — FlashAttention-3 puts the gap at 256×, computed at a 1.83 GHz sustained clock rather than the 1.98 GHz boost. Hold that number; it is the entire reason FlashAttention-3 needs a ping-pong schedule (Part 7).

When threads in a warp take different paths through an if, the warp does not split into two warps. It executes both sides, one after the other, with the non-participating lanes masked off.

The cost model follows mechanically: a two-way branch where both sides are taken costs the sum of both sides, not the max. A 32-way switch where every lane picks a different case costs 32×.

1.4.1 Divergence that costs nothing

Divergence is only expensive when it happens inside a warp. This is fine:

// Every warp takes exactly one side: warp_id is uniform within a warp.
int warp_id = threadIdx.x / 32;
if (warp_id < 4) { produce(); }   // warps 0-3
else            { consume(); }   // warps 4-7

And this is not:

// Lanes within one warp disagree -> both arms run, masked.
if (x[threadIdx.x] > 0.0f) y[threadIdx.x] = fast_path();
else                      y[threadIdx.x] = slow_path();

A useful rule: ask whether the predicate is constant across each group of 32 consecutive threads. If yes, the branch is free. If no, budget for both arms.

1.4.2 Independent thread scheduling, and why __syncwarp() exists

Before Volta, all 32 threads shared one program counter and an active mask. Volta gave each thread its own PC and call stack, which fixed a class of deadlocks (a lane holding a lock while its warp-mates spin) but broke a very common optimization: warp-synchronous programming, where people omitted synchronization inside a warp because "warps are lockstep anyway."

They are no longer guaranteed to be. The convergence optimizer may group lanes back together, but you cannot rely on it. Any intra-warp data exchange must name its participants explicitly:

// WRONG on Volta and later: assumes implicit lockstep.
sdata[tid] += sdata[tid + 16];
sdata[tid] += sdata[tid +  8];

// RIGHT: explicit mask, explicit sync.
float v = sdata[tid];
for (int off = 16; off > 0; off >>= 1)
    v += __shfl_down_sync(0xffffffff, v, off);   // mask names all 32 lanes

Now we can put the pieces together. Consider the simplest possible memory-bound kernel: load a value, do a little arithmetic, store it. One warp's timeline is mostly dead air. The question is how many warps you need before the dead air is completely covered.

The governing relation is Little's Law, which Vasily Volkov put at the centre of GPU performance analysis in his 2010 GTC talk:

For memory: to saturate HBM you must keep bandwidth × latency bytes in flight at all times. Plug in real numbers:

Luo et al. measured 466 cycles on A100 (331 ns at 1.41 GHz) and 479 cycles on H800 PCIe (273 ns at 1.755 GHz). The starred H100/H200 entries reuse that H800 time as an explicit same-GH100-family estimate; they are neither direct measurements nor conversions at those cards' own clocks. Bytes in flight = bandwidth × latency, with decimal SI units: 914,550 B = 915 kB = 893 KiB. B200 latency is unpublished; the original measurement behind the directional row finds it higher than H100's.

Read the last column carefully. On an H100 at full occupancy, one 4-byte load per thread is roughly enough to saturate HBM. On a B200 it is not — you need something closer to a 16-byte vector load per thread in flight. Each generation raises the amount of memory-level parallelism a kernel must expose. That is the structural reason modern kernels use float4 loads, deep software pipelines, and TMA.

Occupancy is resident warps per SM divided by the architectural maximum (64 on Hopper and datacenter Blackwell; 48 on consumer Blackwell and Ampere GA10x). It is an upper bound on how much latency you can hide, and it is decided at launch by three independent constraints:

1. Registers. Each SM has 65,536 32-bit registers. A thread using R registers means a warp uses 32R, so the SM fits ⌊65536 / (32R)⌋ warps. At 32 registers/thread you get all 64 warps; at 64 you get 32; at the 255 maximum you get 8.
2. Shared memory. An SM has up to 228 KB on Hopper, but CUDA reserves 1 KB per resident block. Two blocks asking for 113 KB each consume 228 KB including those reservations; two asking for 114 KB would need 230 KB and do not fit. A block asking for the 227 KB user maximum fits once. A high-performance GEMM or attention kernel often chooses to sit at 1–2 blocks per SM.
3. Block granularity. At most 32 blocks per SM — the same ceiling on datacenter and consumer Blackwell, per NVIDIA's Blackwell Tuning Guide — and at most 2048 threads (1536 on consumer parts). So 64-thread blocks cap you at 32 × 2 = 64 warps only in the best case, and any register or SMEM pressure drops you below that immediately.

The binding constraint is whichever of the three is smallest. Try it:

Two practical notes. First, registers are allocated in granules, so the real allocation rounds up; use cudaOccupancyMaxActiveBlocksPerMultiprocessor rather than doing it by hand when it matters. Second, CUDA reserves 1 KB of shared memory per block, which is why the per-block maximum is 227 KB while the per-SM maximum is 228 KB. A block that needs more than 48 KB must put the excess in dynamic shared memory and explicitly opt in:

// Anything above 48 KB of shared memory must be dynamic AND opted into.
cudaFuncSetAttribute(my_kernel,
    cudaFuncAttributeMaxDynamicSharedMemorySize, 200 * 1024);
my_kernel<<<grid, block, 200 * 1024>>>(...);

// Ask the runtime what occupancy you actually got.
int blocks;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocks, my_kernel, 256, 200*1024);

1.6.1 Why the calculator disagrees with your profiler

The figure above computes theoretical occupancy: the share of warp slots that could be filled if every selected block stayed resident for the whole kernel. Nsight Compute reports a second number next to it — achieved occupancy, actual resident warps averaged over time — and on perfectly healthy kernels the two differ by 20–30%. That gap has three ordinary causes:

• The compiler knows things the calculator does not. Registers are allocated per warp in 256-register granules and rounded up, and nvcc routinely reserves more than the source implies, to keep spills off the critical path. __launch_bounds__ trades that headroom back for residency, and cudaOccupancyMaxActiveBlocksPerMultiprocessor reads the post-compilation truth rather than your guess.
• Blocks do not retire in lockstep. Occupancy touches its ceiling only while every scheduled block is alive. Blocks finish at different times and replacements trickle in, so averaged over the whole kernel the resident count sits below the peak — worst near the tail of the grid.
• Resident does not mean issuing. A block parked at __syncthreads() keeps all of its warps resident and none of them useful. Barrier-heavy kernels show high theoretical occupancy and low issue rates, and no occupancy number will tell you which is happening.

So treat the calculator as a bound on the possible, not a prediction. Whether those resident warps have anything useful to issue is a different question — the subject of the rest of this part.

Here is where a lot of tuning advice goes wrong. Little's Law says you need parallelism, and warps are only one source of it. The other is instruction-level parallelism: independent work within a single thread. A thread that has eight independent loads in flight contributes as much memory-level parallelism as eight threads with one each — while using a fraction of the warp slots.

Volkov's measurements are still the cleanest demonstration:

Vasily Volkov, Better Performance at Lower Occupancy, GTC 2010. Measured on a GTX 480; the principle has aged better than the hardware.

The honest framing: occupancy is a budget for hiding latency you cannot otherwise avoid. If your kernel has enough independent work per thread, you need less of it. If it is a pointer-chasing, low-ILP kernel, you need a lot. Measure the stall reasons, not the occupancy number.

The naive launch shape — one thread per element — is fine for teaching and rarely optimal in practice. The idiom you will see in real code is the grid-stride loop:

__global__ void saxpy(int n, float a, const float* x, float* y) {
    int stride = blockDim.x * gridDim.x;
    for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += stride)
        y[i] = a * x[i] + y[i];
}
// Launch with a grid sized to the DEVICE, not to the data:
int blocks; cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocks, saxpy, 256, 0);
saxpy<<<blocks * numSMs, 256>>>(n, a, x, y);

import triton, triton.language as tl

@triton.jit
def saxpy(x_ptr, y_ptr, a, 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)
    y = tl.load(y_ptr + offs, mask=mask)
    tl.store(y_ptr + offs, a * x + y, mask=mask)

grid = lambda m: (triton.cdiv(n, m['BLOCK']),)
saxpy[grid](x, y, a, n, BLOCK=1024)

# Ground truth to check against.
y = a * x + y

Grid-stride loops decouple the grid size from the problem size, which lets the launch be sized to the machine. That matters because of wave quantization: blocks execute in waves of roughly SMs × blocks_per_SM, and a grid of 133 blocks on 132 SMs takes two waves — the second one 99% empty.

Everything above turns into three questions you will ask of every kernel in this series:

1. How many warps are resident, and what are they waiting for? Waiting on HBM is a bandwidth problem (Parts 2–3). Waiting on SMEM is a bank-conflict or layout problem (Parts 2, 6). Waiting on the SFU is a softmax problem (Parts 4, 7). Waiting on nothing is a launch problem.
2. Is the per-thread work independent? If each thread does one dependent chain, you need occupancy. If each thread owns an 8×8 accumulator tile, you do not — and you should spend the registers instead.
3. Does the grid fill the machine, in whole waves? This is the difference between 76% and 99% utilization on a decode-shaped matmul, and no amount of inner-loop tuning fixes it.

Next: where the data actually lives, why a 4-byte load can cost you 8× the bandwidth you meant to use, and the bank-conflict pattern that every tiled kernel has to solve.

A kernel is a plan for moving data. The arithmetic is almost incidental — an H100 can issue 989 TFLOP/s of BF16 matmul but only read 3.35 TB/s from memory, so for every byte you fetch you had better find 295 floating-point operations to do with it. Everything in this part is about where a byte can live, what it costs to get it, and the two ways you can accidentally throw away 8× or 32× of your bandwidth without any warning from the compiler.

Five tiers, from the fastest thing on the chip to the slowest. The numbers are for one H100 SXM5.

* L2 bandwidth is an estimate, not an NVIDIA specification: Luo et al., Table V measure 3,942.4 B/clock aggregate on H800; applying an H100 1.98 GHz clock gives about 7.8 TB/s. HBM bandwidth is the published NVIDIA H100 figure.

Threads do not read bytes from DRAM. The memory system reads sectors, and a sector is 32 bytes, naturally aligned. Four sectors make a 128-byte cache line. The Programming Guide states the rule directly:

So the question for every load instruction is: how many distinct 32-byte sectors do the 32 lanes of this warp touch? The best case is four — 32 lanes × 4 bytes = 128 contiguous bytes. The worst case is thirty-two, and then you are using one eighth of your bandwidth to move data you will throw away.

Three cases are worth internalising. Write the byte address of lane l as base + 4·l and count the sectors it crosses:

• Stride 1, aligned — lane l reads byte 4·l, so the warp spans bytes 0–127: sectors 0, 1, 2, 3. Four sectors, 100% efficiency. This is what y[i] = a*x[i] + y[i] gives you when i = blockIdx.x*blockDim.x + threadIdx.x.
• Stride 1, misaligned by 4 bytes — lane l reads byte 4·l + 4; the span 4–131 now touches sectors 0 through 4. Five sectors, 80% efficiency. Costs you one extra sector per warp; often partially recovered from L2 because the neighbouring warp already pulled that line.
• Stride 8 floats — lane l reads byte 32·l, and since one sector is exactly 32 bytes every lane opens a fresh one: 32 sectors, 12.5% efficiency. Any larger stride behaves the same. This is the column-major access pattern, and it is the single most common reason a "simple" kernel runs 8× slower than it should.

2.2.1 Why the block shape decides this

Recall from Part 1 that warps are cut along the linear thread index, so threadIdx.x varies fastest within a warp. That makes the following two kernels — identical in every other respect — differ by 8×:

// dim3 block(32, 8);  A is row-major, N columns.

// COALESCED: within a warp, tid.y is fixed and tid.x runs 0..31,
// so the 32 lanes read 32 consecutive floats = 4 sectors.
int row = blockIdx.y * 8 + threadIdx.y;
int colM = blockIdx.x * 32 + threadIdx.x;
float v = A[row * N + colM];

// UNCOALESCED: swap the roles. Now the 32 lanes of a warp read
// addresses N*4 bytes apart -> 32 separate sectors.
float w = A[colM * N + row];

# Triton makes the same distinction, just with block pointers.
offs_m = pid_m * BLOCK_M + tl.arange(0, BLOCK_M)
offs_n = pid_n * BLOCK_N + tl.arange(0, BLOCK_N)

# Coalesced: the LAST axis of the offset tensor is contiguous in memory.
a = tl.load(A + offs_m[:, None] * stride_am + offs_n[None, :] * stride_an)

# If stride_an != 1 you have a strided load. Triton will still
# emit it, and it will still cost you sectors.

A lane can load 1, 2, 4, 8 or 16 bytes in a single instruction, provided the address is naturally aligned to that size. Using the wide forms does two things at once: it cuts the instruction count by 4×, and — more importantly — it quadruples the bytes each thread has in flight.

That second effect is Little's Law again. From Part 1: under the stated latency estimate, an H100 needs ~915 kB (893 KiB) of memory traffic in flight to saturate HBM, which is 3.4 bytes per thread at full occupancy. A B200 needs more than 15. A kernel where each thread has one 4-byte load outstanding is structurally incapable of saturating a B200's memory system, no matter how good its access pattern is.

// 1 instruction, 16 bytes per lane, 128 bytes per QUARTER-warp.
// Requires x to be 16-byte aligned and n divisible by 4.
const float4* x4 = reinterpret_cast<const float4*>(x);
float4 v = x4[i];                 // compiles to LDG.E.128
v.x *= a; v.y *= a; v.z *= a; v.w *= a;

// cudaMalloc always returns >=256-byte-aligned pointers, so the base
// is fine; it is your OFFSETS that break alignment.

One warp issuing LDG.E.128 requests 512 bytes — sixteen sectors, four cache lines — in one instruction. Same efficiency and payload as four separate warp-wide 4-byte coalesced load instructions, one quarter of the instructions, four times the memory-level parallelism per thread. This is why many high-performance kernels use float4 or int4 as their unit of global traffic when alignment and layout permit it.

2.4.1 L1 and shared memory are the same silicon

L1 and shared memory have been one physical array since Volta, split at kernel launch between hardware-managed cache and programmer-managed scratchpad. The array has grown every generation — 96 KB on Volta, 192 KB on A100, 256 KB per SM on Hopper and Blackwell (128 KB on consumer parts) — and so has the scratchpad's maximum share: 164 KB of shared memory per SM on A100, 228 KB on H100 and B200. Selecting the 228 KB carveout leaves 28 KB of L1; because CUDA reserves 1 KB per block, one block can request at most 227 KB of it. That is a real trade, and for a well-tiled kernel it is usually the right one — you are replacing a cache that guesses with a scratchpad that knows.

// Static __shared__ arrays are capped at 48 KB. Above that: dynamic + opt-in.
cudaFuncSetAttribute(kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, 227*1024);

// The carveout is a HINT; percentages round UP to a supported capacity.
cudaFuncSetAttribute(kernel, cudaFuncAttributePreferredSharedMemoryCarveout, 100);

2.4.2 L2 is big, shared, and not uniform

40 MB on an A100, 50 MB on an H100, ~126 MB on a B200 — large enough that for many LLM decode steps the entire activation working set is L2-resident, and only the weights come from HBM. But L2 is physically split into two partitions, each closer to half the GPCs. Crossing that boundary costs measurably more latency and bandwidth, and on a B200 the split is also a die boundary: measured 21 TB/s within a partition against 16.8 TB/s across it.

NVIDIA publishes the 126 MB GB200 L2 capacity in its Blackwell Tuning Guide. The 21/16.8 TB/s partition-bandwidth figures are separate original Chips and Cheese measurements, not NVIDIA specifications.

This is why tile scheduling order matters. A GEMM that walks output tiles in row-major order streams a fresh row of B through L2 for every tile; one that walks them in small 2D groups reuses what is already there. Part 5 measures the difference.

You also get an explicit lever: a persistence window that marks a range of global memory as preferentially retained in L2.

// Set aside part of L2 for persisting accesses (max 75% of L2 on A100/H100).
cudaDeviceSetLimit(cudaLimitPersistingL2CacheSize, size);

cudaStreamAttrValue attr;
attr.accessPolicyWindow.base_ptr  = kv_cache;      // hot, reused every step
attr.accessPolicyWindow.num_bytes = bytes;
attr.accessPolicyWindow.hitRatio  = 0.6f;         // <1.0 avoids LRU thrash
attr.accessPolicyWindow.hitProp   = cudaAccessPropertyPersisting;
attr.accessPolicyWindow.missProp  = cudaAccessPropertyStreaming;
cudaStreamSetAttribute(stream, cudaStreamAttributeAccessPolicyWindow, &attr);

At the instruction level, PTX exposes cache operators for the same idea — ld.global.cg to skip L1 entirely, .cs ("streaming", evict-first) for data you will touch once, .L1::evict_last for data you want kept. In CUDA C++ these surface as __ldcg, __ldcs, __stwt, and cuda::annotated_ptr.

Shared memory is a scratchpad: no tags, no eviction, no misses. It is organised as 32 banks, each 4 bytes wide, each serving one 32-bit word per clock. Successive 32-bit words live in successive banks, so:

A warp-wide access is conflict-free if the 32 lanes hit 32 different banks — or the same word, which is broadcast for free. Anything else serialises:

So an N-way conflict costs exactly N×. The pathological case is 32-way, and you hit it with the most natural code imaginable.

2.5.1 Why the column access is 32-way

With __shared__ float tile[32][32], element tile[r][c] lives at word index 32r + c, so its bank is (32r + c) mod 32 = c. The row index has no effect on the bank. Now run the two natural warp accesses lane by lane, say for row/column 7:

A warp reading a column asks one bank for 32 different words, and gets them one per cycle.

2.5.2 Padding: the cheap fix that breaks later

Declare tile[32][33] instead. Now tile[r][c] is at word 33r + c, bank (33r + c) mod 32 = (r + c) mod 32. Reading column c gives banks c, c+1, c+2, … — all 32 distinct. Conflict gone, for the price of 3.1% of your shared memory.

This one-float pad is right for the hand-written FP32 access above. It does not compose with the standard tensor-core layouts: a 33-float row is 132 bytes, so every row after the first loses the 16-byte alignment expected by ldmatrix, wgmma descriptors and TMA layouts. Padding is not categorically forbidden if the padded stride preserves every required alignment, but the naive +1 fix does not. Which brings us to the usual alternative.

A swizzle permutes addresses within a power-of-two block instead of expanding the block. The trick is to XOR the 16-byte-chunk index with the row index: for any fixed row, c ⊕ r is a permutation of that row's chunks — apply it twice and you are back where you started — so every element still has exactly one home and the array stays dense. And because the permutation moves whole 16-byte chunks and never touches the bits inside a chunk, 16-byte alignment survives.

CUTLASS expresses the whole family with three integers:

M is how many low bits to leave alone, B is how many bits participate, S is the shift. Every swizzle mode that TMA and wgmma accept uses M = 4 (the low 16 bytes are never permuted) and S = 3 (the pattern repeats every 8 rows):

The 128 B mode is the one you want by default: its atom is 8 rows × 128 bytes, which matches both a full cache line on the global side and a full conflict-free warp access on the shared side. In C it is one line:

The Swizzle<B, M, S> parameterisation is CUTLASS’s (cute::Swizzle); the same three non-trivial patterns are the CU_TENSOR_MAP_SWIZZLE_32B/64B/128B values accepted by cuTensorMapEncodeTiled.

The figure below is the important one in this part. It shows a 128-byte-swizzled 8×8 grid of 16-byte chunks, and what happens when a tensor-core load asks for the same logical chunk from all 8 rows — which is exactly what ldmatrix.m8n8 does.

One constraint to know before you meet it: with 128 B swizzle the innermost TMA box dimension must be ≤ 128 bytes, which for BF16 means at most 64 elements. A BLOCK_K = 128 BF16 K-major tile therefore needs a split or reoriented box, a different contiguous dimension, or no swizzle. A smaller non-zero swizzle does not help: the Driver API limits tighten to 64 and 32 bytes for the 64 B and 32 B modes. Violating this constraint is a common cause of CUDA_ERROR_INVALID_VALUE out of cuTensorMapEncodeTiled.

Historically, staging a tile from global to shared memory meant a round trip through registers: LDG into a register, STS out to shared. That burns registers, burns issue slots, and forces the thread to wait.

Ampere added cp.async, which copies global → shared without touching the register file, asynchronously, with group-based completion tracking. Hopper added TMA, which goes further: a single thread describes a multidimensional tile copy, and dedicated hardware generates all the addresses, handles the out-of-bounds edges, applies the swizzle, and signals a barrier when the bytes have landed.

// Each thread copies 16 B global -> shared, bypassing registers.
unsigned smem_addr = __cvta_generic_to_shared(&smem[stage][tid*4]);
asm volatile("cp.async.cg.shared.global [%0], [%1], 16;\n"
             :: "r"(smem_addr), "l"(&gmem[idx]));

asm volatile("cp.async.commit_group;\n");       // close this stage
asm volatile("cp.async.wait_group 2;\n");       // keep 2 stages in flight
__syncthreads();

// One thread issues the whole tile copy; hardware does the rest.
if (threadIdx.x == 0) {
    mbarrier_init(&bar[stage], 1);
    mbarrier_expect_tx(&bar[stage], TILE_BYTES);      // byte-counting barrier
    cp_async_bulk_tensor_2d(smem[stage], tensor_map, k_off, m_off, &bar[stage]);
}
// Everyone waits on the phase bit; no per-thread address math at all.
mbarrier_wait_parity(&bar[stage], phase);

Both exist to enable the same structure: a multi-stage software pipeline, where the copy for tile i+2 is in flight while the math for tile i is running. This is the single most important structural idea in modern GEMM and attention kernels, and it is why shared memory is allocated as a circular buffer of 3–6 stages rather than one tile.

Hopper added a level between block and grid. A thread block cluster is up to 8 blocks (16 non-portably) guaranteed to be co-scheduled on one GPC, and — the actual point — every block in a cluster can load, store and do atomics directly on the shared memory of its peers, over a dedicated SM-to-SM network. This is distributed shared memory (DSMEM).

Measured on an H800: SM-to-SM latency 180 cycles, a 32% reduction versus going through L2, and 3.27 TB/s of throughput at cluster size 2. A size-8 cluster on an H100 gives a kernel 1.8 MB of single-hop, software-managed memory.

Source: Luo et al., Benchmarking and Dissecting the NVIDIA Hopper GPU Architecture, §III-D3; measured on an H800 PCIe.

The practical payoff for GEMM is TMA multicast: one block issues a TMA load and the hardware delivers the same tile into the shared memory of every block in the cluster. Two SMs cooperating on one B-tile load halves the L2 traffic for that operand — worth 704 → 734 TFLOP/s in one benchmark in this H100 GEMM worklog. Treat that delta as configuration-specific, not a universal multicast gain.

__global__ void __cluster_dims__(2, 1, 1) kernel(
    const float* input, float* output) {
    namespace cg = cooperative_groups;
    cg::cluster_group cluster = cg::this_cluster();

    extern __shared__ float smem[];
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    smem[threadIdx.x] = input[i];

    cluster.sync();                     // every peer exists and its SMEM is initialized
    float* peer = cluster.map_shared_rank(smem, cluster.block_rank() ^ 1);
    float v = peer[threadIdx.x];        // ~180 cycles, no L2 round trip
    output[i] = v;
    cluster.sync();                     // no block exits while a peer still reads its SMEM
}

Everything above compresses into six questions to ask of any kernel:

1. Does a warp's global load span 128 contiguous bytes? If not, count the sectors. Fix the block shape or transpose the data before you fix anything else.
2. Are you loading 16 bytes per thread per instruction? If you are loading 4, you have a quarter of the memory-level parallelism you could have, and on Blackwell that is disqualifying.
3. Does a warp's shared-memory access hit 32 distinct banks? If it reads a column of a power-of-two-strided array, it does not.
4. If tensor cores are involved, does the shared layout satisfy their alignment? The naive one-element padding fix breaks the standard ldmatrix, wgmma and TMA layouts; a compatible swizzle usually avoids both the conflict and the extra stride.
5. Is there anything in flight while you compute? A single-stage loop pays full memory latency on every iteration.
6. Are you re-reading from HBM what L2 or SMEM already has? That is a scheduling-order question, not a bandwidth question — Part 5.

Next: how to tell, before writing a line of code, whether a kernel can be fast — and which of the two ceilings it will hit.

Before you optimize a kernel you should be able to state, in one sentence, what is stopping it. There are only three candidate answers, and a single quantity decides between the first two. Getting this right early saves you from the classic failure mode: spending a week on tensor-core utilization for a kernel that was never going to be limited by arithmetic.

Arithmetic intensity is the number of floating-point operations a kernel performs per byte it moves between the GPU and its memory:

The machine has a matching number, the ridge point: peak arithmetic throughput divided by peak memory bandwidth. It is the intensity at which the two ceilings cross.

Below the ridge you are memory-bound: the arithmetic units are idle waiting for data, so you win by moving fewer bytes or by using more of the available bandwidth. Above it you are compute-bound: memory keeps up, so you win by doing the arithmetic more efficiently (better instructions, better tiles, lower precision).

Two properties of this plot are worth stating explicitly, because they are the whole reason the model is useful:

• Horizontal moves are optimizations that change what you touch — tiling, fusion, caching, quantizing weights. They move a kernel right, toward the ridge.
• Vertical moves are optimizations that change how well you use the ceiling you are already under — coalescing, bank-conflict removal, better instruction mix. They move a kernel up toward the roof it is already beneath, and they cannot move it above it.

A kernel at 4 FLOP/byte running at 90% of the bandwidth ceiling is done. There is nothing more to win without restructuring the algorithm.

The plot has a pedigree: it is the Roofline model of Williams, Waterman and Patterson (CACM 2009), drawn originally for multicore CPUs and the GPUs of that era. It has since been promoted into Nsight Compute, which will draw this exact chart for any kernel you point it at.

Horace He's taxonomy, from his post Making Deep Learning Go Brrrr From First Principles, is the practical version of the model: every slow kernel is slow for one of three reasons, and each has a different fix.

The third regime — overhead-bound — is invisible to the roofline model because it is not about the GPU at all. A kernel launch costs a few microseconds of CPU and driver work. A transformer layer with 20 small unfused operations, each taking 3 µs of GPU time and 5 µs of launch overhead, spends most of its wall clock in Python and the driver. This is the normal state of an unoptimized PyTorch model at small batch size, and no kernel tuning fixes it — CUDA graphs, torch.compile, and fusion do.

Here is the whole table for a transformer, worked out. Assume BF16 (2 bytes) unless stated, and count only HBM traffic — the intensity of a well-tiled kernel is set by what crosses the DRAM boundary, not by what crosses the SMEM boundary.

Read the last three rows together, because they are the central fact of LLM serving. Long-sequence or well-batched prefill is typically a compute problem; low-batch decode is typically a bandwidth problem, and they are the same mathematics with different shapes. A single-sequence decode step reads the entire weight matrix to do two FLOPs per element. Its arithmetic intensity is 1. The ridge point is 295. You are using 0.34% of the machine's arithmetic, and that is the optimal outcome for that shape.

For this weight-stationary shape, the main reuse lever is batching — sharing that weight read across more queries. Quantizing the weights is the other direct way to reduce the compulsory bytes.

This is why inference servers use continuous batching, and why FP8 weights improve both sides of the model: half the weight bytes help the memory-bound side, while twice the compute ceiling helps after the knee. Both effects are about 2× in this idealized model, not a compounded 4×; because FP8 doubles both arithmetic intensity and the ridge, the crossover batch stays roughly unchanged. An H200 — the same peak compute as H100 with 43% more bandwidth — is dramatically better at low-batch decode and can also help bandwidth- or capacity-limited training phases, although it does not raise the ceiling for compute-bound training kernels.

A naive matmul computes each output element independently: it reads a full row of A and a full column of B for every one of the MN outputs, so its A/B operand reads move 2MNK elements. Its intensity is 1 FLOP per operand element read — hopeless. The output write is omitted for this first reuse derivation and added back when counting total DRAM traffic.

Now give each thread block a T×T output tile. The block reads TK elements of A and KT of B, and there are (M/T)(N/T) blocks:

In intensity terms, for a BM × BN tile of depth BK and element size s:

Note what is missing from that expression: BK cancels out of this SM-inbound A/B operand model. It still affects shared-memory footprint, pipeline depth, edge handling, and how output traffic is amortized in a complete kernel. The harmonic-mean form is maximised for a fixed accumulator area when BM = BN, so square-ish tiles get more reuse from that fixed area. Real CUTLASS shapes also answer to tensor-instruction geometry, registers, SMEM, occupancy and problem shape; intensity is one constraint, not the entire explanation.

For memory-bound work, the optimization is always the same: touch DRAM fewer times. A chain of k elementwise operations, each launched as its own kernel, reads and writes the whole tensor k times. Fused into one kernel, it reads and writes once.

Concretely, for a Llama-style block at hidden size 8192, batch 8192 tokens, BF16 — one activation tensor is 8192×8192×2 = 134 MB:

The same arithmetic as an audit you can run — four materialized PyTorch operations on one side, one Triton kernel on the other:

# Four high-level eager operations. Profile your PyTorch/CUDA version:
# a native reduction may use more than one internal kernel.
h = torch.nn.functional.rms_norm(x, (x.shape[-1],), weight=None, eps=eps)
s = torch.nn.functional.silu(h)
t = s * g
y = t + x

# DRAM audit, bf16, one activation tensor = 8192×8192×2 B = 134 MB:
#   rmsnorm   read x, write h          268 MB
#   silu      read h, write s          268 MB
#   mul       read s and g, write t    402 MB
#   add       read t and x, write y    402 MB
#   10 passes = 1342 MB ⇒ floor 400 µs at 3.35 TB/s

@triton.jit
def rms_silu_gate_add(X, G, Y, stride, H, eps, BLOCK: tl.constexpr):
    row  = tl.program_id(0)
    cols = tl.arange(0, BLOCK)
    mask = cols < H
    x = tl.load(X + row*stride + cols, mask=mask, other=0.).to(tl.float32)
    g = tl.load(G + row*stride + cols, mask=mask, other=0.).to(tl.float32)
    inv = tl.math.rsqrt(tl.sum(x*x, axis=0) / H + eps)
    h  = x * inv
    tl.store(Y + row*stride + cols, (h * tl.sigmoid(h) * g + x).to(tl.bfloat16), mask=mask)

# 3 passes: read x and g, write y = 402 MB ⇒ floor 120 µs.

The same logic drives GEMM epilogue fusion: the accumulator is already in registers at the end of the mainloop, so applying a bias, an activation, a residual add or a quantization step there costs nothing but a few instructions. Writing it out and reading it back in another kernel costs two full passes over the output. Part 5 shows where the epilogue sits in a real kernel.

And it is the reason FlashAttention exists at all: standard attention writes the N×N score matrix to HBM and reads it back, then softmax writes and reads an equally large probability matrix — four one-way N×N transfers, or two round trips, about 17 GB for seq 8192 and 32 heads, none of it useful output. Part 7.

The ridge point is the single most useful summary of a GPU generation, and it has a surprising history: it roughly doubled from A100 to H100, meaning Hopper made it substantially harder to be compute-bound. Blackwell then held it steady for BF16 while pushing it far out for FP4.

Two consequences worth carrying:

1. Lower precision raises both ceilings. FP8 weights halve the bytes, so bandwidth-bound decode gets about 2×; FP8 tensor cores also double peak compute, so the compute-bound side gets about 2×. Arithmetic intensity and ridge both double, however, so in this idealized weight-stationary model the batch-size crossover stays near the same point rather than doubling.
2. The gap between generations is mostly not in the ridge. B200 gives you 2.3× the BF16 FLOPs and 2.3× the bandwidth of an H100. If your kernel was memory-bound on H100 it is still memory-bound on B200 — just 2.3× faster. Restructuring the algorithm is the only thing that changes which ceiling you are under.

The procedure, before writing any code:

1. Count the FLOPs. For a matmul, 2MNK. For attention, 4N²d per head forward, ×2.5 for backward.
2. Count the compulsory bytes — the inputs and outputs that must cross the DRAM boundary at least once. Ignore intermediates you intend to keep on chip.
3. Divide, and compare to the ridge. Below it, your target is bytes / bandwidth. Above it, your target is FLOPs / peak.
4. Now measure. If you are at 85% of the relevant ceiling, stop. If you are at 20%, the gap is a Parts 1–2 problem: coalescing, bank conflicts, occupancy, or an empty pipeline.

That last step is the discipline the model buys you. "This kernel takes 400 µs" is not actionable. "This fused kernel must move about 402 MB, so it cannot beat roughly 120 µs at peak H100 bandwidth, and it takes 400" is a bug report.

Next: the kernels that live at the bottom of the roofline — elementwise, reductions, softmax, RMSNorm — and how to actually reach the bandwidth ceiling rather than merely aiming at it.

Everything in this part has an arithmetic intensity between 0.2 and 3. For a large enough tensor, Part 3 therefore gives the first target: compulsory bytes divided by achievable bandwidth. Reaching it still requires coalesced traffic, enough work in flight, and little unnecessary movement. A naive elementwise kernel gets maybe 60% of peak; a naive reduction gets 5%.

Start with a copy. Its job is to read n bytes and write n bytes, so on an H100 the floor for a 1 GB tensor is 2×10⁹ / 3.35×10¹² = 597 µs. Whether you get near it depends entirely on how much traffic you keep in flight, which is Little's Law from Part 1:

Two knobs get you there: how many warps are resident, and how many independent loads each thread has outstanding. They multiply.

In code, the two knobs are the vector width and the unroll factor:

// A grid-stride copy with 16-byte accesses and 4 loads in flight per thread.
__global__ void scale4(int n4, float a, const float4* __restrict__ x, float4* __restrict__ y)
{
    int stride = blockDim.x * gridDim.x;
    for (int i = blockIdx.x*blockDim.x + threadIdx.x; i < n4; i += stride*4) {
        float4 a0, a1, a2, a3;                // four INDEPENDENT loads —
        a0 = x[i];                            // the compiler issues all four
        a1 = (i+stride   < n4) ? x[i+stride  ] : a0;   // before waiting on any
        a2 = (i+stride*2 < n4) ? x[i+stride*2] : a0;
        a3 = (i+stride*3 < n4) ? x[i+stride*3] : a0;
        a0.x*=a; a0.y*=a; a0.z*=a; a0.w*=a;   // scale every loaded vector
        a1.x*=a; a1.y*=a; a1.z*=a; a1.w*=a;
        a2.x*=a; a2.y*=a; a2.z*=a; a2.w*=a;
        a3.x*=a; a3.y*=a; a3.z*=a; a3.w*=a;
        y[i] = a0;
        if (i+stride   < n4) y[i+stride  ] = a1;
        if (i+stride*2 < n4) y[i+stride*2] = a2;
        if (i+stride*3 < n4) y[i+stride*3] = a3;
    }
}

# Triton gets the same effect from a large BLOCK_SIZE: the block is
# lowered to vectorized loads and the loads are issued together.
@triton.jit
def scale(x_ptr, y_ptr, a, 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)
    tl.store(y_ptr + offs, a * x, mask=mask)

# BLOCK=4096 with num_warps=8 gives each thread 16 elements in flight.
# This is what triton.autotune searches over — see Part 9.

y = a * x     # dispatches to a hand-tuned elementwise kernel
              # torch.compile will fuse this with its neighbours

A reduction turns n values into one. It is the shape underneath sum, max, norm, softmax, and every loss function, and it is the standard worked example of GPU optimization because the naive version is so bad.

The difficulty is that a reduction is inherently sequential in its dependency structure, so you have to build a tree. How you index that tree determines whether your warps diverge and whether your shared-memory accesses conflict — both of which you now know how to reason about.

The modern answer is the third one. Since Kepler, a warp can exchange registers directly:

// Sum a contiguous set of active lanes starting at lane 0.
// A full warp takes 5 instructions; the result lands in its first lane.
__inline__ __device__ float warpReduceSum(float v, unsigned mask) {
    int lane = threadIdx.x & 31;
    int count = __popc(mask);
    int off = 1;
    while (off < count) off <<= 1;
    for (off >>= 1; off > 0; off >>= 1) {
        float other = __shfl_down_sync(mask, v, off);
        if (lane + off < count) v += other;
    }
    return v;
}

// Two-level: reduce within each warp, stage up to 32 partials in SMEM,
// then have warp 0 reduce those. Works for any block size up to 1024.
__inline__ __device__ float blockReduceSum(float v) {
    __shared__ float s[32];
    int lane = threadIdx.x & 31, wid = threadIdx.x >> 5;
    int lanesThisWarp = min(32, blockDim.x - wid*32);
    unsigned warpMask = lanesThisWarp == 32
        ? 0xffffffffu : ((1u << lanesThisWarp) - 1u);
    v = warpReduceSum(v, warpMask);
    if (lane == 0) s[wid] = v;
    __syncthreads();
    int numWarps = (blockDim.x + 31) >> 5;
    float blockSum = 0.0f;
    if (wid == 0 && lane < numWarps) {
        blockSum = s[lane];
        unsigned finalMask = numWarps == 32
            ? 0xffffffffu : ((1u << numWarps) - 1u);
        blockSum = warpReduceSum(blockSum, finalMask);
    }
    return blockSum;                  // valid in thread 0
}

# In Triton the reduction is a single call; the compiler emits the
# shuffle tree, the shared-memory staging, and the barriers for you.
total   = tl.sum(x, axis=0)
maximum = tl.max(x, axis=0)

# Associative custom reductions are available too:
@triton.jit
def _combine(m0, d0, m1, d1):
    m = tl.maximum(m0, m1)
    return m, d0 * tl.exp(m0 - m) + d1 * tl.exp(m1 - m)

The step count is worth stating exactly, because it is the whole complexity story. Within a full warp: log₂32 = 5 shuffle steps. Within a block: those 5, plus one shared-memory round trip, one barrier, and up to 5 more shuffles for the warp partials. Across the grid: each block writes its single partial to global memory, and a second tiny kernel over the B partials finishes the job — or each block does one FP32 atomicAdd, which saves the launch but can make the result run-to-run nondeterministic because floating-point addition is not associative and block arrival order can change. The two-kernel path adds a logical 8B bytes for writing and rereading the partials; atomics instead perform contended read-modify-writes whose physical traffic depends on the cache and memory system. Tree depth is O(log N) at every level, but only the first level touches the full input; that is the difference between a tuned reduction sitting on bandwidth and a naive one sitting at 5% (Harris 2007).

One variant to know: built as a butterfly with __shfl_xor_sync, the same 5 steps leave the sum in every lane instead of only lane 0 — which is what you want when all threads consume the result, as in the norms of §4.5 (those broadcast through shared memory instead). Same depth either way; the choice is only about where the answer has to end up.

Softmax needs a maximum for numerical stability, a sum of exponentials for normalization, and then a division. Written naively that is three passes over the data:

Milakov and Gimelshein's observation is that the first two passes can be fused, if you are willing to retroactively correct the running sum every time the maximum increases:

That bolded factor is the whole trick. When the running max jumps from m_j−1 to a larger m_j, every term already accumulated in d was scaled by the wrong exponent; multiplying by exp(m_j−1 − m_j) fixes all of them at once. And because the exponent is always ≤ 0, the correction never overflows.

The state is a pair (m, d), and the combining operator is associative and commutative:

which means it drops straight into the shuffle tree from §4.2. Every thread keeps its own partial (m, d), and one warp reduction combines them. That is the entire parallel online softmax.

The payoff is memory accesses per element: 4 for the safe 3-pass version, 3 for the online version, and the paper measures a 1.15–1.3× end-to-end speedup for softmax alone. That is a modest win here — and a decisive one in attention, where the "array" is a score matrix you cannot afford to store at all.

From Part 1's instruction table: the special-function unit produces 16 results per clock per SM. On an H100 that works out to about 3.9 trillion exp results/s, against 989.5 TFLOP/s of BF16 matmul — unlike units, but a useful throughput ratio for the operation mix that FlashAttention-3 rounds to 256×.

Practical consequences:

• Use exp2, not exp. The hardware instruction is MUFU.EX2. Computing e^x costs an extra multiply by log₂e — so fold that constant into the scale you were already applying. Attention kernels multiply QK^T by 1/√d · log₂e and then call exp2.
• Prefer __expf/-use_fast_math semantics where accuracy allows; the accurate expf is a multi-instruction sequence.
• Overlap it with matmul. This is exactly what FlashAttention-3's ping-pong schedule does: while one warpgroup runs its GEMM on the tensor cores, the other runs its softmax on the SFU. Two different units, so they genuinely run concurrently. Part 7.

RMSNorm is the modern default (Llama, Mistral, Gemma, DeepSeek) because it drops the mean-centering entirely:

One reduction instead of two, no mean subtraction, same normalization effect in practice. It is entirely memory-bound: for hidden size H it does O(H) FLOPs on 2H·s bytes, an intensity around 0.75.

// One block per row. H must fit in the block's registers + a reduction.
template<int BLOCK>
__global__ void rmsnorm(const __nv_bfloat16* __restrict__ x,
                        const __nv_bfloat16* __restrict__ g,
                        __nv_bfloat16* __restrict__ y, int H, float eps)
{
    const int row = blockIdx.x;
    x += (size_t)row * H;  y += (size_t)row * H;

    float acc = 0.f;
    for (int i = threadIdx.x; i < H; i += BLOCK) {
        float v = __bfloat162float(x[i]);
        acc = fmaf(v, v, acc);                 // accumulate in FP32 — always
    }
    acc = blockReduceSum(acc);

    __shared__ float inv;
    if (threadIdx.x == 0) inv = rsqrtf(acc / H + eps);
    __syncthreads();

    for (int i = threadIdx.x; i < H; i += BLOCK)   // second pass hits L1/L2
        y[i] = __float2bfloat16(__bfloat162float(x[i]) * inv * __bfloat162float(g[i]));
}

@triton.jit
def rmsnorm(X, G, Y, stride, H, eps, BLOCK: tl.constexpr):
    row  = tl.program_id(0)
    cols = tl.arange(0, BLOCK)
    mask = cols < H
    x = tl.load(X + row*stride + cols, mask=mask, other=0.).to(tl.float32)
    inv = tl.math.rsqrt(tl.sum(x*x, axis=0) / H + eps)
    g = tl.load(G + cols, mask=mask, other=0.).to(tl.float32)
    tl.store(Y + row*stride + cols, (x*inv*g).to(tl.bfloat16), mask=mask)

rms = torch.rsqrt(x.float().pow(2).mean(-1, keepdim=True) + eps)
y   = (x.float() * rms * g.float()).to(x.dtype)

4.5.1 If you need the mean too: Welford

LayerNorm needs variance, and the textbook formula Var = E[x²] − E[x]² is numerically dangerous: for data with a large mean and small spread it subtracts two nearly equal large numbers. Welford's algorithm avoids that by updating incrementally:

And it merges in parallel — the merge below is Chan, Golub and LeVeque's update — which is what makes it usable on a GPU: the merge is associative in exact arithmetic, so it goes through the same shuffle tree:

For the large-tensor cases modelled in this part, bandwidth is the limiting roof, so the main structural optimization is to do more per transferred byte. Three levels of it, in increasing order of payoff:

1. Fuse elementwise chains. RMSNorm → SiLU → multiply → residual as one kernel instead of four is a 4× win, as computed in Part 3. This is what torch.compile does automatically, and it is why it helps so much on small models.
2. Fuse into the GEMM epilogue. The accumulator is already in registers when the mainloop ends. Bias, activation, residual add, and quantization all cost a handful of instructions there, versus two full passes over the output tensor as separate kernels. Part 5.
3. Fuse across the algorithm. This is FlashAttention: not "fuse these three elementwise ops" but "restructure so the intermediate never exists". Part 7.

4.6.1 What Triton does and does not do for you

Triton handles shared-memory allocation, swizzled layouts, software pipelines, reduction trees, and tensor-core instruction selection. For a masked tl.load, the mask makes out-of-bounds lanes safe and the compiler can vectorize when the pointer layout, alignment and divisibility permit; coalescing still depends on the layout you supplied and should be checked. You choose block sizes, num_warps, num_stages, and the algorithm. Current Triton APIs can express persistent scheduling and explicit warp-specialized producer/consumer pipelines, but the compiler does not invent that algorithmic restructuring for you. Those schedules remain programmer-visible, which is why the fastest attention kernels still spend so much code on them.

For everything in this part, that gap does not matter. A Triton elementwise, reduction or norm kernel will match a good CUDA one, because the ceiling is memory bandwidth and both reach it. Write these in Triton.

Next: the op at the other end of the roofline. Matmul, from a naive triple loop to something within 10% of cuBLAS, one measured step at a time.

Matmul is the op worth understanding in detail, for three reasons. It is 90–95% of the FLOPs in a transformer. It sits far above the roofline ridge, so it is the one place where arithmetic efficiency is the whole game. And its optimization ladder is the canonical demonstration of every idea in Parts 1–3, in order, each with a measurable payoff.

We are computing C = A B with A of shape M×K, B of K×N, C of M×N. The bare product is 2MNK FLOPs against a compulsory (MK + KN + MN)·s bytes: read A and B, write C. BLAS's more general C ← αAB + βC adds another MN·s read when β ≠ 0. For 4092³ in FP32, the bare product is 137.0 GFLOP against 134 MB read plus 67 MB written = 201 MB; Boehm's BLAS-form benchmark also reads C, for 268 MB total. The A6000's 38.7-TFLOP/s peak gives a 3.54 ms arithmetic floor, while Boehm's measured 23.25-TFLOP/s cuBLAS result takes 5.89 ms; his simplified 30-TFLOP/s estimate is 4.57 ms. Even the larger 268 MB transfer is only about 0.34 ms at peak bandwidth. A good GEMM is decisively compute-bound. Every rung below is a step toward actually achieving that.

__global__ void sgemm_naive(int M, int N, int K, const float* A, const float* B, float* C)
{
    int x = blockIdx.x * blockDim.x + threadIdx.x;   // row of C
    int y = blockIdx.y * blockDim.y + threadIdx.y;   // col of C
    if (x < M && y < N) {
        float acc = 0.f;
        for (int k = 0; k < K; ++k)
            acc += A[x*K + k] * B[k*N + y];
        C[x*N + y] = acc;
    }
}

C = A @ B      # the thing we are trying to match

Problem one: it re-reads everything. Each of the MN threads reads a whole row of A and a whole column of B, so the kernel moves 2MNK elements — an arithmetic intensity of 1 FLOP per element, which Part 3 says is hopeless.

Problem two: half the accesses are uncoalesced. The variable named x is derived from threadIdx.x, so it varies fastest within a warp — and it indexes the row of both A and C. The 32 lanes of a warp therefore read 32 addresses K×4 bytes apart: 32 separate sectors instead of 4.

Fixing the mapping so that consecutive lanes read consecutive columns takes the kernel from 309 to 1,986 GFLOP/s. GMEM throughput goes from 15 GB/s to 110 GB/s. Nothing else changed.

The next rung stages tiles in shared memory. Each block owns a 32×32 output tile and walks K in chunks of 32: load a 32×32 slab of A and of B into SMEM, synchronize, do 32 multiply-adds per thread out of SMEM, synchronize, advance.

__shared__ float As[32*32], Bs[32*32];
float acc = 0.f;
for (int kb = 0; kb < K; kb += 32) {
    int aCol = kb + tx, bRow = kb + ty;
    As[ty*32 + tx] = (row < M && aCol < K) ? A[row*K + aCol] : 0.f;
    Bs[ty*32 + tx] = (bRow < K && col < N) ? B[bRow*N + col] : 0.f;
    __syncthreads();
    for (int k = 0; k < 32; ++k)
        acc += As[ty*32 + k] * Bs[k*32 + tx];
    __syncthreads();                            // before overwriting the tile
}
if (row < M && col < N) C[row*N + col] = acc;

Why so little? Because the arithmetic intensity is still low. With a 32×32 tile each result costs K/16 global accesses, down from 2K — a 32× reduction in DRAM traffic — but each thread now does two shared-memory loads per FMA. The bottleneck simply moved from HBM to SMEM. Nsight shows 66% occupancy, so occupancy is not the problem either.

This is the step that matters. Instead of one accumulator per thread, each thread holds an 8×8 block of accumulators in registers, and the inner loop becomes an outer product: load 8 values of A and 8 of B from shared memory, and do 64 FMAs.

// BM=BN=128, BK=8, TM=TN=8  →  256 threads, each owning an 8×8 tile of C.
float acc[TM][TN] = {0.f};     // 64 registers of accumulator per thread
float regA[TM], regB[TN];

for (int kb = 0; kb < K; kb += BK) {
    loadTilesToSmem(As, Bs, A, B, kb);
    __syncthreads();

    for (int k = 0; k < BK; ++k) {
        // 16 SMEM loads ...
        for (int i = 0; i < TM; ++i) regA[i] = As[k*BM + threadRow*TM + i];
        for (int j = 0; j < TN; ++j) regB[j] = Bs[k*BN + threadCol*TN + j];
        // ... feeding 64 FMAs, all independent -> deep ILP
        for (int i = 0; i < TM; ++i)
            for (int j = 0; j < TN; ++j)
                acc[i][j] += regA[i] * regB[j];
    }
    __syncthreads();
}

# Triton expresses the same thing as a block-level dot product.
# You choose the tile; the compiler chooses the register blocking.
acc = tl.zeros((BLOCK_M, BLOCK_N), dtype=tl.float32)
for k in range(0, tl.cdiv(K, BLOCK_K)):
    a = tl.load(a_ptrs, mask=offs_k[None, :] < K - k*BLOCK_K, other=0.0)
    b = tl.load(b_ptrs, mask=offs_k[:, None] < K - k*BLOCK_K, other=0.0)
    acc = tl.dot(a, b, acc, input_precision="ieee")  # match FP32 SGEMM semantics
    a_ptrs += BLOCK_K * stride_ak
    b_ptrs += BLOCK_K * stride_bk

Notice what this does to the numbers from Part 1. Each thread now uses 64 registers for accumulators alone, plus operand registers and addresses — well over 100. At 128 registers per thread the SM fits 65536 / (32 × 128) = 16 warps: a third of the A6000's 48-warp maximum, a quarter of Hopper's 64. That low occupancy is the deliberate purchase. This is Volkov's tradeoff in its natural habitat: you spend occupancy to buy instruction-level parallelism, and the 64 independent FMAs in the inner loop cover the latency that the missing warps would have covered.

Vectorized access (Part 2, §2.3) is the next 10%: transpose the A tile as it is stored into shared memory so that both the global loads and the shared loads can be 128-bit. Boehm measures 18,237 GFLOP/s, 78.4%, and notes the SMEM transpose alone is worth about 3%.

Warptiling inserts a level between the block tile and the thread tile. It exists because a warp is a real hardware unit: the 32 threads of a warp issue together and their shared-memory accesses are coalesced together, so giving a warp a contiguous rectangle of the output tile makes those accesses regular. The hierarchy becomes:

That is 21,779 GFLOP/s, 93.7% of cuBLAS. It is also exactly the hierarchy CUTLASS uses, and exactly the level at which tensor cores plug in — in Part 6 the innermost level stops being an FMA and becomes an mma.sync or wgmma instruction, but the three enclosing levels are unchanged.

Look again at the structure of the kernel: load tile → barrier → compute → barrier → load next tile. The tensor cores or FMA pipes are idle for the whole load, and the memory system is idle for the whole compute. Perfect overlap changes a serial cost L+C into roughly max(L,C): a gain of up to 2×, reached only when load and compute take equal time.

The fix is double buffering, and then multi-stage pipelining — the structure from Part 2, §2.7. The exact operand path changes with the instruction generation:

• Ampere: cp.async fills circular shared-memory stages; ldmatrix double-buffers register fragments while mma.sync consumes the other set.
• Hopper: TMA fills circular shared-memory stages and the common WGMMA SS form reads A and B there through descriptors, with no ldmatrix stage. The WGMMA RS form can instead keep A in registers while B stays in shared memory.

Two problems remain, and neither is fixed by touching the mainloop.

5.6.1 L2 rasterization

CUDA does not guarantee the order in which blocks execute, and a wave need not be exactly one block per SM. But schedulers commonly draw nearby block IDs into the same period of execution, so the mapping is still a useful locality heuristic. With row-major tile IDs, a wave can span a long row of the output: its blocks reuse one row-panel of A but touch many different column-panels of B. Re-mapping nearby IDs onto a compact 2D region gives both operands a better chance of reuse in L2; measure the result rather than depending on a launch-order guarantee.

# Grouped ordering — five lines, worth >10% on some hardware.
pid = tl.program_id(axis=0)
num_pid_m = tl.cdiv(M, BLOCK_SIZE_M)
num_pid_n = tl.cdiv(N, BLOCK_SIZE_N)
num_pid_in_group = GROUP_SIZE_M * num_pid_n
group_id = pid // num_pid_in_group
first_pid_m = group_id * GROUP_SIZE_M
group_size_m = min(num_pid_m - first_pid_m, GROUP_SIZE_M)
pid_m = first_pid_m + ((pid % num_pid_in_group) % group_size_m)
pid_n = (pid % num_pid_in_group) // group_size_m

// CUTLASS calls this a threadblock swizzle; it maps consecutive
// block indices onto "packed two-dimensional regions" of the output
// to maximise last-level-cache reuse.
//   cutlass/gemm/threadblock/threadblock_swizzle.h

5.6.2 Not enough tiles: split-K and Stream-K

The other problem is the shape of LLM decode. With M = 1 and a 128×256 tile, an 8192×8192 weight matrix produces 32 output tiles for 132 SMs. Three quarters of the GPU is idle, and no inner-loop optimization changes that.

The available parallelism is along K, and it is usable because the partial sums can be reassociated. Because floating-point addition is not associative, that changes rounding; it becomes run-to-run nondeterministic only if the implementation also uses an order that can vary, such as contending atomics. A fixed reduction tree can be deterministic. Split-K partitions the reduction dimension across blocks, each computing a partial product, followed by a reduction kernel. Stream-K generalizes it: instead of assigning tiles to blocks, it assigns an even share of the total inner-loop iterations to a fixed set of blocks, then fixes up the partial sums.

The Stream-K paper (Osama et al., 2023) reports up to 14× over data-parallel CUTLASS with the same blocking, and 6.7× over cuBLAS, measured across 32,824 GEMM shapes — from a single tile configuration per precision. The mechanism is the quantization problem from Part 1. Colfax's worked example: on an H100 PCIe with 114 SMs, going from 114 output tiles to 115 roughly halves device utilization — and Stream-K removes that cliff.

When the mainloop ends, the whole output tile is sitting in registers. Work done there is usually much cheaper than a separate kernel, because it avoids writing 128×128 floats to HBM and reading them back. It is not literally free: extra instructions, registers, shared memory or barriers can reduce throughput and occupancy. Standard fusions:

• bias add, scaling, and the residual add
• activation — GELU, SiLU, and the SwiGLU gate multiply
• quantization back down to FP8/FP4, including computing the per-block scale (Part 8)
• the transpose or layout change needed by the next op

CUTLASS expresses these as epilogue visitor trees so that arbitrary elementwise DAGs can be attached without writing a new kernel. In Triton you just write the arithmetic after the k loop, before tl.store. On Hopper the epilogue usually stages the result back through shared memory with stmatrix so that a TMA store can write it out in one swizzled transfer — Part 6.

A tensor core is a fixed-function matrix-multiply unit. You hand it a small tile of A, a small tile of B and an accumulator, and it produces D = A·B + C in a few cycles. That is the whole idea, and it is worth roughly 15× the FP32 throughput of the same SM.

The complexity is entirely in the interface. The unit does not read "a matrix" — it reads specific bits out of specific registers of specific lanes, or specific bytes at specific shared-memory offsets. Almost every difficult thing in a modern GEMM or attention kernel is a consequence of getting data into that exact shape without wasting bandwidth on the way.

The instruction has been redesigned in every architecture, and the direction of travel is consistent: the operands move further away from the threads, and the issuing scope gets wider. That is not aesthetics. Each generation roughly doubled tensor-core throughput while the register file stayed at 256 KB per SM, so the register file stopped being able to hold or feed the operands.

Start with the one you can hold in your head: mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32. One warp computes a 16×8 output from a 16×16 A and a 16×8 B. Every one of those 128 output elements, and all 384 input elements, lives in a named register of a named lane.

Two derived indices do all the work:

Write g and t for those. Then, straight from PTX §9.7.15.5.8:

Reading a table of index formulas is not the same as seeing it. Pick a lane:

Two structural facts fall out of the picture, and both matter downstream:

• Four consecutive lanes own one row of A, two elements each. Lanes 4g…4g+3 cover columns 0–7 of row g. Four lanes × 2 bf16 elements = 16 bytes. Hold that number.
• The accumulator is already in the right shape to be an A operand for the next matmul — almost. C/D uses the same (g, 2t) pattern as A's first register. This is why FlashAttention can feed P = softmax(S) straight into the second GEMM without a shared-memory round trip.

// One warp: D[16x8] = A[16x16] * B[16x8] + C[16x8], bf16 in, fp32 accumulate.
uint32_t a[4], b[2];      // each holds 2 packed bf16
float    c[4] = {0,0,0,0};

asm volatile(
  "mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32 "
  "{%0,%1,%2,%3}, {%4,%5,%6,%7}, {%8,%9}, {%0,%1,%2,%3};\n"
  : "+f"(c[0]), "+f"(c[1]), "+f"(c[2]), "+f"(c[3])
  : "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]), "r"(b[0]), "r"(b[1]));

Now the problem. Your tile is sitting in shared memory in some sane row-major order. The fragment layout above is not row-major, not column-major, and not any stride pattern a normal load naturally produces — lane 5 needs elements (1,2), (1,3), (9,2), (9,3), (1,10), (1,11), (9,10), (9,11). Ordinary shared-memory loads require multiple strided scalar or vector loads per lane plus a pile of index arithmetic; the exact instruction count depends on their width.

ldmatrix is a warp-collective gather that produces exactly that layout in one instruction:

The calling convention is unusual and worth stating precisely. Each lane supplies an address, and that address is the start of a matrix row, not the lane's own data. For .x4, lanes 0–7 give the eight row addresses of the first 8×8 matrix, lanes 8–15 the second, and so on. The hardware then redistributes: lane t receives row t/4, columns 2(t\,\%\,4) and 2(t\,\%\,4)+1 — identical to the mma.sync A-fragment mapping. That is the entire point of the instruction.

Hopper adds stmatrix (sm_90), the inverse, used in epilogues to push accumulators back into shared memory in a TMA-friendly layout. The matching transpose, movmatrix — an 8×8 register-to-register transposition across the warp, handy in attention backward passes where you need S^T without a shared-memory round trip — is older than people assume: it has existed since Turing (sm_75+, PTX §9.7.15.5.17).

Here is the collision, stated as plainly as possible.

An ldmatrix.m8n8.b16 operand is 8 rows × 16 bytes. Store your tile row-major with a K dimension of 64 bf16 values, and each row is 128 bytes, so the eight addresses the warp supplies are 128 bytes apart. Bank = (byte / 4) mod 32, and 128 bytes is exactly 32 banks — so all eight addresses land on the same four banks. An 8-way conflict, on every single operand load, for the rest of your kernel's life.

Padding can rotate the bank mapping, but it does not inherently destroy alignment: adding 16 bytes makes the stride 144 bytes, still a multiple of 16. Smaller arbitrary padding can violate the alignment rules, and any padding costs footprint and complicates the descriptor/layout contract. Production kernels therefore usually permute the layout with an XOR that preserves every 16-byte chunk and matches the layouts TMA and WGMMA describe. From Part 2:

Every mode that TMA and wgmma accept uses M = 4 and S = 3: the low 16 bytes are never permuted and the pattern repeats every 8 rows. So the 128-byte mode reduces to

and the eight addresses of an ldmatrix operand land in eight different bank groups. The figure below is the same explorer as Part 2, repeated here because this is where it earns its keep:

Three consequences worth memorising:

• The 8×16 B no-swizzle atom is exactly one ldmatrix.m8n8 operand. The swizzled atoms are the same 8 rows widened to 32, 64 or 128 bytes.
• Prefer 128 B swizzle when the tile permits it. Its atom covers one 128-byte contiguous region on the global side and a full conflict-free warp access on the shared side.
• But 128 B swizzle caps the innermost TMA box at 128 bytes — 64 elements in bf16. A BLOCK_K = 128 bf16 K-major tile needs two TMA boxes, or a narrower swizzle. Violating this rule is one common cause of CUDA_ERROR_INVALID_VALUE from cuTensorMapEncodeTiled.

wgmma.mma_async changes three things at once. It is issued by 128 threads (four contiguous warps whose first warp-rank is a multiple of 4). It is asynchronous — it returns immediately and you wait on a group. And its B operand, and optionally its A operand, are read directly from shared memory by the tensor core, with no ldmatrix and no registers involved.

What you pass instead of a register is a 64-bit matrix descriptor: a packed structure holding the shared-memory base address, two strides, and the swizzle mode. The hardware walks the tile itself.

// SS form: both operands described by shared-memory descriptors.
wgmma.fence.sync.aligned;                       // order prior writes vs. the async read
fence.proxy.async;                              // generic proxy -> async proxy

wgmma.mma_async.sync.aligned.m64n256k16.f32.bf16.bf16
      {d0,...,d127}, a-desc, b-desc, 1, 1, 1, 0, 0;
//                       scale-d ^  ^ scale-a  ^ scale-b  ^ trans-a ^ trans-b

wgmma.commit_group.sync.aligned;                // close the batch
wgmma.wait_group.sync.aligned 1;                // keep 1 group in flight

// The layout atoms CuTe uses to keep SMEM indexing and the descriptor consistent.
using SmemLayoutAtom = GMMA::Layout_K_SW128_Atom<bfloat16_t>;
auto smem_layout = tile_to_shape(SmemLayoutAtom{}, Shape<_128,_64>{});

// Swizzle<3,4,3> is the 128B mode: 8 rows x 128 B, chunk ^= row.
//   include/cute/atom/mma_traits_sm90_gmma.hpp

The accumulator cost is worth noting: for m64nNk16 with FP32 accumulate, each of the 128 threads holds N/2 registers. At N=256 that is 128 registers per thread for the accumulator alone — half the architectural maximum. That is the pressure that eventually pushed the accumulator out of the register file entirely.

If the tensor core reads its operands straight out of shared memory, then the kernel's remaining job is to keep shared memory full. Hopper's answer is the Tensor Memory Accelerator: a single thread issues one instruction describing a multidimensional tile copy, and dedicated hardware generates every address, detects out-of-bounds elements, applies the swizzle, and signals completion. Out-of-bounds loads are filled with zero or the configured OOB-NaN value; stores are discarded rather than clamped to an edge element.

The descriptor is built on the host:

CUtensorMap map;
cuTensorMapEncodeTiled(&map,
    CU_TENSOR_MAP_DATA_TYPE_BFLOAT16, /*rank*/ 2, A_gmem,
    globalDim, globalStrides,            // strides: multiples of 16 B, < 2^40
    boxDim, elementStrides,              // boxDim[0]*sizeof(elem) <= swizzle size
    CU_TENSOR_MAP_INTERLEAVE_NONE,
    CU_TENSOR_MAP_SWIZZLE_128B,          // must match the SMEM layout atom
    CU_TENSOR_MAP_L2_PROMOTION_L2_128B,
    CU_TENSOR_MAP_FLOAT_OOB_FILL_NONE);
// Pass it into the kernel as: const __grid_constant__ CUtensorMap map

// One elected thread arms the barrier and issues the copy.
if (elect_one()) {
  mbarrier.arrive.expect_tx.shared::cta.b64 _, [bar], TILE_BYTES;
  cp.async.bulk.tensor.2d.shared::cluster.global.mbarrier::complete_tx::bytes
      [smem], [tensorMap, {k_off, m_off}], [bar];
}
// Everyone waits on the phase bit. try_wait may return early, hence the loop.
wait: mbarrier.try_wait.parity.shared::cta.b64 p, [bar], phase;
      @!p bra wait;

auto tma = make_tma_copy(SM90_TMA_LOAD{}, gmem_tensor, smem_layout);
// ... in the kernel
copy(tma.with(barrier, mcast_mask), tAgA(_,_,k), tAsA(_,_,stage));

6.6.1 mbarriers count bytes, not just threads

A __syncthreads() counts arrivals. An mbarrier also counts transactions — bytes that asynchronous engines have promised to deliver. The handshake is:

1. A thread does mbarrier.arrive.expect_tx [bar], TILE_BYTES, which both arrives and adds TILE_BYTES to the barrier's transaction count.
2. The TMA engine decrements that count as bytes land in shared memory.
3. The barrier flips phase only when both the pending arrivals and the transaction count reach zero (PTX §9.7.14.16.8).

Consumers wait on the phase parity — a single bit that alternates each time the barrier completes. For a circular buffer of S stages, the barrier at stage s has completed k/S phases by iteration k, so the parity to wait for is (k/S) & 1. That one line is most of the bookkeeping in a Hopper mainloop.

TMA has two further tricks that matter for GEMM. Multicast: with a cluster (Part 2, §2.8), one cp.async.bulk.tensor….multicast::cluster delivers the same tile into the shared memory of every CTA named in a 16-bit mask — one HBM read feeding several SMs. And im2col mode, which unrolls convolution windows during the copy, so a convolution can use a plain GEMM mainloop.

Once loads are issued by one thread and the math is issued by a warpgroup, it stops making sense for every warp to do both. Hopper kernels split the block into producer and consumer warpgroups: the producer does nothing but issue TMA copies and manage barriers; the consumers do nothing but wgmma and the epilogue.

The producer needs almost no registers. The consumers need as many as they can get. So the ISA lets a warpgroup hand registers back:

The register pool is per-CTA. .inc blocks until enough registers are free, the value must be between 24 and 256 and a multiple of 8, and every warp of the warpgroup must execute the same instruction (PTX §9.7.20.5). Done right, a producer warpgroup costs you 4 warps' worth of slots and almost no register file.

The ping-pong part is specific to attention and is covered in Part 7; for a plain GEMM the two consumer warpgroups simply work on different halves of the output tile. FlashAttention-3 measures each contribution separately: 582 TFLOP/s with warp specialization but no GEMM–softmax pipelining, 570 with pipelining but no warp specialization, and 661 with both.

tcgen05.mma takes the trend to its conclusion. It is issued by a single thread. B always comes from shared memory; A comes from shared memory or from Tensor Memory; and the accumulator lives only in Tensor Memory. The MMA itself is asynchronous, like everything else at this generation: completion is tracked on an mbarrier by tcgen05.commit (PTX §9.7.17.12.1) — the same handshake TMA uses.

TMEM is a dedicated 256 KB scratchpad per SM, organised as 128 lanes × 512 columns of 32-bit cells. It is allocated dynamically at runtime by one warp, in units of 32 columns, power-of-two sized, and it must be explicitly freed before the kernel exits. Access is restricted by warp: warp 0 of a warpgroup can only touch lanes 0–31, warp 1 lanes 32–63, and so on — so reading back a full accumulator requires a whole warpgroup.

Assembling everything: a Hopper GEMM mainloop, one iteration, warp-specialized, with a multi-stage circular buffer.

// ---- setup ----
// SMEM: S >= 2 stages x (A tile + B tile), laid out with Layout_K_SW128_Atom
// mbarriers: full[S] (producer -> consumer), empty[S] (consumer -> producer)

if (warpgroup_id == PRODUCER) {
    setmaxnreg_dec(40);
    for (k = 0; k < num_k_tiles; ++k) {
        s = k % S;
        wait(empty[s], parity_of(k/S - 1));       // stage free?
        if (elect_one()) {
            expect_tx(full[s], TILE_BYTES);
            tma_load(A_smem[s], tensor_map_A, {k*BK, m0}, full[s]);
            tma_load(B_smem[s], tensor_map_B, {k*BK, n0}, full[s]);
        }
    }
} else {                                        // CONSUMER
    setmaxnreg_inc(232);
    wgmma_fence();                                  // once, after ordinary acc/A-reg writes
    for (k = 0; k < num_k_tiles; ++k) {
        s = k % S;
        wait(full[s], parity_of(k/S));               // bytes landed?
        wgmma_mma_async(acc, desc(A_smem[s]), desc(B_smem[s]));
        wgmma_commit_group();
        if (k > 0) {
            wgmma_wait_group<1>();                   // group k-1 is complete
            old = (k - 1) % S;
            arrive(empty[old]);                      // only now may TMA overwrite it
        }
    }
    wgmma_wait_group<0>();                       // drain the last group before reading acc
    if (num_k_tiles > 0)
        arrive(empty[(num_k_tiles - 1) % S]);      // release the final operand stage
    epilogue(acc);        // stmatrix -> SMEM -> TMA store
}

# Compact algorithmic form; it does not contractually select TMA or warp specialization.
@triton.autotune(configs=[
    triton.Config({'BLOCK_M':128,'BLOCK_N':256,'BLOCK_K':64,'GROUP_M':8},
                  num_stages=3, num_warps=8),
], key=['M','N','K'])
@triton.jit
def matmul(...):
    ...
    acc = tl.zeros((BLOCK_M, BLOCK_N), dtype=tl.float32)
    for k in range(0, tl.cdiv(K, BLOCK_K)):
        a = tl.load(a_ptrs); b = tl.load(b_ptrs)
        acc = tl.dot(a, b, acc)
        a_ptrs += BLOCK_K * stride_ak; b_ptrs += BLOCK_K * stride_bk

6.9.1 So when do you need the CUDA version?

For a plain GEMM: almost never. Triton's generated Hopper kernels are competitive, and cuBLAS/CUTLASS is right there.

You reach for the low level when the schedule is the optimization — when you need a specific producer/consumer split, a ping-pong between two consumer warpgroups, an in-kernel FP8 transpose, or TMEM allocated in a particular shape. That is the FlashAttention-3 territory in Part 7, and it is the reason the fastest attention kernels are still CUTLASS.

Attention is where every idea in this series meets at once. The online-softmax recurrence from Part 4, the tiling from Part 5, the tensor-core layouts and warp specialization from Part 6 — FlashAttention is those four things assembled into one kernel. And then, at inference time, it turns into something else entirely.

The operation itself is three lines:

Count the arithmetic first. Per head, forward: 2N²d for QK^T plus 2N²d for PV = 4N²d FLOPs. The compulsory memory is Q, K, V, O — four N×d tensors. So the arithmetic intensity is

Seven times above the H100 ridge point. Attention should be comfortably compute-bound.

It is not, in the naive implementation, because S and P are materialised. Those are N×N matrices, written to HBM and read back — twice, in most implementations, once for the softmax and once for the second matmul. And N² grows a great deal faster than Nd.

The obstacle is that softmax needs the whole row before it can normalize anything, and the whole row is exactly what we are refusing to materialise. Part 4 already dissolved that: the online-softmax recurrence lets you consume a row in pieces, keeping only a running maximum m and a running sum ℓ.

FlashAttention adds one thing: the output accumulator has to be rescaled too. Here is the forward pass, verbatim in structure from FlashAttention-2, Algorithm 1:

That bolded line is the entire difference from Part 4. Whenever a new block reveals a larger maximum, the accumulated output is scaled by the same factor as the accumulated sum, so the two stay consistent. Both factors are e to a non-positive power, so neither can overflow.

Two FA2-specific refinements are visible in the last two lines. Dividing by ℓ is deferred out of the inner loop, because a division is a non-matmul FLOP and those are 16× more expensive than matmul FLOPs on an A100. And only the single scalar L = m + \log ℓ is stored for the backward pass, not both m and ℓ.

# The inner loop, in the shape Triton's tutorial 06 uses.
for start_n in range(0, hi, BLOCK_N):
    k = tl.load(K_block_ptr)
    qk = tl.dot(q, k) * qk_scale                     # S tile — stays in registers
    if CAUSAL:
        qk = tl.where(offs_m[:, None] >= (start_n + offs_n[None, :]), qk, -1.0e6)

    m_new = tl.maximum(m_i, tl.max(qk, 1))            # running max
    qk    = qk - m_new[:, None]
    p     = tl.math.exp2(qk)                         # exp2, not exp — see Part 4
    alpha = tl.math.exp2(m_i - m_new)                # THE rescale factor

    l_i = l_i * alpha + tl.sum(p, 1)
    acc = acc * alpha[:, None]                       # rescale the accumulator
    v   = tl.load(V_block_ptr)
    acc = tl.dot(p.to(v.dtype), v, acc)              # second GEMM, straight from registers
    m_i = m_new
    K_block_ptr = tl.advance(K_block_ptr, (0, BLOCK_N))
    V_block_ptr = tl.advance(V_block_ptr, (BLOCK_N, 0))

acc = acc / l_i[:, None]                             # normalize once

# What it must equal. Note this allocates the N×N matrix.
S = (q @ k.transpose(-2, -1)) * (1.0 / math.sqrt(d))
if causal:
    S = S.masked_fill(mask, float('-inf'))
P = torch.softmax(S, dim=-1)
O = P @ v

Notice what does not appear: any store of p to HBM. The tile stays on chip and feeds the second GEMM; the compiler and kernel handle whatever register-layout conversion the chosen dtype and MMA instruction require. That qualification matters: FlashAttention-3's FP8 wgmma path, for example, explicitly permutes values because its accumulator and A-operand layouts differ.

The 2022 version proved the idea and reached 25–40% of an A100's peak. Getting from there to 75% of an H100's peak took two more papers, and each one is a lesson from an earlier part of this series.

7.3.1 FA2: parallelism and warp partitioning

• Fewer non-matmul FLOPs. On an A100, 312 TFLOP/s of matmul against 19.5 TFLOP/s of FP32 non-matmul — every non-matmul FLOP costs 16×. Hence deferring the division and storing only L.
• Parallelize over sequence length, not just over batch×heads. With long sequences and small batches there simply were not enough thread blocks to fill the GPU — the wave-quantization problem from Part 1.
• Split-Q instead of split-K. FA1 split K and V across the four warps of a block, so each warp computed a slice of QK^T and then had to write intermediates to shared memory, synchronize, and sum. FA2 splits Q across warps and keeps K and V shared. In the paper's words: "There is no need for communication between warps."

Result: about 2× over FA1, reaching 50–73% of theoretical peak on an A100 (FA1 managed 25–40%).

7.3.2 FA3: asynchrony, and the exponential problem

On Hopper the bottleneck moves. From Part 4: the SFU produces 3.9 TFLOP/s of exp against 989.5 TFLOP/s of BF16 matmul — a 254× gap, which the paper rounds to 256×. In an attention forward pass with head dimension 128 there are roughly 512× more matmul FLOPs than exponentials, so the exponentials take about half as long as the matmuls themselves. They are on a different execution unit, so the fix is to run them concurrently.

• Producer/consumer warp specialization with TMA and wgmma, exactly the structure in Part 6 §6.7. The producer warpgroup drops to 40 registers and issues TMA loads into an s-stage circular buffer; consumer warpgroups take 232 registers and run the two GEMMs.
• Ping-pong scheduling. bar.sync instructions force warpgroup 1's GEMMs to be scheduled ahead of warpgroup 2's, so WG1's softmax runs during WG2's GEMMs and then the roles swap. Measured: 570 → 620–640 TFLOP/s for FP16, head dim 128, sequence 8192. A second overlap runs inside each consumer warpgroup: a 2-stage register pipeline issues block j+1's QK^T while block j's softmax is still in flight — that is the "GEMM–softmax pipelining" row of the ablation below.
• FP8, with two accuracy devices: per-block quantization (one scale per B_r×d block, fused into the preceding rotary embedding for free), and incoherent processing — multiply Q and K by a random orthogonal matrix M before quantizing. Since MM^T = I, (QM)(KM)^T = QK^T, so the result is unchanged, but outliers get spread across dimensions.

Peak reported by FlashAttention-3: 740 TFLOP/s FP16 (75% of theoretical) and close to 1.2 PFLOP/s in FP8. Speedup over FA2: 1.5–2.0× forward, 1.5–1.75× backward.

With a causal mask, query i only attends to keys j ≤ i. In a tiled kernel that means roughly half the (i,j) tile pairs are entirely masked out, and you can skip computing them altogether — not mask them, skip them. Only the tiles on the diagonal need an actual elementwise mask applied.

FA2 measures this at 1.7–1.8× versus non-causal. It also creates a load-imbalance problem: block i does work proportional to i, so a naive assignment leaves half the SMs finishing early. Production kernels reverse or interleave the block order so that long and short rows are paired.

At inference, after the prompt has been processed, each step has one query row per sequence. The arithmetic changes character completely:

Decode attention is not a matmul problem. It is a KV-cache read, at an arithmetic intensity of 1 against a ridge point of 295. And FlashAttention's parallelism — batch × heads × query blocks — loses its third factor entirely.

The fix is the same one as split-K in Part 5: find parallelism along the reduction dimension. Flash-Decoding splits K and V along the sequence, runs an independent FlashAttention against each chunk, stores one extra log-sum-exp scalar per chunk, and then combines them — using the very same (m, ℓ) merge operator from Part 4, applied one level up.

Flash-Decoding reports near-constant time up to ~32k tokens, and up to 8× end-to-end on long-context decoding — with the attention step itself up to 50× faster than FlashAttention. Note the middle column at 131k: FlashAttention-2 is slower than eager PyTorch there, because it has traded parallelism for memory efficiency in exactly the regime where parallelism is what you need.

Everything a decode step reads is the KV cache, so its size is the whole story:

For Llama-3-70B (80 layers, 8 KV heads, head dim 128, BF16) that is 320 KiB per token. At 8k context and batch 32 the cache alone is 85.9 GB — more than an 80 GB H100 — and every decode step must re-read all of it on top of the 140 GB of weights.

7.6.1 Grouped-query attention, and what batching cannot amortise

With g query heads sharing one KV head, each KV element is loaded once and consumed g times, so the intensity of decode attention is

All still one to two orders of magnitude below the ridge point. And here is the asymmetry that shapes serving systems:

7.6.2 PagedAttention

The second KV problem is not bandwidth but fragmentation. A request's cache grows unpredictably, so a serving system that pre-allocates a contiguous maximum-length buffer wastes most of it. The vLLM paper (Kwon et al., 2023) measured the fraction of allocated KV memory actually holding token state: 20.4–38.2% for existing max-length allocators, 96.3% for PagedAttention.

The mechanism is virtual memory, applied to the cache: split it into fixed-size blocks, keep a per-request block table mapping logical to physical blocks, and let the blocks be scattered. Because every block has the same size there is no external fragmentation at all, and the waste that remains is internal — the unfilled slots of one final block per request, under 4% in the paper's measurements. Blocks can then be shared — between beam-search branches, parallel samples, or requests with a common prefix. The reported end-to-end gain is 2–4× throughput at the same latency.

The kernel pays for it with an extra indirection: every K/V access goes through the block table. vLLM's cache layout is shaped around keeping the resulting loads 16-byte aligned:

That trailing x exists purely so a thread group's 16-byte vector load lands contiguously — the rule from Part 2 §2.3, reappearing as a tensor layout decision.

• Training / prefill: use FlashAttention-4 on supported Hopper and Blackwell systems, with FA3 or FA2 as compatibility fallbacks. Alternatively use F.scaled_dot_product_attention, which selects among PyTorch's documented backends, including FlashAttention-2 where eligible; it does not promise to dispatch to external FA3 or FA4 packages.
• Decode: a split-KV kernel — Flash-Decoding, vLLM's paged kernel, or SGLang's. The right choice depends on your batching and cache-sharing strategy more than on the kernel.
• Writing your own: reasonable in Triton if you have an unusual mask, a custom score modification, or a fused pre/post-op. Tutorial 06 is a working starting point. On standard shapes, compare against the current FlashAttention implementation before claiming a win.

Next: the last big lever on both of those — using fewer bits.

Every other technique in this series moves a workload along the roofline. Low precision moves the roofline itself. On an H100, going from BF16 to FP8 halves the bytes a decode step must read and doubles the peak FLOP/s — 989.5 to 1979 TFLOP/s, dense. On a B200, FP4 does it again: 9 PFLOPS dense (18 with sparsity) against 4.5 for FP8.

The catch is that the numbers still have to be right. A BF16 value has 8 exponent bits and covers the same range as FP32; an FP8 E4M3 value has 4, and tops out at 448. One activation outlier can make an entire tensor unrepresentable. Everything below is a way of shrinking the group of numbers that has to share a single scale.

All of these are the same idea — a sign bit, some exponent bits, some mantissa bits — with the budget split differently. Exponent bits buy range; mantissa bits buy precision.

An FP8 E4M3 tensor can only hold values up to 448, so before you cast anything you multiply by a scale chosen so the largest element lands near the top of the range. The GEMM then runs in FP8 and exposes an FP32 accumulator/output interface; on Hopper, however, the tensor core retains fewer bits internally during FP8 accumulation, as §8.4 shows. The result is divided by the product of the two input scales.

// Quantize two FP32 values to E4M3 and back, with a per-tensor scale.
// s = amax / 448 comes out of an amax reduction; here we pass s and 1/s.
// sm_89+ (PTX ISA 8.1). .rn = round to nearest even; .satfinite = clamp
// to +-448 on overflow instead of producing inf -- E4M3 has none.
__device__ float2 qdq_e4m3(float2 v, float s, float s_inv)
{
    unsigned short q;                    // two packed e4m3 values
    // First source lands in the HIGH half, so the operands are swapped
    // to keep v.x in the low byte.
    asm("cvt.rn.satfinite.e4m3x2.f32 %0, %2, %1;"
        : "=h"(q) : "f"(v.x*s_inv), "f"(v.y*s_inv));
    unsigned int p;
    asm("cvt.rn.f16x2.e4m3x2 %0, %1;" : "=r"(p) : "h"(q));   // exact in f16
    float2 r = __half22float2(*reinterpret_cast<__half2*>(&p));
    r.x *= s;  r.y *= s;                 // the dequantization multiply
    return r;
}

# The same round trip in Triton: the compiler emits the same cvt
# instructions, and casts to fp8 saturate to +-448 — unlike eager
# PyTorch, whose float8_e4m3fn cast produces NaN on overflow.
xq = (x * s_inv).to(tl.float8e4nv)     # quantize: E4M3, clamped
xd = xq.to(tl.float32) * s             # dequantize, never leaves registers

That works well until the tensor contains an outlier. LLM activations reliably do: a handful of channels carry values one or two orders of magnitude larger than the rest. Scale for the outlier and many ordinary values may underflow or collapse into the bottom few representable values.

The central scaling question is therefore: how few numbers can be made to share one scale?

NVIDIA's taxonomy of the options, in increasing order of granularity:

• Per-tensor, delayed — the scale comes from a rolling history of amax. Smooths transient outliers, but a large one dominates the history for as long as the window lasts.
• Per-tensor, current — the scale comes from this pass's amax. Adapts instantly, no historical stability.
• Generic block scaling — user-chosen 1×128 or 128×128 blocks with FP32 scales. This is what DeepSeek-V3 uses.
• Microscaling — 32 elements per scale, in hardware. NVIDIA's MXFP8 pre-training recipes report less than 0.5% difference in validation perplexity against BF16 across a full pre-training run, on 843M and 8B dense models.

The OCP Microscaling (MX) specification — written jointly by AMD, Arm, Intel, Meta, Microsoft, NVIDIA and Qualcomm — fixes the block structure so hardware can implement it: 32 consecutive elements plus one shared scale in UE8M0 — eight exponent bits, no mantissa, so the scale is a pure power of two and applying it is an exponent add.

NVIDIA's Blackwell-native NVFP4 tightens both knobs: blocks of 16 instead of 32, and a two-level scale — an E4M3 microscale per block, which unlike a power of two can express fractions, plus a single FP32 scale per tensor that keeps those microscales inside E4M3's own range.

On Blackwell the grouping and the scale application happen inside the tensor core — you supply a scale-factor tensor in its required tiled layout alongside the data, and tcgen05.mma.block_scale does the rest. On Hopper there is no hardware block scaling, which is why DeepGEMM's README notes that SM90 wants FP32 scale factors while SM100 wants packed UE8M0.

NVIDIA reports accuracy degradation of 1% or less across seven evaluations when DeepSeek-R1-0528 is requantized from FP8 to NVFP4.

Here is a failure mode that is easy to miss, and DeepSeek-V3 documented it precisely: the Hopper FP8 tensor core does not accumulate in full FP32 internally.

Their fix is a nice illustration of using two units at once: every 128 elements of K — four wgmma instructions — the partial result is copied out to FP32 registers and accumulated on the CUDA cores at full precision. In this schedule, the dequantization multiply can overlap work on CUDA cores with tensor-core MMA. And because two wgmma operations are typically in flight, one warpgroup can be doing the promotion while the other is doing its MMA.

8.4.1 The rest of DeepSeek-V3's recipe

DeepSeek's other choices are worth recording because they contradict the usual advice:

• E4M3 everywhere, including gradients — not the conventional E4M3-forward / E5M2-backward split. Their argument: fine-grained grouping already shares exponent bits among small groups, so the extra exponent bit is not needed.
• Different granularity for activations and weights: 1×128 tiles for activations (per token, per 128 channels), 128×128 blocks for weights.
• These are not interchangeable. Appendix B.2 reports that 128×128 block quantization of activation gradients destabilises training — the backward Dgrad computation is specifically sensitive. Use tiles there.

Measured outcome: relative training-loss error below 0.25% at 16B/1.33T tokens and 230B/0.9T tokens, and DeepGEMM reaching up to 1550 FP8 TFLOP/s on an H800.

There are two fundamentally different reasons to quantize, and they help in opposite regimes.

Weight-only (W4A16 — GPTQ, AWQ) stores weights in 4 bits and dequantizes them in-kernel back to FP16 before the matmul. The tensor cores still run at FP16 rate. That sounds pointless until you remember Part 3: at small batch the decode step can be bandwidth-bound, so 4-bit weights mean up to 4× less weight traffic. Its dequantization cost can be hidden while arithmetic resources and tensor-core throughput have headroom; it is not literally free in every kernel.

W8A8 (FP8 or INT8) quantizes activations too, so the tensor cores genuinely run at 2× the FP16 rate. That is what helps prefill and large-batch decode. Its cost is activation outliers — exactly the problem §8.2 is about.

AWQ's central observation is worth stating on its own: "protecting only 1% salient weights can greatly reduce quantization error" — and the salient channels are identified by activation magnitude, not weight magnitude. Since it needs no backpropagation or reconstruction, it avoids the reconstruction-overfitting failure mode the paper contrasts with GPTQ.

1. Prefer FP32 accumulation for accuracy-sensitive reductions wherever the hardware path provides it, regardless of the storage type. At FP8 on Hopper, remember that an FP32 accumulator interface does not guarantee full internal FP32 accumulation; periodic promotion is one measured remedy.
2. Quantize the KV cache before you quantize the weights, if you serve long contexts. Part 7 showed KV traffic overtaking weight traffic; FP8 KV halves the dominant term.
3. Pick the granularity from the tensor, not from a default. Weights are well-behaved and tolerate coarse blocks; activations have outliers and want 1×128 tiles or finer; activation gradients want tiles specifically.
4. Fuse quantization into the producing kernel's epilogue. Computing a per-block amax and scaling costs a reduction over data that is already in registers. As a separate pass it costs two full trips to HBM.
5. Measure end-task accuracy, not RMSE. The FlashAttention-3 ablation in Part 7 is the cautionary example: two techniques were credited for the FP8 accuracy result, and the ablation shows one of them was doing essentially all the work.

One part left: how to tell whether any of this actually helped.

The single most common failure in kernel work is optimizing something that was never the bottleneck. Parts 1–8 gave you mechanisms; this part is about pointing them at the right thing, and about knowing when to stop.

"This kernel takes 400 µs" is not information. "This kernel moves 268 MB, so it cannot beat 80 µs, and it takes 400" is a bug report with a number attached to it. From Part 3, the procedure is:

1. Count the FLOPs and the compulsory bytes.
2. Compute both ceilings: bytes / bandwidth and FLOPs / peak.
3. The larger is your ideal runtime. If measured throughput reaches roughly 85% of the relevant bandwidth or compute ceiling — equivalently, runtime is no more than ideal time / 0.85 — stop.
4. If you are not, the gap has a cause, and Nsight will name it.

Nsight Compute exposes thousands of counters. About a dozen matter for the kernels in this series, and each maps to something you have already seen.

Two of them deserve special attention because they are the direct measurements of Part 2's two rules:

• Sectors per request. A naturally aligned, coalesced 32-lane, 4-byte load touches 4 32-byte sectors; misalignment can add another. If this instruction reports 32, it generated 8× as many sector transactions as the aligned ideal. That is strong evidence of a lane-to-address problem, not a promise of an 8× whole-kernel speedup or a one-line fix.
• Shared-memory bank conflicts. Zero is the ideal for ordinary loads and stores. A non-zero count is evidence to inspect the stride and instruction layout; padding, swizzling, or a different access mapping are possible remedies, and the right one depends on the instruction.

Two or three profiling passes cover the whole list:

# Full section set on one launch of the kernel, after warm-up — the first
# launches are distorted by autotuning, caching and clock ramp-up:

ncu --set full -k regex:sgemm --launch-skip 2 --launch-count 1 ./app

# Targeted pass when you already know what you are looking for:
# one replay instead of dozens, seconds instead of minutes.

ncu --metrics gpu__dram_throughput.avg.pct_of_peak_sustained_elapsed,\
l1tex__average_t_sectors_per_request_pipe_lsu_mem_global_op_ld.ratio ./app

This is the highest-information view in the profiler. The Warp State Statistics section counts, per issued instruction, how many cycles the warps spent unable to issue — attributed by reason. Each reason supplies a hypothesis; correlate it with the responsible instruction and the memory or pipeline metrics before choosing a remedy.

9.3.1 Reading a breakdown: Part 5's naive matmul

Profile the first kernel of Part 5 — 4092³, 309 GFLOP/s, 1.3% of cuBLAS — with --set full. Speed Of Light shows both headline throughputs in low single digits: nothing is saturated, so this is not yet a bandwidth or arithmetic problem. Achieved occupancy sits close to its theoretical limit, so warp supply is not it either (Part 1). Then Warp State Statistics names the disease: Long Scoreboard dwarfs every other reason, and Sectors per Request reads 32, not 4.

For this kernel, source correlation gives the pair one clear reading: warps wait on the global loads whose column-major walk down B generates eight times the ideal sector traffic. Remap threads so lanes read consecutive columns. The profiler confirms the diagnosis — sectors drop to 4, and the same kernel now measures 1,986 GFLOP/s. Re-profile before celebrating: Long Scoreboard no longer stands alone, while Short Scoreboard and Wait suggest inspecting the MIO dependencies and instruction-level parallelism; the next measured rung stages tiles through shared memory. Form a hypothesis from the top reason, confirm it against the responsible instruction and companion counters, change one thing, and re-measure.

Part 5 had the data point: the best GEMM configuration on an A6000 is BM=BN=128, BK=16, TM=TN=8; on an A100 it is BM=BN=64, BK=16, TM=TN=4, and using the A6000 numbers there costs 6%. The search space is a few hundred valid combinations and the objective is not smooth — occupancy is a staircase, wave quantization is a sawtooth, and shared memory is a cliff.

So you search. In Triton this is one decorator:

@triton.autotune(
    configs=[
        triton.Config({'BLOCK_M':128,'BLOCK_N':256,'BLOCK_K':64,'GROUP_M':8}, num_stages=3, num_warps=8),
        triton.Config({'BLOCK_M':64, 'BLOCK_N':256,'BLOCK_K':32,'GROUP_M':8}, num_stages=4, num_warps=4),
        triton.Config({'BLOCK_M':128,'BLOCK_N':128,'BLOCK_K':32,'GROUP_M':8}, num_stages=4, num_warps=4),
        # ... a few dozen more
    ],
    key=['M', 'N', 'K'],      # re-tune when the shape changes
)
@triton.jit
def matmul_kernel(...):
    ...

# The knobs, and what each one trades off:

BLOCK_M, BLOCK_N   # arithmetic intensity (Part 3) vs registers and
                   # wave quantization. Square often balances reuse for
                   # square GEMMs; rectangular shapes can favor rectangles.
BLOCK_K            # shared-memory footprint per stage → pipeline depth.
                   # Does NOT affect intensity — it cancels out.
GROUP_M            # L2 rasterization order (Part 5, §5.6.1)
num_warps          # threads per block → occupancy staircase (Part 1)
num_stages         # circular buffer depth (Part 2, §2.7).
                   # Fit the device's per-block limit and target residency.

Three practical warnings. First, key matters: if your shapes vary continuously, you will re-tune constantly. Bucket them. Second, autotuning finds the best point in the space you gave it — if every config in your list has the same num_stages, no amount of searching will discover pipelining. Third, Triton runs the kernel repeatedly for every configuration: if it mutates an input, use reset_to_zero, restore_value, or hooks so trials do not accumulate state. And when the space is CUTLASS-sized rather than decorator-sized, use the tool that ships with it: cutlass_profiler sweeps CUTLASS's whole GEMM zoo over your shapes and reports the winner, which is how you pick a tile to pin in code.

Part 3's third regime is invisible on a roofline and dominates small-model, small-batch inference. Every kernel launch costs a few microseconds of CPU and driver work. Twenty small ops per transformer layer, each running for 3 µs and costing 5 µs to launch, spends most of its wall clock in Python and the driver — and the GPU timeline is mostly gaps.

The three tools, in the order you should reach for them:

1. torch.compile — fuses adjacent elementwise work automatically, which removes launches and DRAM round trips (Part 3). Usually the largest single win on an unoptimized model, for the least effort.
2. CUDA Graphs — record a whole sequence of launches once and replay it as a single submission. Removes per-launch CPU cost for a fixed sequence of kernels with fixed shapes. This is what makes small-batch decode viable; it is also why serving frameworks pad batch sizes to a fixed set of buckets.
3. Hand fusion — write the fused kernel. Reach for this when the fusion crosses a boundary the compiler will not cross, such as into a GEMM epilogue.

An honest list, because the answer is "usually don't":

• Plain GEMM. cuBLAS and CUTLASS have absorbed more engineering than you can match. Boehm's 94% took six weekends and was a learning exercise, not a replacement.
• Standard attention. FlashAttention-4 is the current implementation for supported Hopper and Blackwell systems; FA3's historical 75% of H100 peak explains the engineering in Part 7. Use the maintained implementation that supports your hardware and dtype.
• Anything torch.compile will fuse. Try it first; it is one line and it handles the common cases.

The cases where writing one is right:

• A fusion that crosses a library boundary — quantize-and-store in a GEMM epilogue, rotary embedding fused into the attention prologue, an optimizer step fused into the gradient reduction.
• A shape or mask the library does not have — unusual attention biases, MoE routing, sparsity patterns, custom score modifications.
• A schedule that is the optimization — a specific producer/consumer split, a particular TMEM allocation. That is the CUTLASS-level work in Part 6, and it is where the remaining performance actually lives.
• Understanding. Which is most of why you read this.

And the practical middle ground: write it in Triton first. For simple, coalesced streaming kernels, Triton often approaches the same bandwidth ceiling as hand-written CUDA. Memory-bound does not mean every implementation reaches that ceiling: layout, cache behavior, irregular accesses and scheduling still matter. Drop to CUDA when measurement shows that lower-level schedule or layout control is the missing piece.

Nine parts, compressed into the question each one answers.

Everything else was a consequence.

Where to go next

• The PTX ISA. Long, dry, and the only place several things in Part 6 are written down.
• CUTLASS and CuTe, plus Colfax's tutorials, which are the best available writing on Hopper-era kernel structure.
• Programming Massively Parallel Processors (Hwu, Kirk, El Hajj) — the textbook this series compresses. Parts 1–5 are its early chapters with measurements attached.
• Yifan Yang's blog — the swizzle post that prompted this series, and several others in the same style.
• Simon Boehm's matmul walkthrough, with the code for every rung of Part 5's ladder.
• The Triton tutorials — read 03 (matmul) and 06 (attention) against Parts 5 and 7.
• Read a real kernel end to end. FlashAttention's Hopper source, or DeepGEMM, which is small enough to finish.
• FlashInfer's blog and vLLM's — production notes on decode attention, KV-cache layout, and the serving tricks of Part 7, written by the people maintaining them.