GPU Tensor Cores — Hardware Matrix Math for Deep Learning

The throughput gap is the whole story

A conventional FP32 lane performs one multiply-add per cycle. An MMA instruction performs an entire small matrix product — hundreds of multiply-adds — in one issue, because the multipliers and adders are wired as a fixed array rather than as independent lanes. The result is that the matrix pipe's peak rate sits roughly an order of magnitude above the same SM's vector FP32 rate, and each step down the operand ladder (TF32, bf16/fp16, fp8, int8/int4) roughly doubles it again.

Memory did not follow. HBM bandwidth has improved substantially generation over generation, but nowhere near the multiple that matrix throughput has. The two curves diverging is the single most important fact about modern GPU performance work.

The roofline balance point — peak FLOP/s divided by peak bytes/s, the arithmetic intensity at which a kernel stops being bandwidth-limited — used to sit around ten FLOPs per byte on pre-tensor-core parts. Measured against the matrix pipe it now sits in the hundreds. Almost no kernel written naively clears that bar. Elementwise ops do not, normalisation does not, attention as commonly written does not, and neither does a matmul whose operands stream in from HBM once per use.

So the first number to compute for any kernel is not its FLOP count. It is FLOPs divided by bytes of HBM traffic, compared against that balance point. The architecture overview sets out the roofline in general; the tensor-core-specific version of it is simply that the ridge moved so far right that being compute-bound is now the exception you have to engineer for.

Only a GEMM has enough reuse to feed one

Arithmetic intensity is a property of the shape, not of the hardware. Work it out for a dense matmul C = A × B with dimensions M×N×K. The FLOP count is 2 × M × N × K. If every operand element were touched exactly once, the traffic would be M×K + K×N + M×N elements. For a square problem that is 2N³ FLOPs against 3N² elements: intensity grows linearly with N. A large GEMM is essentially the only common shape whose intensity can be pushed above the balance point at all.

Everything else falls short structurally. An elementwise op reads one element, does a handful of FLOPs and writes one element — intensity below one, permanently memory-bound, and no tensor core anywhere will change that. The only lever is to move fewer bytes, which is what kernel fusion is for.

Matrix-vector is the case that surprises people. With M = 1, the FLOP count is 2NK and the weight matrix alone is NK elements, so intensity is about two FLOPs per element regardless of how large the matrix is. Autoregressive decode at batch size one is therefore a bandwidth exercise: the GPU is reading weights, not multiplying them. That is why batching, KV-cache layout and weight quantization move the needle there while matmul tuning does not.

Attention sits in between: QK⊤ and the value product are real GEMMs, but the head dimension is small and materialising the S×S score matrix costs quadratic traffic — the observation FlashAttention is built on. Before optimising anything, classify the shape.

The tile hierarchy is where reuse is manufactured

That intensity calculation assumed each operand element is fetched from HBM exactly once. Nothing gives you that for free. The reuse has to be constructed by blocking, and a GEMM kernel is essentially three nested tiles: the CTA tile (the output block one thread block computes), the warp tile it is subdivided into, and the instruction tile, which is fixed by the hardware and the operand dtype and is not yours to choose.

The reuse arithmetic is worth doing once by hand. A thread block computing a BM×BN output tile, stepping K in chunks of BK, stages BM×BK + BK×BN operand elements per step and performs BM×BN×BK multiply-accumulates on them. The ratio is (BM×BN)/(BM+BN): 64 MACs per staged element at 128×128, 42.7 at 128×64, 32 at 64×64. Bigger tiles are strictly better for bandwidth.

They stop growing because of capacity, not because of diminishing returns. The accumulator tile lives in registers for the whole K loop: a 128×128 fp32 accumulator is 16,384 registers per block, a large fraction of an SM's register file. Shared memory must simultaneously hold several pipeline stages of staged operands. Both budgets are fixed per SM and both bound how many blocks can be resident, so tile size trades directly against occupancy and against pipeline depth.

// An illustrative GEMM tile configuration (CUTLASS-style naming)
ThreadblockShape = 128 x 128 x 32   // BM x BN x BK  -- your choice
WarpShape        =  64 x  64 x 32   // 4 warps, 2x2 over the CTA tile
InstructionShape =  16 x   8 x 16   // fixed by hardware + operand dtype
Stages           = 4                // shared-memory pipeline depth

