How a Matmul Kernel Gets Fast

Ten kernels, from 1.3% of the hardware to 93.7%, with every byte and every multiply accounted for.

August 2026 · a companion to Simon Boehm's CUDA matmul worklog · sequel to The CUDA Programming Model

This post is the second half of a pair. The first one took a photograph of a dog and followed it down into the silicon: grids, blocks, warps, schedulers, the clock, and the 128-byte trays the memory system serves. It never optimized anything. This one does nothing but optimize. We take a single problem, multiplying two 4092 × 4092 matrices, and rewrite the kernel ten times, from a version that reaches 1.3% of the machine to one that reaches 93.7%.

The route is Simon Boehm's, from his excellent worklog on optimizing a CUDA matmul kernel, and all the measurements quoted here are his, taken on an NVIDIA A6000.1 1A6000: 84 SMs, 768 GB/s of memory bandwidth, about 30 TFLOPs/s of fp32 arithmetic, 48 KB of shared memory per block. Boehm's code is on GitHub. What this post adds is the thing I wanted when I first read it: a picture for every step, and an explicit answer, in bytes and instructions, to the question why did that make it faster? Each kernel below gets the same treatment. What changed, what it costs, and the arithmetic that turns the change into a number of FLOPs.

1The problem, and what cuBLAS is

The operation is called SGEMM, and the name is an acronym stack worth unpacking because it explains the other name too. GEMM is GEneral Matrix Multiply. The S in front means single precision, that is, ordinary 32-bit floats. GEMM belongs to a catalogue of routines called BLAS, Basic Linear Algebra Subprograms, a standard list of operations agreed on in the 1970s that every numerical library implements: dot products, matrix-vector products, matrix-matrix products. cuBLAS is simply NVIDIA's own implementation of that catalogue for their GPUs. It is closed source, it ships with CUDA, and when PyTorch multiplies two matrices, this is very often the code that actually runs.

So cuBLAS is not an algorithm or a technique. It is the professionally tuned competitor: the number to beat, or rather to approach. In the scoreboard below, cuBLAS is 100% by definition, and every other row is a kernel you could write yourself in an afternoon or a week.

The full operation GEMM defines is C = αAB + βC, where α and β are scalars.2 2The α and β exist so the routine can also scale and accumulate into an existing C, which is what makes it a useful building block. They add two multiplies per output element and change none of the analysis; ignore them. We will treat it as plain C = AB. Here is what that means for one output element, on small matrices where you can check the arithmetic by eye.

Scaling up to the real problem: with 4092 × 4092 matrices, C has 16.7 million elements, each needing 4092 multiplies and 4092 adds. That is 2 × 4092³ ≈ 137 billion floating-point operations for one matmul. The machine we are running on can do about 30 billion of them per millisecond. The question of this entire post is why the first kernel takes 445 milliseconds to do something the hardware could do in 4.6.

2Two clocks: matmul should be compute-bound

Before optimizing anything it is worth knowing what the ceiling is and which resource we expect to hit it. There are only two candidates: the arithmetic units, and the memory system. So compute both times.

The arithmetic clock. 137 GFLOP at 30 TFLOP/s takes 4.6 milliseconds. That is the floor. No arrangement of the code beats it.

The memory clock. How many bytes must move, at minimum? A, B and C are each 4092² × 4 bytes ≈ 67 MB. You must read all of A, all of B, and read and write C: about 268 MB in total. At 768 GB/s that takes 0.35 milliseconds.

Compare them and the conclusion is emphatic. The arithmetic takes thirteen times longer than the unavoidable memory traffic. Matmul is inherently a compute-bound problem, and a perfectly written kernel should be sitting there waiting on the FP32 units with the memory system idling. That is the good news. The bad news is the third bar.

3Kernel 1: the naive version, 1.3%

The obvious kernel gives one thread to each element of C, and that thread walks the dot product in a loop.

// one thread per output element
const uint x = blockIdx.x * blockDim.x + threadIdx.x;
const uint y = blockIdx.y * blockDim.y + threadIdx.y;

float tmp = 0.0;
for (int k = 0; k < K; ++k)
    tmp += A[x * K + k] * B[k * N + y];   // 2 loads, 2 flops
