Kernels, tiling, coalescing and fusion

Before this lesson: lesson 1 (grid, blocks, warps, SMs, shared memory) and lesson 2 (arithmetic intensity and the roofline). Time: about 65 minutes, plus the notebook.

---

Why this matters #

• torch.compile and JAX/XLA speed up models mainly by fusing chains of elementwise ops and reductions (lesson 13). To predict what they gain, you count bytes, as in this lesson.
• FlashAttention (lesson 4) is tiling plus fusion applied to attention.
• Matrix multiply libraries (cuBLAS, CUTLASS) reach a large fraction of the peak only because of tiling at several levels.
• Layouts matter in science code too: an FFT along a non-contiguous axis, or a channels-first versus channels-last tensor, decides whether loads are coalesced.

1. Predict first #

2. The picture #

One kernel, one round trip to HBM. In eager PyTorch, each operation is its own kernel. The kernel reads its inputs from HBM, computes, and writes its output back to HBM. In a chain of four elementwise ops, the data crosses the path between HBM and the chip eight times: four reads and four writes. A fused kernel reads the input once, keeps the three intermediate results in registers, and writes only the final output: two crossings.

Every kernel also has a fixed cost. Before the GPU can run a kernel, the CPU must prepare and launch it. In eager PyTorch this costs a few microseconds per operation (Horace He measured about 280,000 PyTorch operations per second on tiny tensors, that is, about 3.6 µs each). The CPU queues kernels while the GPU runs earlier ones, so in a long chain each kernel takes about

where $o$ is the launch cost per kernel. Small tensors are overhead-bound ($o$ wins); large tensors are memory-bound (bytes/BW wins).

Interactive figure: loading…

Reuse needs fast memory. In $C = AB$, row $i$ of $A$ is needed by every output in row $i$ of $C$, and column $j$ of $B$ by every output in column $j$. A naive kernel fetches them again from HBM for every output. A tiled kernel loads a $T \times T$ tile of $A$ and a $T \times T$ tile of $B$ into shared memory once, and uses each loaded value $T$ times. The reads drop by the factor $T$.

Interactive figure: loading…

The hardware reads in chunks. HBM is not read one number at a time. The smallest unit of a read is a 32-byte sector. When the 32 threads of a warp load 32 neighbouring fp32 values (128 bytes), the hardware fetches 4 sectors and every byte is used. When the 32 threads load values that are far apart, each value costs its own sector: 32 sectors (1024 bytes) for 128 useful bytes.

Interactive figure: loading…

3. Definitions and main results #

The CUDA programming model #

A kernel is a function that runs on the GPU in many threads at once. In CUDA C++ it is marked __global__. The launch kernel<<<grid, block>>>(...) starts grid blocks of block threads each. Every thread reads its own coordinates: threadIdx (its position in the block), blockIdx (the position of its block in the grid) and blockDim (the number of threads per block). The simplest kernel, vector add, gives each thread one element:

__global__ void add(const float* x, const float* y, float* out, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // global index of this thread
    if (i < n) out[i] = x[i] + y[i];                  // guard: the last block may stick out
}
// launch: enough blocks of 256 threads to cover n elements
add<<<(n + 255) / 256, 256>>>(x, y, out, n);

Step by step, for $n = 1000$ and blocks of 256 threads:

1. The launch starts $\lceil 1000/256 \rceil = 4$ blocks, so $4 \times 256 = 1024$ threads.
2. Thread 5 of block 2 computes $i = 2 \cdot 256 + 5 = 517$ and writes out[517].
3. Threads with $i \ge 1000$ (the last 24 threads of block 3) do nothing. Without the guard if (i < n) they would write outside the arrays.

The rules of lesson 1 decide what a kernel may do:

1. Blocks are independent. They may run in any order, on any SM. There is no barrier between the blocks of one kernel; the end of the kernel is the barrier for the whole grid.
2. The threads of one block run on the same SM and share shared memory (__shared__ arrays, up to 227 KB per block on an H100). Other blocks cannot see it.
3. __syncthreads() is a barrier for the threads of one block: no thread passes it until all threads of the block have reached it. Put it between writing shared memory and reading what other threads wrote. Every thread of the block must reach it.
4. The block runs as warps of 32 consecutive threads (threadIdx.x varies fastest).