// per K-step, per thread block:
//   operand elements staged : 128*32 + 32*128   = 8,192
//   multiply-accumulates    : 128*128*32        = 524,288
//   reuse                   : 64 MACs per staged element
//   smem for operands       : 4 stages * 8,192 * 2 B = 64 KB
//   fp32 accumulators       : 128*128 = 16,384 registers per block

Read the comment block as the actual design constraint: pick BM/BN for reuse, then discover how many stages the remaining shared memory affords.

Pipeline depth is a capacity calculation

Write the K loop naively — copy a tile into shared memory, synchronise the block, run the MMAs, synchronise, repeat — and the matrix pipe is idle for the entire copy. An HBM load takes several hundred cycles; the MMAs for one K-step take far fewer. A monolithic loop therefore leaves the most expensive unit on the chip idle most of the time, and it will do so no matter how well the tile is chosen.

The fix is to overlap: issue the copy for stage k+1 while the pipe computes stage k, with shared memory holding several stages in a circular buffer. The copy mechanism is a separate topic — asynchronous copy lets a load land in shared memory without staging through registers or stalling the issuing warp, and on newer parts the Tensor Memory Accelerator moves whole tiles from a descriptor with one thread issuing.

What belongs here is the sizing rule. Stages needed is roughly ceil(load latency / compute time per stage). Too shallow and the pipe still stalls; deeper always hides more latency, but every stage costs a full operand buffer out of an SM's fixed shared memory, which is the same budget the tile size is drawing on. Depth and tile size are two claims on one pool, and that tension is most of what a GEMM autotuner is searching.

The other structural answer is warp specialisation: dedicate some warps to producing (copies only) and the rest to consuming (MMAs only), coupled by barriers rather than by a block-wide __syncthreads(). It maps cleanly onto hardware where the copy path and the matrix pipe are genuinely separate issue resources, and it removes the whole-block synchronisation that a monolithic loop forces on every iteration.

Wave quantization - a slightly larger shape can cost a whole extra wave

A GEMM decomposes into ceil(M/BM) × ceil(N/BN) thread blocks, and those blocks are scheduled onto SMs in waves. If the device can host S blocks concurrently and the grid has G of them, the kernel takes ceil(G/S) waves, and the last wave is almost always partial. Utilisation is G / (S × ceil(G/S)), and when G is just over a multiple of S that fraction is brutal.

Concretely, and taking a device that happens to host 64 such blocks at once: with 128×128 tiles, a 1024×1024 output is 8×8 = 64 blocks and fits in exactly one wave. Grow M to 1152 and it becomes 9×8 = 72 blocks. Those eight extra blocks do not cost one-eighth more time — they cost an entire additional wave, during which most SMs sit idle. The wall clock for that shape can be indistinguishable from a considerably larger one that happens to fill the second wave properly. Substitute your own device's block capacity for the 64; the arithmetic is what matters, not the constant. Benchmark curves that look like staircases rather than lines are showing you wave quantization, not measurement noise.

Tile quantization is the finer-grained sibling. M = 129 against BM = 128 produces two row-tiles, the second of which is one row wide and 1/128 utilised; the hardware still does the full tile's work. Padding M to 256 costs nothing extra in time and is simply more honest about what you are paying for.

The practical rules follow directly. Pad batch and sequence dimensions to multiples of the tile rather than to whatever the data happened to be. Prefer shapes whose block count lands near a multiple of the machine's block capacity. Accept that for small problems a heuristic will pick a smaller tile — more blocks quantize more finely, at the cost of reuse. And when M and N are both small but K is enormous there are not enough output tiles to fill the machine at all; that is what split-K is for, partitioning the reduction across blocks and paying for a separate reduction over the partial sums in exchange for parallelism.

The epilogue, and where the accumulators go

When the K loop ends, the output tile is sitting in registers as accumulators. What happens next decides a surprising share of end-to-end time, and it is the part of a GEMM that gets designed last and costs most.

Store the tile and let a separate kernel add the bias, and another apply the activation, and the output makes a full HBM round trip per elementwise op. By the intensity argument above, each of those trips runs at intensity near one — pure bandwidth. A compute-bound body followed by three memory-bound tails can easily spend more time in the tails. This is the general case for fusion, but the GEMM instance of it is special because the data is already in registers and has never touched memory.