C[x * N + y] = tmp;

It is correct, it is readable, and it reaches 309 GFLOPs, which is 1.3% of what cuBLAS gets and about 1% of the hardware. To see why, stop thinking about the thread and start thinking about the warp: 32 consecutive threads that issue their loads together, as one request.

In this code threadIdx.x varies fastest within a warp, and it feeds x, which is the row index. So the 32 threads of a warp own 32 consecutive rows of C, all in the same column. Now ask what they read on the very first pass of the loop, when k = 0.

Count the bill for that one k-step, for the whole warp. A costs 32 sectors × 32 bytes = 1024 bytes. B is free, since all 32 threads want the identical element and the hardware broadcasts it from one 32-byte sector. Total: 1056 bytes moved to perform 64 FLOPs (32 threads, one multiply and one add each). Boehm's profiler agrees with the diagnosis: the kernel achieves 15 GB/s of useful memory throughput on a GPU that can do 768.

4Kernel 2: memory coalescing, 8.5%

The fix does not touch the arithmetic, the loop, the block size, or the number of threads. It swaps which index each thread takes.

// same threads, same work, different assignment
// row: threadIdx.x DIVIDED by the block width
const uint x = blockIdx.x * BLOCKSIZE + (threadIdx.x / BLOCKSIZE);
// col: threadIdx.x MODULO the block width  ← the one that matters
const uint y = blockIdx.y * BLOCKSIZE + (threadIdx.x % BLOCKSIZE);

Now threadIdx.x feeds y, the column index, so a warp's 32 threads own 32 consecutive columns of one row of C. Re-run the same question. At k = 0, all 32 threads want the same element of A, which is a free broadcast. And they want 32 consecutive elements of a row of B, which in row-major storage are 32 consecutive floats: one unbroken run of 128 bytes.

This is what coalescing means, and it is worth stating precisely because the word gets used loosely. The hardware does not merge requests because it is clever; it merges them because it was always going to fetch whole 32-byte sectors, and the only question is whether the 32 threads' addresses fall inside a few sectors or thirty-two of them. Coalescing is not an optimization the GPU performs. It is a property your access pattern either has or lacks.

So here, concretely, is how a memory change becomes a FLOPs number. The kernel is limited by the rate the memory system can feed it. Feeding 64 FLOPs cost 1056 bytes and now costs 160, a factor of 6.6. If the kernel is purely memory-limited, its speed should rise by that same factor. Measured: 309 → 1986 GFLOPs, a factor of 6.4. Boehm's profiler shows the matching cause: useful throughput climbs from 15 GB/s to 110 GB/s. Two index lines, a 6× speedup, and the prediction lands within 3% of the measurement.

5The memory hierarchy, and what shared memory is

Kernel 2 fixed how efficiently we ask. The next kernels attack how often we ask, and for that we need one piece of hardware we have not used yet.

An SM has three places to keep numbers, and they differ by orders of magnitude in both speed and size. Registers are the fastest: each thread's private handful, read in the same cycle they are used. Shared memory is a scratchpad of about 48 KB per block, physically inside the SM, visible to every thread in the block, roughly 20–30× faster to reach than DRAM. And global memory is the HBM stacks off-chip: gigabytes, but hundreds of cycles away.

6Kernel 3: caching a tile, 12.8%

Kernel 2 made every request efficient. It did nothing about the number of requests, and the number is enormous. To see how enormous, keep following one specific block: the one computing the top-left 32 × 32 tile of C, with 1024 threads, one per output.

First, count what that block truly needs. Its 32 rows of output depend on rows 0–31 of A, all 4092 columns of them. Its 32 columns of output depend on columns 0–31 of B, all 4092 rows. So the block needs 2 × 32 × 4092 = 261,888 distinct numbers to do its job.

Now count what it actually fetches. Each of the 1024 threads runs the k-loop 4092 times and loads two values per iteration. That is 1024 × 4092 × 2 = 8,380,416 loads. Divide: every number the block needs is being dragged out of DRAM 32 times over.