Each launch also has a fixed launch cost: the CPU time to prepare and submit the kernel (a few µs per operation in eager PyTorch), plus a small time on the GPU to start it.

Triton: one program per tile #

Triton is a Python language for GPU kernels; torch.compile also generates Triton code. You write the code of one program instance, which processes a whole block (tile) of data with vector operations. The compiler decides how threads and warps share the work, makes the loads coalesced, and places data in shared memory. The same vector add:

import triton
import triton.language as tl

@triton.jit
def add_kernel(x_ptr, y_ptr, out_ptr, n, BLOCK: tl.constexpr):
    pid = tl.program_id(axis=0)                  # which tile this program instance owns
    offsets = pid * BLOCK + tl.arange(0, BLOCK)  # BLOCK element indices at once
    mask = offsets < n                           # the last tile may stick out
    x = tl.load(x_ptr + offsets, mask=mask)
    y = tl.load(y_ptr + offsets, mask=mask)
    tl.store(out_ptr + offsets, x + y, mask=mask)

# launch: one program instance per tile of BLOCK elements
# grid = (triton.cdiv(n, 1024),)
# add_kernel[grid](x, y, out, n, BLOCK=1024)

• tl.program_id(0) plays the role of blockIdx.x. There is no threadIdx and no __syncthreads(): the compiler maps the block operations to threads and adds the barriers.
• tl.arange(0, BLOCK) is a whole vector of indices; BLOCK must be a compile-time power of 2.
• The mask switches off the positions beyond n in the last tile. It is the guard of the CUDA kernel.
• You still choose the tile shape (often by autotuning: timing a few candidates) and keep the fastest-moving index contiguous in memory.

Both kernels need a GPU. In this lesson, read them; the notebook simulates what they do.

Coalescing #

Global memory is accessed in 32-byte sectors; a 128-byte cache line has 4 sectors. The 32 addresses of one warp-wide load are combined: the load costs as many sectors as its addresses touch. A load is coalesced when consecutive threads access consecutive addresses, so that the warp uses every byte of every sector it fetches.

Example: a row-major fp32 matrix with 4096 columns. Element $(i, j)$ is at byte $4(4096\,i + j)$.

• Thread $t$ reads $(i, t)$, a piece of a row: 32 consecutive floats = 128 bytes = 4 sectors. Efficiency 100%.
• Thread $t$ reads $(t, j)$, a piece of a column: the addresses are 16 KB apart, one sector each. 32 sectors = 1024 bytes for 128 useful bytes. Efficiency 12.5%: 8× the traffic.

So the layout of a tensor and the mapping of threads to elements decide the real bytes moved. The byte counts of lesson 2 assume coalesced access.

Shared-memory banks #

Shared memory is split into 32 banks. Successive 4-byte words go to successive banks:

Each bank serves one word per cycle. If the 32 threads of a warp access 32 different banks, the access takes one pass. If $k$ threads access different words in the same bank, the hardware splits the access into $k$ passes: a $k$-way bank conflict. Threads that read the same word get it in one pass (a broadcast), with no conflict.

Take float tile[32][32], stored row by row: element tile[i][j] is word $32i + j$.

• Row read (thread $t$ reads tile[r][t]): word $32r + t$, bank $t$. All banks differ: no conflict.
• Column read (thread $t$ reads tile[t][c]): word $32t + c$, bank $c$ for every $t$. All 32 threads hit the same bank: a 32-way conflict, 32 passes.

The padding trick. Declare float tile[32][33]. Now tile[i][j] is word $33i + j$, and the column read hits bank $(33t + c) \bmod 32 = (t + c) \bmod 32$. As $t$ runs over $0, \dots, 31$, these are 32 different banks: no conflict. In words: the extra (unused) word at the end of each row shifts the next row by one bank, so a column runs diagonally across the banks. The cost is 32 unused words per tile.

Interactive figure: loading…

The tiled transpose #

A transpose writes out[j][i] = in[i][j]. A naive kernel reads rows of in (coalesced) but writes columns of out: each thread of a warp writes to its own sector. The tiled transpose moves the strided access into shared memory, where it is cheap:

#define TILE 32
__global__ void transpose(const float* in, float* out, int rows, int cols) {
    __shared__ float tile[TILE][TILE + 1];           // +1 column: no bank conflicts
    int x = blockIdx.x * TILE + threadIdx.x;         // column of `in`
    int y = blockIdx.y * TILE + threadIdx.y;         // row of `in`
    if (y < rows && x < cols)
        tile[threadIdx.y][threadIdx.x] = in[y * cols + x];   // coalesced read of a row
    __syncthreads();                                 // the whole tile is in shared memory
    x = blockIdx.y * TILE + threadIdx.x;             // column of `out` (= row of `in`)
    y = blockIdx.x * TILE + threadIdx.y;             // row of `out`
    if (y < cols && x < rows)
        out[y * rows + x] = tile[threadIdx.x][threadIdx.y];  // coalesced write of a row
}
// launch: one block of 32 x 32 threads per 32 x 32 tile
dim3 block(TILE, TILE), grid((cols + TILE - 1) / TILE, (rows + TILE - 1) / TILE);
transpose<<<grid, block>>>(in, out, rows, cols);

1. Each block reads one $32 \times 32$ tile of in, row by row (coalesced), into shared memory.
2. __syncthreads(): now the whole tile is in shared memory.
3. Each block writes the transposed tile to out, row by row (coalesced). To do this it reads the tile by columns, tile[threadIdx.x][threadIdx.y]. The padding removes the 32-way conflict of this column read.

Mark Harris measured these kernels on an NVIDIA Tesla K20c (2012), with the same tiles but blocks of $32 \times 8$ threads that copy 4 elements each. A plain copy reached 136 GB/s. The naive transpose reached 55 GB/s, the tiled one 98 GB/s, and the tiled and padded one 144 GB/s: as fast as a copy.

Fusion: counting HBM traffic #

Count element transfers per kernel: elements read plus elements written. For the unfused version, add up over all kernels. For the fused version, count only the external inputs and the final outputs; the intermediate results stay in registers. For a chain of $k$ unary elementwise ops on $N$ elements:

If each op also reads one extra input tensor of $N$ elements (like $y = x \cdot w_i$), the counts are $3kN$ and $(2 + k)N$.

Then predict the time with $t_\text{kernel} \approx \max(o, \text{bytes}/\mathrm{BW})$:

• Memory-bound (large $N$): the speed-up is the ratio of bytes. Unary chain: $2kN / 2N = k$. Chain with extra inputs: $3k/(2 + k)$, so 2× for $k = 4$.
• Overhead-bound (small $N$): the speed-up is the ratio of kernel counts, $k / 1 = k$, whatever the bytes.
• The border: one unary bf16 kernel moves $4N$ bytes, which takes as long as a 3.6 µs launch when $4N / 3.35 \times 10^{12} = 3.6 \times 10^{-6}$, that is, $N \approx 3 \times 10^6$ on an H100.

Fusion also works around a matrix multiply. A GEMM epilogue applies the bias and the activation to each output tile before it is written, which saves one write and one read of the output. And fusion is why FlashAttention exists: unfused attention writes the $s \times s$ score matrix to HBM and reads it back several times (lesson 4).

Theorem (tiled matmul traffic) #

Before the general statement, count a small case. Take $M = N = K = 4$ and tiles of $T = 2$. $C$ has $2 \times 2 = 4$ output tiles. For each output tile, the kernel steps through $K$ in $4/2 = 2$ steps, and each step loads one $2 \times 2$ tile of $A$ and one of $B$: $4 + 4 = 8$ elements. Total: $4 \text{ tiles} \times 2 \text{ steps} \times 8 = 64$ reads. The naive kernel reads a row of $A$ and a column of $B$ (4 + 4 elements) for each of the 16 outputs: $128$ reads. The ratio is $T = 2$.

Theorem. Let $C = AB$ with $A$ of size $M \times K$, $B$ of size $K \times N$, all in HBM. Count traffic in elements.

(a) The naive algorithm (each $C_{ij}$ from row $i$ of $A$ and column $j$ of $B$) reads $2MNK$ elements and writes $MN$.

(b) The tiled algorithm uses output tiles of at most $T \times T$. For each output tile it steps through the $K$ dimension in pieces of at most $T$, reads one tile of $A$ and one tile of $B$ per step, and accumulates. It writes each $C$ tile once. It reads