Hence every serious GEMM implementation carries an epilogue stage that runs on the accumulators before the single store: bias add, activation, residual add, per-channel dequantisation scale, and the cast to the output dtype. It is why cuBLASLt exposes epilogue enumerations rather than expecting you to chain kernels, and why CUTLASS makes the epilogue a template parameter you can supply.

Two things go wrong here. First, ordering: the epilogue is where the wide accumulator is narrowed, so casting before applying a scale throws away exactly the protection the wide accumulator existed to provide — see the accumulator argument in tensor core architecture. Second, split-K: partial sums from different blocks must be reduced before an epilogue can be applied at all, so a split-K kernel usually defers its epilogue into a second reduction kernel, which quietly reintroduces the round trip you were avoiding.

Which rung of the software stack to stand on

There is a ladder between "call torch.matmul" and "emit mma.sync", and most performance work fails by starting too far down it.

Descending the ladder costs engineering time and portability in the same motion: a tile configuration hand-tuned for one architecture becomes a liability on the next, because the instruction shapes, shared-memory budget and copy machinery all shift. For a worked example of what the bottom rung buys and costs, the Marlin INT4 kernel is a hand-built GEMM whose entire design — layout, register-level unpacking, double buffering — exists to keep the matrix pipe busy while dequantizing on the fly. The honest trigger for descending is a measurement, not a hunch. Compute achieved throughput for your shape as 2×M×N×K divided by measured kernel time and compare it against the device's matrix-pipe rate for that dtype. At a large fraction of peak you are done. At a small fraction, find out why before writing a kernel — wave quantization, an ineligible shape, or a memory-bound epilogue are all far cheaper to fix. The cuBLASLt heuristic API returning several ranked algorithms you can benchmark is frequently the entire win, at a fraction of the cost of a CUTLASS port.

Fed or starved - reading the verdict off a profiler

Two numbers, taken together and in this order, classify almost any matmul-heavy kernel: tensor pipe utilisation and DRAM throughput as a fraction of peak. Either one alone is misleading; the pair is diagnostic.

High pipe, low DRAM — compute-bound and near the ceiling. The only remaining levers are fewer FLOPs or a narrower operand dtype. Tiling work here is finished.

Low pipe, high DRAM — memory-bound. The pipe is starved because bytes cannot arrive fast enough. Larger tiles for more reuse, fused epilogues, quantized weights, or restructuring so tiles stay on-chip. Tuning the MMA path itself is wasted effort in this quadrant, and it is where most wasted optimisation weeks are spent.

Low pipe, low DRAM — the interesting one, and the most misread. Nothing is saturated, so the kernel is latency-bound or parallelism-starved: pipeline depth too shallow to cover the load latency, a grid too small to fill the machine, wave quantization leaving SMs idle in the tail, or a block-wide synchronisation serialising the loop. Check the grid size against block capacity before anything else.

High pipe, high DRAM is rare and means the kernel is balanced against the device.

Compute achieved FLOP/s independently rather than trusting a single counter, and plot it on a roofline if the tool offers one. Note what the evidence does not tell you: a kernel symbol naming a bf16 tensor path tells you what was selected, not what fraction of peak was achieved. See GPU profiling for capture mechanics, replay behaviour and how much of what you are reading is distortion from the tool itself.

What silently drops you off the matrix pipe

The failure mode that costs the most time produces no error at all. The library quietly selects a non-tensor kernel, the numbers come out correct, and the kernel is several times slower with nothing in the log to explain it. The usual disqualifiers are fp32 inputs with TF32 not enabled, misaligned pointers or leading dimensions, a K that is not a multiple of the instruction's K, an accidental upcast on one operand, and layout or transpose combinations with no tuned kernel behind them. Tensor core architecture has the full table and the reason each rule exists.

What matters for the data path is that eligibility and efficiency are separate questions and are usually confused. A shape can be perfectly eligible for the matrix pipe and still run at a small fraction of peak because it quantizes badly across waves, or because its epilogue costs a round trip, or because the pipeline is two stages deep when it needed four. Conversely a shape can be beautifully tiled and never touch the pipe at all because someone left allow_tf32 off. Check eligibility first — it is a yes/no question and cheap to answer — then, and only then, start reasoning about whether the pipe is being fed.