Thirty-two is not a coincidence, and seeing where it comes from is the whole insight. Take one element of A, say A[0][7]. Which threads want it? Every thread computing an output in row 0 of the C tile, because they all walk row 0 of A. There are 32 such threads, one per column of the tile. In kernel 2 each of them goes and gets it separately.

So the fix writes itself: fetch each value once, put it somewhere every thread in the block can reach, and let all 32 of them read it from there. Registers cannot do this, since a register is private to one thread. This is exactly what shared memory is for.

Why we cannot simply preload everything. The obvious version, load the whole row strip of A and column strip of B and then compute, does not fit. The A strip alone is 32 × 4092 × 4 bytes = 512 KB, and B needs another 512 KB. A block gets 48 KB of shared memory. We are over budget by a factor of twenty-one.

The way out is to notice that the k-loop does not need all of k at once. It accumulates. So chop the k direction into chunks, and hold only one chunk at a time.

Who loads what. Here the numbers line up neatly. A chunk of A is 32 × 32 = 1024 elements, and the block has exactly 1024 threads. So each thread copies precisely one element of A and one of B from global memory into shared memory. No loops, no leftovers, and the loads are coalesced because neighbouring threads copy neighbouring elements.

What a barrier is

Here is the thing shared memory forces on us that registers never did: the threads have to wait for each other. To see why, remember that thread (5, 9) is about to read row 5 of As — thirty-two values. Did thread (5, 9) write those values? No. It wrote exactly one element. The other thirty-one were written by thirty-one other threads.

And threads do not run in lockstep. The block's 1024 threads are 32 warps, and the SM's schedulers run whichever warps are ready, in whatever order suits them. Some warps are hundreds of cycles ahead of others. So it is entirely possible for thread (5, 9) to finish its own copy, race onward, and read As[5][17] before the thread responsible for As[5][17] has written it. It would read whatever junk was left in that shared-memory location from the previous chunk, and the answer would be wrong — sometimes, unpredictably, depending on scheduling.

__syncthreads() is the fix, and it does exactly one thing: every thread in the block stops on that line until all 1024 of them have reached it. It is a meeting point. The fast threads sit and wait; when the last straggler arrives, all of them are released together. Because every thread's copy statement is before the barrier and every thread's read is after it, the barrier guarantees the whole tile is complete before anybody looks at it.

__shared__ float As[32*32], Bs[32*32];   // 8 KB, lives in the SM
float tmp = 0.0;                       // running total, in a register

for (int chunk = 0; chunk < K; chunk += 32) {

    // 1) every thread copies exactly one element of each
    As[threadRow*32 + threadCol] = A[...];
    Bs[threadRow*32 + threadCol] = B[...];

    __syncthreads();   // 2) nobody reads until all have written

    // 3) 32 multiply-adds, every operand from shared memory
    for (int d = 0; d < 32; ++d)
        tmp += As[threadRow*32 + d] * Bs[d*32 + threadCol];

    __syncthreads();   // 4) nobody overwrites until all have read
}
C[...] = tmp;           // after 128 chunks, tmp is the answer

One chunk, step by step

Let us walk a single pass of that loop with nothing left out. Follow thread (5, 9), which owns the output at row 5, column 9 of the block's tile of C.

Step 1 — every thread copies one element. Thread (5, 9) issues two loads from global memory: one element of A, one element of B. These are the slow ones, hundreds of cycles each. It writes what comes back into As[5][9] and Bs[5][9]. Every other thread does the same for its own position, so between them the 1024 threads fill both 32 × 32 tiles completely.

Step 2 — the barrier. All 1024 threads meet at __syncthreads(). When the last one arrives, As and Bs are complete and correct.

Step 3 — 32 multiply-adds. Thread (5, 9) now needs the dot product of row 5 of As with column 9 of Bs. It loops d from 0 to 31:

d = 0:   tmp += As[5][0]  * Bs[0][9]
d = 1:   tmp += As[5][1]  * Bs[1][9]
d = 2:   tmp += As[5][2]  * Bs[2][9]
  …
d = 31:  tmp += As[5][31] * Bs[31][9]

Thirty-two reads of As, thirty-two reads of Bs, thirty-two multiply-adds — and every single one of those 64 reads hits shared memory, which is on the SM, not DRAM. That is the whole point of step 1 having happened.