elements, with equality when $T$ divides $M$ and $N$, and writes $MN$.

(c) On a GPU, the $A$ tile and the $B$ tile live in shared memory, and the $C$ tile is held in registers, spread over the threads of the block. With $e$ bytes per element and $S$ bytes of shared memory per block, the tile size is limited by

Proof of (a) and (b).

1. Naive: each of the $MN$ outputs reads $K$ elements of $A$ and $K$ of $B$. That is $2MNK$ reads, and one write per output.
2. Tiled: split the rows of $C$ into $\lceil M/T \rceil$ row sets $I_i$, the columns into $\lceil N/T \rceil$ column sets $J_j$, and the $K$ dimension into $\lceil K/T \rceil$ sets $K_k$. Only the last set of each split can be smaller than $T$.
3. For the output tile $(i, j)$ and the step $k$, the kernel reads the $A$ tile $I_i \times K_k$ ($|I_i|\,|K_k|$ elements) and the $B$ tile $K_k \times J_j$ ($|K_k|\,|J_j|$ elements).
4. Reads of $A$: $\sum_i \sum_j \sum_k |I_i|\,|K_k| = \big(\sum_i |I_i|\big)\big(\sum_k |K_k|\big)\left\lceil N/T \right\rceil = MK \left\lceil N/T \right\rceil$. The summand does not depend on $j$, so the sum over $j$ is a factor $\lceil N/T \rceil$; the row sets add up to $M$ and the $K$ sets add up to $K$.
5. In the same way, the reads of $B$ are $KN \lceil M/T \rceil$.
6. If $T$ divides $M$ and $N$, the total is $MK \cdot N/T + KN \cdot M/T = 2MNK/T$: $T$ times less than naive.
7. Each $C$ tile stays on the chip as an accumulator over all steps $k$ and is written once: $MN$ writes. $\blacksquare$

Part (c) is a budget, not a proof: the two input tiles take $2T^2$ elements of shared memory. Real kernels also keep extra tiles in flight (see the box below), so they use less than this maximum.

This is close to the best possible. Hong and Kung (1981) proved that any order of the $2MNK$-FLOP matmul with a fast memory of $S$ elements needs $\Omega(MNK/\sqrt{S})$ transfers. Tiling reaches this bound up to a constant factor, because $T \approx \sqrt{S/2}$.

Intensity of the tiled matmul. Per step, a block loads $2T^2$ elements and does $2T^3$ FLOPs ($T^2$ outputs, $T$ multiply-adds each). In bf16 that is

To pass the H100 ridge point of 295 FLOP/byte with HBM traffic alone, $T$ would have to be about 590. But two bf16 tiles of $590 \times 590$ need 1.4 MB, and one block has at most 227 KB.

Real matrix multiply kernels (cuBLAS, CUTLASS) still reach a large fraction of the peak, for two reasons:

• L2 reuse across blocks. The 132 SMs of an H100 work on many output tiles at the same time, and neighbouring tiles need the same rows of $A$ or columns of $B$. When these loads hit in the 50 MB L2 cache, HBM sees the traffic of one large "super-tile". Example: 132 blocks with $128 \times 256$ output tiles, arranged as $12 \times 11$ tiles, cover a $1536 \times 2816$ region of $C$. By the matmul intensity theorem of lesson 2 (for a long $K$ dimension), its intensity is $1536 \cdot 2816 / (1536 + 2816) \approx 994$ FLOP/byte in bf16, far above the ridge. This is why the launch order of the blocks matters; the Triton matmul tutorial groups blocks for this reason.
• Register tiling. Inside a block, each warp and each thread computes a small sub-tile of $C$ in registers. This raises the intensity of the shared-memory traffic in the same way, one level down.

Answers to section 1. (1) For $N = 10^8$: up to 4×. Unfused, the chain moves $8N$ elements = 1.6 GB in bf16, 478 µs at 3.35 TB/s; fused, $2N$ = 0.4 GB, 119 µs. For $N = 10^4$: the fused kernel moves 40 KB, about 12 ns of bandwidth time, far below a launch cost of a few µs. So the four launches against one decide, and the gain is again about 4×, but for a different reason (fewer launches, not fewer bytes). (2) By the tile size, 128×: from $2MNK$ to about $2MNK/128$ reads.