Step 4 — the second barrier, so that no thread starts overwriting As and Bs while a slower neighbour is still on step 3.

Step 5 — go round again. The block copies in the next 32 columns of A and next 32 rows of B, and thread (5, 9) repeats steps 2 to 4 on the new tile, adding into the same register. Its tmp is never reset and never leaves the thread.

This is the part you put your finger on: row 5 does get used again and again. Thread (5, 9) needs all 4092 values of row 5 of A, and it receives them in 32-wide bites, one per chunk. The window slides right 128 times until the row is exhausted.

Here is that accumulation happening on real numbers — one actual thread, one actual dot product of length 4092, built up 32 terms at a time.

And here is why the sharing matters so much within each of those windows. Row 5 of As is not read only by thread (5, 9). It is read by thread (5, 0), thread (5, 1), and every thread along row 5 of the tile — 32 of them, exactly as in kernel 2. The sharing never went away. What changed is where they read it from.

How that turns into FLOPs

Now the accounting, per thread, for the same 32 steps of k, so the two kernels are directly comparable.

In kernel 2, doing 32 multiply-adds required 64 trips to DRAM, and the thread spent nearly all its time waiting for them. In kernel 3 the same 32 multiply-adds require 2 trips to DRAM plus 64 reads of shared memory, which is on-chip and roughly twenty to thirty times quicker to reach, and which does not consume DRAM bandwidth at all.

Since this kernel is limited by how fast operands can be supplied, cutting the DRAM traffic by 32× is what raises the FLOPs number. Across the whole block that is 8,380,416 global loads becoming 262,144 — and 262,144 is essentially the 261,888 distinct values the block genuinely needed, which is to say each value is now fetched once.

What it bought, in numbers. Per chunk the block now issues 2048 global loads and performs 65,536 shared-memory reads against them. Each value fetched from DRAM is therefore read 32 times before being discarded, which is precisely the redundancy we started with, now paid once instead of 32 times.

Which makes the result faintly disappointing. 1986 → 2980 GFLOPs: a 32× cut in global memory traffic bought 1.5× in speed. Understanding why is the subject of the next section, and it is the hinge the whole rest of the post turns on.

7Why kernel 3 disappointed: the bottleneck moved

A 32× cut in global memory traffic bought 1.5× in speed. That gap is not a disappointment to be shrugged at; it is information. It tells us that global memory was no longer the thing we were waiting on. Something else had become the limit while we were busy fixing DRAM.

To find it, use the only rule that ever matters here. A kernel's running time is set by whichever resource runs out first, so ask, for each resource in turn, how much of it one multiply-add costs. There are three candidates: DRAM bandwidth, shared-memory service, and instruction issue.

Count the instructions. Here is kernel 3's inner loop as the machine sees it. One line of C++ becomes three machine instructions.

tmp += As[threadRow*32 + d] * Bs[d*32 + threadCol];

// LDS   — load a word from shared memory  (As)
// LDS   — load a word from shared memory  (Bs)
// FFMA  — the multiply-add itself

Now the hardware side, which we already know from part one. An SM has four warp schedulers, and each can issue one instruction per clock. That is the total budget: four instructions per cycle, for everything. If two out of every three instructions a warp issues are shared-memory loads, then at most a third of the SM's issue capacity is going into arithmetic — and the 128 FP32 lanes, which could accept a multiply-add every cycle, sit idle waiting to be given work.

The second limit points the same way. Shared memory is built from 32 banks, each able to serve one 4-byte word per cycle: about 32 words per cycle for the whole SM. At two words per multiply-add, that supply can sustain roughly 16 multiply-adds per cycle. The lanes can do 128. So even ignoring instruction issue entirely, operand supply alone caps this kernel near an eighth of the machine's arithmetic peak. Kernel 3 measures 2,980 GFLOPs against a 30,000 peak, which is a tenth — right up against that wall.3 3This is a deliberately simplified model: it ignores broadcasts, bank conflicts and latency hiding. It is not meant to predict performance, only to show which resource the kernel is pressed against.