What cuBLAS and CUTLASS add

• Rectangular tiles and a K step. Real kernels use output tiles such as $128 \times 128$ or $128 \times 256$ and step through $K$ in smaller pieces (for example 32 or 64), so each step loads a $128 \times 64$ tile of $A$ and a $64 \times 256$ tile of $B$. The traffic count is the same as in the theorem, with $T$ replaced by the output tile sizes.
• Double buffering (software pipelining). While the tensor cores multiply the tiles of step $k$, the next tiles are already being copied into a second buffer in shared memory. On Ampere GPUs (A100) this uses asynchronous copy instructions (cp.async); the kernel keeps 2–4 stages in flight. This hides the HBM latency of lesson 1, at the price of 2–4× the shared memory per tile.
• Hopper (H100). The Tensor Memory Accelerator (TMA) copies whole tiles between global and shared memory, started by one thread, so the other threads do not compute addresses. The wgmma instructions let a warpgroup (4 warps, 128 threads) run one large matrix multiply on the tensor cores. The operands are read from shared memory (one of them may also come from registers).
• Swizzling. Instead of padding, libraries often store tiles in a permuted ("swizzled") order, which also spreads a column over all banks but wastes no memory.

4. Worked example #

A. Fusing RMSNorm. Eager RMSNorm, written as w * (x * torch.rsqrt(x.pow(2).mean(-1, keepdim=True) + eps)) (bf16, no fp32 casts), runs on $x$ of shape $16{,}384 \times 4096$ on an H100: $N = 2^{26}$ elements, $h = 4096$ features per row. Count element transfers per kernel (reads plus writes):

Kernel	Reads	Writes
x.pow(2)	$N$	$N$
.mean(-1)	$N$	$N/h$
+ eps, rsqrt (2 kernels)	$2N/h$	$2N/h$
x * r	$N + N/h$	$N$
w * (...)	$N + h$	$N$
total	$\approx 4N$	$\approx 3N$

1. Unfused: about $7N$ elements $= 7 \times 2^{26} \times 2$ bytes $= 940$ MB, so 281 µs at 3.35 TB/s.
2. Fused: read $x$ once and write the output once: $2N = 268$ MB, so 80 µs.
3. Predicted speed-up: $7N / 2N = 3.5\times$.
4. How the fused kernel works: it needs the mean of a row before it can scale that row. So one program instance (or one warp) takes whole rows. A row of 4096 bf16 values is 8 KB; it stays in registers or shared memory. The kernel reads it from HBM once, reduces it, and scales it. This is the pattern of the Triton fused-softmax and layer-norm tutorials.
5. At decode time the same norm may run on 8 tokens ($N = 32{,}768$): 0.04 µs of bandwidth time. Now the six launches decide: about $6 \times 3.6 = 21.6$ µs unfused against 3.6 µs fused, about 6×.

B. Tiled matmul traffic. $M = N = K = 4096$ in bf16 on an H100. At the peak of 989 TFLOP/s, the $2 \cdot 4096^3 = 1.37 \times 10^{11}$ FLOPs take 139 µs.

1. Naive: $2MNK = 1.37 \times 10^{11}$ element reads = 275 GB, so 82 ms: about 600× the compute time.
2. Tiled with $T = 128$: $2MNK/128 = 1.07 \times 10^9$ reads = 2.15 GB, so 641 µs. The intensity is $T/2 = 64$ FLOP/byte: still memory-bound, if every tile came from HBM.
3. The minimum (each matrix crosses HBM once): $3 \times 4096^2 \times 2$ bytes = 101 MB, 30 µs.
4. Tuned kernels still reach a large fraction of the peak at this size, through L2 reuse and register tiling (section 3).

5. In code #

NumPy checks the sector and bank arithmetic of section 3, and a tiled matmul counts its reads. The notebook for this lesson has these cells and the exercises.

import numpy as np

lane = np.arange(32)                           # the 32 threads of one warp