Boehm's profiler confirms the diagnosis directly. The dominant reason warps are stalled is MIO throttle — warps queueing for the shared-memory pipeline — not anything to do with DRAM, and not the arithmetic units being busy.

8Kernel 4: 1D blocktiling, 36.5%

We now have a precise target: get below two shared-memory loads per multiply-add. So ask the question that solves it. When can one load serve more than one multiply-add?

Only when two multiply-adds want the same number. So look at the definition of the output and find where operands repeat:

C[i][j]  =  sum over k of  A[i][k] * B[k][j]

Take two outputs in the same column, C[i][j] and C[i+1][j]. At step k the first needs A[i][k]·B[k][j] and the second needs A[i+1][k]·B[k][j]. Two different values of A, but the same value of B. Take two outputs in the same row instead, and it is the other way round: the same A value, two different B values.

So give one thread a column of eight outputs. Thread now owns C[i..i+7][j], and at each step of the inner loop it needs eight different values of A and one value of B, which it loads once into a register and reuses eight times.

float threadResults[8] = {0.0};       // 8 running totals, in registers

for (uint d = 0; d < BK; ++d) {
    float Btmp = Bs[d * BN + threadCol];   // ONE load…
    for (uint i = 0; i < 8; ++i)          // …reused 8 times
        threadResults[i] +=
            As[(threadRow * 8 + i) * BK + d] * Btmp;
}

The count. Eight loads from As, one from Bs, nine in total, producing eight multiply-adds: 1.125 loads per FMA, down from 2.000. The instruction mix flips from two-thirds loads to a bit over half, and the lanes start getting fed.

The block has to be reshaped to match. If every thread now owns eight outputs, a block of the old 1024 threads would cover 8192 outputs, and the shared-memory tile it needs would not fit. Boehm's kernel 4 uses a 64 × 64 tile of C with BK = 8, which works out neatly:

• 64 × 64 = 4096 outputs per block, at 8 per thread, gives 512 threads.
• The As tile is 64 × 8 = 512 elements and the Bs tile is 8 × 64 = 512 — so once again each thread copies exactly one element of each, the same cooperative load as before.
• Shared memory: (64 + 64) × 8 × 4 bytes = 4 KB, comfortably inside the 48 KB budget.

And because the block tile grew from 32 × 32 to 64 × 64, the global traffic improves as well, independently of everything above: reads per output fall from K/16 to K/32. A bigger tile means more reuse of each value fetched from DRAM.

Measured: 2,980 → 8,474 GFLOPs, a 2.8× step. The load ratio alone predicts about 1.8×; the rest comes from the halved DRAM traffic and from those eight independent accumulators giving the scheduler something to overlap. And the obvious question now asks itself: if stretching a thread's territory in one direction helped this much, what happens if we stretch it in two?

9Kernel 5: 2D blocktiling, 68.7%

Look back at the sharing we found in section 8. There were two of them: outputs in a column share their B operand, and outputs in a row share their A operand. Kernel 4 exploited one and ignored the other. Its thread owns a tall thin strip, 8 × 1, so it gets column-sharing on B and nothing at all on A — which is exactly why it still pays 8 separate A loads for 8 multiply-adds.

So make the thread's territory square. Give it an 8 × 8 patch of C, and both sharings come into play at once.

Count that again, because it is the whole kernel. The patch has 64 outputs. Output (i, j) needs A-value number i and B-value number j. Every one of the 8 A-values pairs with every one of the 8 B-values, and each pairing is one of the 64 outputs. So 16 loads produce 64 multiply-adds.

The 8 × 8 grid of products formed by pairing every element of one vector with every element of another has a name: it is an outer product. Here it is with real numbers.

Now the mechanics, in three phases, because the structure of the inner loop matters as much as the count.

// phase 1 + 2: stage the operands into registers
for (uint i = 0; i < TM; ++i)
    regM[i] = As[(threadRow * TM + i) * BK + d];
for (uint i = 0; i < TN; ++i)
    regN[i] = Bs[d * BN + threadCol * TN + i];

// phase 3: 64 multiply-adds that touch no memory at all
for (uint m = 0; m < TM; ++m)
    for (uint n = 0; n < TN; ++n)
        threadResults[m * TN + n] += regM[m] * regN[n];

That staging step is what people mean when they say a kernel "caches in registers". Shared memory was a cache we filled by hand for the block; the register file is now a second cache we fill by hand for the thread. Each kernel in this post has added one more level of the hierarchy to the scheme.

The count: 16 loads, 64 multiply-adds, 0.25 loads per FMA. That is 4.5× better than kernel 4 and 8× better than kernel 3. Recall from section 7 that shared memory can supply roughly 32 words per cycle while the lanes can absorb 128 multiply-adds per cycle: at a quarter of a load per FMA, the supply finally matches the demand, and shared memory stops being the thing holding the kernel back.

The block reshapes again, for the same reason as before. With 64 outputs per thread, Boehm's kernel 5 uses a 128 × 128 tile of C and BK = 8:

• 128 × 128 = 16,384 outputs at 64 per thread gives 256 threads per block.
• The As tile is 128 × 8 = 1024 elements against 256 threads, so this time each thread copies four elements of A and four of B, not one.
• Shared memory: (128 + 128) × 8 × 4 bytes = 8 KB.
• Global reads per output fall again, from K/32 to K/64, because the block tile doubled in both directions.

Stepping back, kernels 3, 4 and 5 are one idea applied three times, and these two charts are the whole story of them.

Measured: 8,474 → 15,972 GFLOPs, and we are past two thirds of cuBLAS with a kernel that is still only about forty lines long. Which raises the obvious question: if a bigger thread tile keeps helping, why stop at 8 × 8? The answer is the subject of the next section.

10Occupancy: the Tetris, revisited

Kernel 5 asks each thread to hold 64 running totals, plus regM and regN, plus its indices. Those 64 accumulators must live in registers for the whole kernel — that is the entire point of them — so the thread needs at least 64 + 8 + 8 = 80 registers, and in practice a little more. Registers are a fixed, shared resource, so this is not free. To see what it costs, we need the packing game from part one, now played with a real kernel.

What occupancy means. An SM on this GPU can host at most 1536 threads, which is 48 warps. Occupancy is simply the fraction of those 48 slots that are actually filled. It matters for the reason part one gave: resident warps are what a scheduler switches to when the current warp stalls, so they are the raw material of latency hiding.

How many blocks fit. An SM accepts blocks until one of three budgets runs out — threads, registers, or shared memory — and whichever runs out first decides. For kernel 3 the numbers are exactly known, and they work out like this:

• Threads: 1536 ÷ 1024 = 1 block.
• Registers: 65,536 ÷ (1024 threads × 37 registers) = 1 block.
• Shared memory: 102,400 ÷ 9,216 bytes = 11 blocks.

The smallest wins: one block, which is 1024 threads, which is 32 warps of the 48 available — 66% occupancy. Shared memory was never close to being the constraint; threads and registers were, together.

So here is the uncomfortable fact. Kernel 5 has lower occupancy than kernel 3, and is five times faster. If occupancy were the thing to maximize, that could not happen. So it is not the thing to maximize, and it is worth being precise about why.

Latency hiding can come from two different places. The first is what part one described: many warps, so that when one stalls the scheduler switches to another. Call it thread-level parallelism. The second is many independent instructions inside a single warp, so that while one multiply-add is still in flight the next can already be issued, because it does not depend on the first. Call it instruction-level parallelism.

Kernel 3 has essentially none of the second kind. Its inner loop has one accumulator, and each iteration must wait for the previous one, so its only defence against latency is having lots of warps. Kernel 5 has 64 accumulators, and the 64 multiply-adds in its inner loop are completely independent of one another — the scheduler can keep many of them in flight from a single warp. It buys with ILP what kernel 3 had to buy with warps, and it gets a much better deal, because those same registers also eliminated seven eighths of the shared-memory traffic.

This also answers the question kernel 5 left hanging. Why not a 16 × 16 thread tile, with 256 accumulators? Because 256 registers per thread is beyond the 255 a thread may have at all, and long before that limit the register file would allow so few resident warps that nothing could hide even the shared-memory latency. The 8 × 8 tile sits near the sweet spot of that trade, which is precisely the kind of thing an autotuner is good at finding — as we will see in section 12.