# Coalescing: which 32-byte sectors does one warp-wide load touch?
cols = 4096                                    # row-major fp32 matrix, element (i, j) at byte 4*(cols*i + j)
row_load = 4 * (cols * 7 + lane)               # thread t reads (7, t): part of a row
col_load = 4 * (cols * lane + 7)               # thread t reads (t, 7): part of a column
for name, addr in [("row", row_load), ("column", col_load)]:
    sectors = np.unique(addr // 32).size
    print(f"{name:6s}: {sectors:2d} sectors, {32 * sectors:4d} bytes fetched, {128 / (32 * sectors):.1%} used")
# row   :  4 sectors,  128 bytes fetched, 100.0% used
# column: 32 sectors, 1024 bytes fetched, 12.5% used

# Bank conflicts: 32 banks of 4 bytes, word w is in bank w % 32.
def conflict_ways(words):                      # the most distinct words that share one bank
    banks = words % 32
    return max(np.unique(words[banks == k]).size for k in np.unique(banks))

for width in (32, 33):                         # float tile[32][width]; thread t reads tile[t][5]
    print(f"tile[32][{width}] column read: {conflict_ways(lane * width + 5)}-way")
# tile[32][32] column read: 32-way
# tile[32][33] column read: 1-way

# Tiled matmul: count the elements read from "HBM" and check the result.
def tiled_matmul(A, B, T):
    M, K = A.shape; N = B.shape[1]
    C, reads = np.zeros((M, N)), 0
    for i in range(0, M, T):
        for j in range(0, N, T):
            acc = np.zeros((min(T, M - i), min(T, N - j)))   # the C tile stays "on chip"
            for k in range(0, K, T):
                a, b = A[i:i+T, k:k+T], B[k:k+T, j:j+T]      # load one tile of A and one of B
                reads += a.size + b.size
                acc += a @ b
            C[i:i+T, j:j+T] = acc                             # write the C tile once
    return C, reads

rng = np.random.default_rng(0)
A, B = rng.standard_normal((64, 64)), rng.standard_normal((64, 64))
for T in (1, 4, 16):
    C, reads = tiled_matmul(A, B, T)
    print(T, np.allclose(C, A @ B), reads, 2 * 64**3 // T)
# 1 True 524288 524288
# 4 True 131072 131072
# 16 True 32768 32768

torch.compile fuses the RMSNorm of section 4. On a CPU it generates one fused C++ kernel (named cpp_fused_add_mean_mul_pow_rsqrt_0) instead of six eager kernels, and the byte count predicts the gain. This block needs PyTorch; it is not in the notebook.

import time
import torch

def rmsnorm(x, w, eps=1e-6):              # eager: six kernels (pow, mean, add, rsqrt, mul, mul)
    r = torch.rsqrt(x.pow(2).mean(-1, keepdim=True) + eps)
    return w * (x * r)

def best_time(fn, repeats=10):
    fn()                                  # warm-up (and compilation)
    times = []
    for _ in range(repeats):
        t0 = time.perf_counter(); fn(); times.append(time.perf_counter() - t0)
    return min(times)

torch.manual_seed(0)
x, w = torch.randn(8192, 4096), torch.randn(4096)   # fp32, N = 2^25 elements (128 MiB)
fused = torch.compile(rmsnorm)
print(torch.allclose(fused(x, w), rmsnorm(x, w), atol=1e-5))   # True
n_bytes = 4 * x.numel()
t_eager, t_fused = best_time(lambda: rmsnorm(x, w)), best_time(lambda: fused(x, w))
print(f"eager    {t_eager * 1e3:5.1f} ms, {7 * n_bytes / t_eager / 1e9:3.0f} GB/s for 7N transfers")
print(f"compiled {t_fused * 1e3:5.1f} ms, {2 * n_bytes / t_fused / 1e9:3.0f} GB/s for 2N transfers")
print(f"speed-up {t_eager / t_fused:.1f}x (the byte count predicts 3.5x)")

On the Apple-silicon laptop where this lesson was checked (PyTorch 2.15, CPU), eager RMSNorm took 20–26 ms and the compiled one 6–9 ms: a speed-up of 3.0–3.3×. Both ran at about the same bandwidth (30–50 GB/s), so the byte count predicts the time. The first call takes several seconds, because it compiles. Your times will differ; the ratio will be similar.

6. Where it shows up #

• torch.compile (TorchInductor) fuses chains of elementwise ops and reductions into generated Triton kernels on GPUs (C++ on CPUs); JAX/XLA fuses in the same way (lesson 13).
• FlashAttention is tiling plus fusion for attention: the score matrix never goes to HBM (lesson 4).
• Matrix multiply libraries (cuBLAS, CUTLASS) tile at three levels: block tiles in shared memory, warp tiles, and per-thread tiles in registers, with tensor-core instructions at the bottom.
• Fused optimizers. torch.optim.AdamW(..., fused=True) updates all parameters in a few kernels, instead of several elementwise kernels per parameter tensor.
• Inference engines fuse norms, rotary embeddings and activations, because decode steps work on few tokens and are often overhead-bound (lesson 12).
• ML for science. PDE residuals and physics-informed losses (finite differences, spectral penalties) are long elementwise chains on large grids, so fusion cuts their traffic by about the length of the chain. The data layout matters too: channels-first or channels-last, or an FFT along a non-contiguous axis of a weather grid, decides whether accesses are coalesced.

7. Common mistakes #

• Predicting a fusion gain from FLOPs. Memory-bound kernels are timed by bytes: count reads and writes per kernel.
• Forgetting that a reduction reads its whole input. mean(-1) reads $N$ elements even though it writes only $N/h$.
• Forgetting the launch cost. For small tensors, the number of kernels, not the bytes, sets the time.
• __syncthreads() in a branch that not all threads take. The behaviour is undefined (often a hang).
• Expecting a barrier between blocks. In a normal launch, blocks of one kernel cannot wait for each other; split the work into two kernels.
• Counting only one tile in the shared-memory budget. The $A$ tile and the $B$ tile must both fit ($2T^2 e$ bytes), and pipelined kernels need several of each.
• Assuming any read of $n$ bytes costs $n$ bytes. A strided or misaligned warp access fetches whole 32-byte sectors.
• Fusing without measuring. A fused kernel that needs too many registers loses occupancy (lesson 1), and one that recomputes expensive values can be slower.

8. Check yourself #

1. Count element transfers for out = silu(a) * b (a SwiGLU gate; a, b and out have $N$ elements each) in eager mode (two kernels) and fused. What is the predicted speed-up for a large $N$?
   Answer

   Eager: `silu` reads $a$ and writes a temporary $t$ ($2N$); `*` reads $t$ and $b$ and writes `out` ($3N$): $5N$ in total. Fused: read $a$ and $b$, write `out`: $3N$. Speed-up $5/3 \approx 1.7\times$. With more external inputs, fusion gains less.
2. A warp reads 32 consecutive bf16 values from one row of a row-major bf16 matrix with 4096 columns, and another warp 32 values from one column. How many sectors does each load touch, and what fraction of the fetched bytes is used?
   Answer

   Row: $32 \times 2 = 64$ bytes = 2 sectors, 100% used. Column: 32 sectors = 1024 bytes for 64 useful bytes, 6.25% used (16× the traffic).
3. Count the reads of the tiled matmul for $M = K = N = 1000$ and $T = 128$. Compare with $2MNK/T$ and with the naive algorithm.
   Answer

   $\lceil 1000/128 \rceil = 8$, so reads $= 10^6 \cdot 8 + 10^6 \cdot 8 = 16{,}000{,}000$. The estimate $2MNK/T$ gives 15,625,000. Each element is read $\lceil 1000/128 \rceil = 8$ times, not $7.8$ times, so the exact count is 2.4% higher. Naive: $2 \times 10^9$ reads, 125× more.
4. What is the largest square tile for which two fp32 tiles fit in 227 KB (232,448 bytes) of shared memory, and what is the intensity of one step with it?
   Answer

   $2 \cdot T^2 \cdot 4 \le 232{,}448$, so $T^2 \le 29{,}056$ and $T = 170$ ($170^2 = 28{,}900$). One step loads $2T^2$ fp32 elements ($8T^2$ bytes) and does $2T^3$ FLOPs: $T/4 = 42.5$ FLOP/byte.
5. A warp reads the diagonal tile[t][t] of float tile[32][32]. Is there a bank conflict? And for tile[t][0]?
   Answer

   Diagonal: word $33t$, bank $33t \bmod 32 = t$: all different, no conflict. `tile[t][0]`: word $32t$, all in bank 0: a 32-way conflict.

Recall cardsAnswer in your head first, then reveal. Add the cards to your review deck to see them again just before you would forget.

1. What is the smallest unit in which a GPU reads global memory (HBM)?

   Show answer

2. When is a warp's load coalesced?

   Show answer

3. Compute the efficiency of a warp load where each thread reads one fp32 value from a different row of a row-major matrix with 4096 columns.

   Show answer

4. How many banks does shared memory have, and which bank holds word $w$?

   Show answer

5. Why does padding float tile[32][32] to tile[32][33] remove the conflict of a column read?

   Show answer

6. What does a kernel launch cost in eager PyTorch, and when does it set the time?

   Show answer

7. How many element transfers do an unfused and a fused chain of $k$ unary elementwise ops on $N$ elements need?

   Show answer

8. How many elements does a tiled matmul read from HBM, compared with the naive one?

   Show answer

9. What limits the tile size $T$ of a matmul kernel?

   Show answer

Add these cards to my reviewOpen the review →

9. Interview angle #

• "Write a CUDA kernel for vector add, then for a transpose. What makes the transpose fast?" — The global index and the bounds guard; a shared-memory tile so that both the read and the write are coalesced; padding to 33 columns against bank conflicts.
• "What is kernel fusion, and when does it help?" — It removes HBM round trips of intermediate results and kernel launches. It helps memory-bound chains (by the ratio of bytes) and small tensors (by the ratio of launches), not compute-bound matmuls.
• "How does tiling reduce the memory traffic of a matmul?" — From $2MNK$ to about $2MNK/T$ reads; the intensity per step is $T/2$ FLOP/byte in bf16; this is optimal up to a constant (Hong–Kung).
• "What is a bank conflict?" — Several threads of a warp access different words in the same bank, so the access is split into passes. Fix it with padding or a swizzled layout.
• "Why is a strided access slow even if the total data is small?" — The hardware fetches whole 32-byte sectors; a stride wastes the unused part of each sector.

10. Practice #

The notebook for this lesson runs in Colab, with NumPy only. It has five exercises, each with a check cell and a folded solution:

1. sectors_touched(addresses): the sectors and the efficiency of one warp-wide load.
2. bank_conflict_degree(tile_width, access) for row, column and diagonal reads.
3. A naive and a tiled matmul on a CountingMemory that counts reads from "HBM", checked against A @ B, $2MNK$ and $2MNK/T$.
4. The fusion time model: the time of an unfused and a fused chain, and the size $N$ where a kernel changes from overhead-bound to memory-bound.
5. A tiled transpose in NumPy blocks that counts contiguous and strided accesses.

11. Watch and read #

Lecture 5: Going Further with CUDA for Python ProgrammersJeremy Howard · GPU MODE lecture 5 · 78 min

Shared memory and a tiled matrix multiply, built step by step from Python.

Lecture 14: Practitioners Guide to TritonUmer Adil · GPU MODE lecture 14 · 82 min

Triton's programming model (one program per tile, masks), vector add and a tiled matmul.

• GPU MODE, lecture notes and code. Several lectures follow the book by Hwu, Kirk and El Hajj, Programming Massively Parallel Processors (4th edition, 2022), whose chapters on memory coalescing and shared-memory tiling match this lesson.
• Triton, tutorials: "Vector Addition", "Fused Softmax" (its traffic count of the naive softmax) and "Matrix Multiplication" (the section on L2 cache optimizations).
• Mark Harris, An Efficient Matrix Transpose in CUDA C/C++ (NVIDIA technical blog): the three transpose kernels of section 3, with measurements.
• NVIDIA, CUDA C++ Best Practices Guide: the sections "Coalesced Access to Global Memory" and "Shared Memory and Memory Banks".
• Simon Boehm, How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: from a naive kernel to block tiling, register tiling and vectorized loads, with the speed of each step.
• Horace He, Making Deep Learning Go Brrrr From First Principles: overhead, bandwidth and compute, and why fusion helps.
• Hong and Kung (1981), "I/O complexity: The red-blue pebble game", STOC: the lower bound on matmul traffic.