11Kernel 6: vectorized memory access, 78.4%

From here the remaining gains are smaller and more specialist. We are close enough to the arithmetic limit that the wins come from shaving instructions rather than traffic.

A thread that loads four consecutive floats can do it as four 32-bit load instructions, or as one 128-bit load instruction. The bytes are identical; the instruction count is four times lower, and instruction issue is now part of what limits us. CUDA exposes this through the float4 type.

12Kernel 9: autotuning, 84.8%

By now the kernel has five free parameters: the block tile BM and BN, the k-chunk BK, and the per-thread tile TM and TN. Every choice trades against the others. Bigger tiles mean more reuse but more registers and more shared memory; more threads mean better latency hiding but a smaller share of the register file each.

Nobody can reason their way to the optimum, because it depends on the specific SM's register file size, its shared memory, its scheduler count, and the matrix shape. So you stop reasoning and start measuring: generate every legal configuration, compile it, run it, keep the winner.

The lesson is not "tuning gives 5%". It is that the optimum is a property of the specific chip, not of the algorithm, which is exactly why a portable library cannot simply hard-code one kernel. Hold that thought for cuBLAS.

13Kernel 10: warptiling, 93.7%

Until now our hierarchy has had two levels: the block tile, which maps to an SM, and the thread tile, which maps to a thread's registers. But the first post established that the hardware has a third level in between: the warp, the group of 32 threads that a single warp scheduler issues instructions for. Kernel 10 makes that level explicit in the code.

Why does merely relabelling who-computes-what buy 10%? Three reasons, all of them consequences of the warp being the real unit of execution. A warp now reads a compact, predictable region of shared memory, so its 32 lanes hit different memory banks instead of colliding. The values a warp needs stay resident in registers across more of the inner loop. And because each warp's work is independent of the others', the scheduler has cleaner runs of non-conflicting instructions to issue. 19,721 → 21,779 GFLOPs.

14cuBLAS, and the last 6%

Ten kernels in, we are at 93.7% of a library that has had two decades of engineers on it. The remaining gap is instructive, because it is not one clever trick.

cuBLAS is not one kernel. The library ships hundreds of megabytes of compiled code because it contains many specialised kernels and picks between them at runtime based on the matrix shapes. Our single kernel is tuned for one shape; cuBLAS carries a tuned kernel for practically every shape, which is the autotuning result of section 12 taken to its logical conclusion.

It changes strategy for awkward shapes. When the matrices are small or thin, there is not enough independent work to fill 84 SMs. cuBLAS then splits the K dimension across several blocks, has each compute a partial sum, and adds the partials in a second pass. More total work, better machine utilisation, faster wall clock.

It overlaps the loads with the arithmetic. Our kernels load a tile, sync, compute it, sync, load the next. cuBLAS prefetches: while computing on tile n it is already fetching tile n+1 into other registers, so the arithmetic units rarely stop. This is double buffering, and it is the main structural idea we skipped.

And even cuBLAS, at 23,250 GFLOPs, is only 78% of the chip's 30,000. The rest goes to the things no amount of software removes: the k-loop's own bookkeeping, the tails at the matrix edges, scheduling gaps, and the fact that the fp32 pipeline cannot be perfectly saturated by real code. The one way past it is to stop doing fp32 general arithmetic at all and use the tensor cores, which is a different instruction set and a different post.

15The whole ladder in one view

Reading down the third column, the same idea appears in different clothes each time. Kernel 2 is about the shape of a request. Kernel 3 is about not asking twice. Kernels 4 and 5 are about getting more arithmetic out of each thing you asked for. Kernel 6 is about asking in fewer instructions. Kernels 9 and 10 are about fitting the request to the specific machine. Underneath all of them is the sentence from section 2: the arithmetic never changes, only the cost of feeding it.

If you want to go further, the natural next steps are Boehm's own repository, where every kernel here is a compilable file you can profile yourself, and NVIDIA's CUTLASS, which is what a production-grade version of kernel 10 looks like when generalised to every shape and precision. And if any of the hardware words in this post felt like vocabulary rather than machinery, the first post builds all of them from a photograph of a dog.