Kernel Engineering and the Dimensions of Kernel Optimization

Lesson 3 shrank the numbers themselves — fewer bits per value, the same information where it mattered. Lesson 4 leaves precision alone and asks a different question: now that we know what to compute and how small each value is, how do we make the arithmetic itself — the loop that multiplies and accumulates millions of times — run as close to what the silicon can actually deliver as possible?

1. From Graph to Kernel: Where the Real Arithmetic Happens

Every optimization in Lessons 1 through 3 happened above the loop that does the multiplying. The roofline model told you whether a kernel was starved for data or starved for ALUs. Graph-level optimization decided which tensor layout to use, how big a tile to carve out of a giant matmul, and which operators to weld together so an intermediate never round-trips through DRAM. Quantization decided how many bits each value gets. None of those three lessons wrote a single instruction that actually multiplies two numbers together.

Kernel engineering is where that stops being true. A kernel, in this sense, is the actual low-level compute routine that implements one operator — the concrete nested loop, running on real registers and a real cache hierarchy, that a graph compiler eventually lowers a MatMul or Conv2D node into. Everything upstream of this point was deciding what to compute and how much data to move; kernel engineering is where you decide the literal order of instructions that do it.

Take matrix multiplication, since it is both the simplest possible case and the one every other kernel in a neural network eventually reduces to in some form (convolution via im2col or implicit GEMM, attention projections, feed-forward layers). For C = A × B with A an M×K matrix, B a K×N matrix, and C the resulting M×N matrix, the textbook algorithm is three nested loops:

for (int i = 0; i < M; ++i)
    for (int j = 0; j < N; ++j)
        for (int k = 0; k < K; ++k)
            C[i][j] += A[i][k] * B[k][j];

This is mathematically complete. Compile it, run it, and it produces the exact right answer for every input. It is also, on any real CPU, GPU, or NPU built in the last two decades, dramatically slower than what the same silicon is capable of — often by one, sometimes two orders of magnitude. Not because the hardware is being asked to do extra arithmetic (it computes exactly 2 × M × N × K FLOPs either way), but because this particular arrangement of the same arithmetic makes catastrophically poor use of registers, caches, and vector units.

Closing that gap — turning a mathematically-correct triple loop into a kernel that saturates real hardware — is a discipline with its own name, its own literature, and its own set of largely independent levers. This post walks through nine of them. None of the nine changes what gets computed; every one of them changes how efficiently the exact same arithmetic gets executed.

#	Dimension	What problem it solves
1	Loop ordering	Which memory-access pattern the same arithmetic produces — sequential and cache-friendly, or scattered and cache-hostile
2	Blocking / tiling	Keeping working sets small enough to survive in a fast memory tier across many reuses (covered in depth in Lesson 2)
3	Cache reuse	Maximizing how many times data already sitting in L1/registers gets used before it's evicted, inside the inner loop itself
4	SIMD	Processing multiple data elements with a single instruction instead of one at a time
5	Register blocking	Keeping an entire output sub-tile resident in registers for the full reduction, avoiding repeated L1 traffic
6	Prefetching	Bringing data into cache before the pipeline actually stalls waiting for it
7	Alignment	Avoiding the extra cycles many SIMD load/store paths charge for misaligned addresses
8	Threading	Splitting the same arithmetic across cores without forcing them to fight over data
9	Specialized instructions	Using ISA-level fused and matrix-oriented instructions that do more useful work per instruction issued

The rest of this post takes each one in turn. Dimensions 4 (SIMD) and 2 (blocking) get comparatively short treatment here — SIMD because it earns a full dedicated post next (Lesson 4, Part 2), and blocking because Lesson 2 already derived it in depth as a graph-level technique; here it just needs to be placed correctly among its eight siblings. The other seven get the depth they deserve.

2. Loop Ordering: Same Math, Nine Different Memory-Access Patterns

2.1 The same arithmetic, six possible loop orders

The triple loop in Section 1 has three indices — i, j, k — and nothing in the mathematics of matrix multiplication says they have to be nested in that particular order. Any of the six permutations (ijk, ikj, jik, jki, kij, kji) computes the identical result, because addition is commutative and every (i, j, k) triple gets visited exactly once regardless of nesting order. What changes between permutations is purely the sequence in which memory addresses get touched — and that sequence is the single biggest factor in whether the kernel spends its time doing arithmetic or waiting on memory.

This is worth sitting with, because it is a genuinely strange fact if you haven't met it before: six pieces of code compute bit-for-bit the same output, cost the exact same number of floating-point operations, and yet can differ in wall-clock time by an order of magnitude. The FLOP count on a spec sheet or in a profiler's "operations" column tells you nothing about which of these six you're looking at.

2.2 The classic pair: ijk vs. ikj

The two most instructive orders to compare are ijk (the naive form from Section 1, with the reduction index k innermost) and ikj (with k in the middle and j innermost):

// ijk: reduction index (k) is innermost
for (int i = 0; i < M; ++i)
    for (int j = 0; j < N; ++j)
        for (int k = 0; k < K; ++k)
            C[i][j] += A[i][k] * B[k][j];
 
// ikj: reduction index (k) is in the middle, j is innermost
for (int i = 0; i < M; ++i)
    for (int k = 0; k < K; ++k)
        for (int j = 0; j < N; ++j)
            C[i][j] += A[i][k] * B[k][j];

Both loops execute M × N × K multiply-add operations. The difference is entirely in what the innermost loop touches as its index sweeps from 0 to its bound, and that difference is a direct consequence of how the matrices sit in memory.

Assume the standard C/C++ convention: matrices are stored row-major, meaning the elements of a single row are contiguous in memory, and moving to the next row means jumping forward by an entire row's worth of elements. For an N-wide matrix stored row-major, element B[k][j] sits at linear address base + (k * N + j) * elementSize — this is exactly the stride arithmetic the companion memory-hierarchy post in this site's Computer Architecture series works through for cache-line addressing, and it's the same row-major convention this series' Lesson 2 used when deriving NCHW strides.

2.3 Tracing the innermost loop's address stream

Under ijk, the innermost loop varies k while i and j are held fixed. Look at what each operand does as k increments by 1:

A[i][k]: address advances by exactly one element (elementSize bytes) — contiguous, stride-1, excellent locality.
B[k][j]: address advances by N * elementSize bytes, because incrementing the row index of a row-major matrix jumps a full row's worth of elements forward. For a 1024-wide float32 matrix, that's a 1024 * 4 = 4096-byte jump on every single step of the innermost loop.
C[i][j]: doesn't depend on k at all, so a compiler will typically hoist it into a single register accumulator for the duration of the inner loop and write it back once.

The B[k][j] access is the problem. A 4096-byte stride means every iteration of the innermost loop lands on a completely different cache line than the one before it — there is no spatial locality to exploit, because the next value you need is nowhere near the value you just used.

Under ikj, the innermost loop varies j while i and k are held fixed:

A[i][k]: doesn't depend on j at all — it's a single scalar value for the entire inner loop, naturally hoisted into a register and reused unchanged across every iteration.
B[k][j]: address advances by exactly one element per step — contiguous, stride-1.
C[i][j]: address advances by exactly one element per step — also contiguous, stride-1.

Every operand that varies in the ikj inner loop moves through memory the way it's actually laid out. Nothing jumps a row at a time; everything sweeps a row left to right, exactly matching how row-major storage puts data next to itself.

2.4 Quantifying the difference in cache-line terms

Put real numbers on it. Take a 64-byte cache line (the standard line size on essentially every modern x86 and ARM core) and float32 elements (4 bytes each), so one cache line holds 16 consecutive floats.

ijk inner loop over k (N = 1024 columns):
  B[k][j] address step = N * 4 bytes = 4096 bytes  ← jumps 64 lines every step
  Distinct cache lines touched to sweep k = 1024 (one new line per access)
  Useful bytes per fetched line = 4 out of 64       ← 6.25% line utilization
 
ikj inner loop over j (N = 1024 columns):
  B[k][j] address step = 4 bytes                     ← contiguous
  Distinct cache lines touched to sweep j = 1024 / 16 = 64
  Useful bytes per fetched line = 64 out of 64       ← 100% line utilization

To service one full sweep of 1024 B elements, ijk pulls in 1024 separate 64-byte cache lines and uses only 4 bytes out of each one; ikj pulls in 64 cache lines and uses every byte. That's a 16x reduction in the number of cache-line fetches needed to touch the same 1024 elements of B — exactly the 16-elements-per-line reuse factor you'd expect, made concrete instead of asserted.

It helps to see the two address streams side by side rather than just the summary numbers. Here is what the first few steps of each innermost loop actually touch, in units of the 64-byte cache line each address falls into:

One honest caveat worth naming, because it's the kind of detail that separates a real understanding from a slogan: ikj isn't a free lunch on the C side. Because C[i][j] now varies with the innermost j loop but the whole row gets revisited once per value of k, a full row of C is read-modify-written K times instead of being held in a single register for the entire reduction the way it was under ijk. In practice this cost is small — one row of C is only N * 4 bytes (4 KB for N = 1024), which sits comfortably in L1 and stays resident across the K repeated sweeps — but it's a real tradeoff, not a strict win on every axis. The B stride disaster ikj avoids is far larger than the C re-scan cost it accepts, which is why ikj-style orderings (and their kij sibling, which is nearly identical in spirit) are the ones real BLAS and NPU kernel authors build from, but it's worth knowing why the tradeoff nets positive rather than treating it as a rule with no exceptions.

2.5 Why this matters before any of the other eight dimensions

Loop ordering is listed first for a reason: every other dimension in this post — blocking, register allocation, SIMD, prefetching — operates on top of whatever access pattern the loop order already establishes. A perfectly vectorized inner loop that strides 4096 bytes between SIMD lanes' worth of data is still going to stall on cache misses; a loop with excellent locality but no vectorization at all is leaving throughput on the table but at least isn't fighting the memory system. Get the loop order wrong and every downstream optimization is polishing a fundamentally bad access pattern.

This is also precisely the transformation compilers call loop interchange — reordering a loop nest without changing its semantics, purely to improve the resulting memory-access pattern. It's one of the standard entries in the loop-nest-optimization toolkit alongside the tiling and interchange machinery described in Agner Fog's optimization manuals and in Intel's optimization reference manual, both cited at the end of this post.

3. Blocking / Tiling: The Same Idea, at Kernel Granularity

Lesson 2 derived tiling in full — how to size a tile against an SRAM budget, how arithmetic intensity scales with tile dimension, and the two failure modes of picking a tile too small or too large. That derivation doesn't need repeating here; it needs placing. Blocking is one of the nine dimensions a kernel author has to get right, sitting directly upstream of the next several sections.

The connection to loop ordering from Section 2 is direct: the ikj access pattern is good for a single pass over an entire row, but real matrices are far too large to keep any single dimension resident in L1 or L2 in full. Blocking restructures the loop nest itself — splitting each of M, N, and K into an outer "which tile" loop and an inner "within this tile" loop — so that a T×T (or more generally MC×KC/KC×NC) sub-block of A, B, and C is small enough to actually survive in a chosen cache tier across the many reuses that tile is going to receive. Everything Section 2 established about which loop order to use within a tile still applies once you're inside it; blocking is the layer that decides how big "within a tile" is allowed to be before the working set stops fitting where you want it to live.

Production GEMM libraries take this multi-level: OpenBLAS, following the GotoBLAS lineage, blocks separately for L3 (a KC × NC panel of B), L2 (an MC × KC panel of A), and L1 (a KC × NR micro-panel of B streamed during the innermost compute), stacking three tiers of the same idea from Lesson 2 rather than picking just one:

DRAM ──► L3 tile:  KC x NC panel of B          (megabytes, blocked for L3)
              │
              ▼
         L2 tile:  MC x KC panel of A          (hundreds of KB, blocked for L2)
              │
              ▼
         L1 tile:  KC x NR micro-panel of B    (tens of KB, blocked for L1)
              │
              ▼
      register tile: MR x NR block of C        (Section 6 — kept in registers)

The next dimension — cache reuse — is what happens inside whichever tile blocking has already sized correctly.

4. Cache Reuse: Locality Inside the Inner Loop

It's worth being precise about how this differs from blocking, because the two are easy to conflate and the task of a kernel author genuinely is different at each level.

Blocking answers the question "which chunk of the tensor is resident in a given memory tier right now?" — it's a decision about SRAM/L2-scale working sets, tens to hundreds of kilobytes, made once per tile.

Cache reuse, as its own dimension, answers a narrower question: "given that a tile is already resident, in what order — and how many times — does the inner loop touch each byte before it's done with that byte?" This operates at L1/register scale, kilobytes at most, and it's decided by the fine-grained structure of the innermost loops, not by the outer tiling scheme.

Section 2's ikj example is, underneath the loop-ordering framing, already a cache-reuse story: the entire reason ikj beats ijk is that it reuses A[i][k] across all N iterations of the inner loop instead of touching it once and moving on, and reuses each fetched cache line of B across 16 consecutive elements instead of one. Cache reuse is the property that access pattern is buying you; loop ordering is one of the tools that produces it.

The general principle generalizes past this one example: whenever a piece of data is going to be used more than once inside a bounded region of the computation, the kernel should be structured so that all of those uses happen while the data is still resident in the fastest tier that can hold it, rather than being separated by enough other work that the data gets evicted and has to be re-fetched. A concrete second example: in a blocked matmul's inner compute, a single row of the A micro-panel is typically reused against every column of the B micro-panel it's paired with before moving to the next row — the loop order inside the micro-kernel is chosen specifically to maximize how many multiply-adds happen per element loaded, which is exactly the arithmetic-intensity metric Lesson 1's roofline model formalized, now applied at the scale of a single resident tile instead of the whole kernel's DRAM traffic.

Register blocking, next, is what happens when you push this same idea to its logical extreme: instead of asking "how many times can I reuse data before it leaves L1," you ask "can I avoid L1 traffic for this data altogether by keeping it in registers for the computation's full duration."

5. SIMD: A Preview

SIMD — Single Instruction, Multiple Data — is the idea that one instruction can operate on several data elements packed into a single wide register simultaneously, rather than one scalar value at a time. A 256-bit AVX2 register holds eight float32 values; a single vfmadd instruction on that register performs eight independent multiply-adds in the time a scalar instruction performs one. NEON on ARM, AVX/AVX-512 on x86, and the vector units inside most NPU compute lanes are all instances of the same underlying idea: widen the data path, and get the compiler or kernel author to keep it full.

SIMD interacts with everything covered so far. It needs the loop order from Section 2 to have already produced contiguous, stride-1 access — a SIMD load instruction wants 8 or 16 consecutive elements sitting next to each other in memory, and a loop order that scatters those elements across a 4096-byte stride defeats vectorization before it starts, no matter how wide the register is. It needs the tile from Section 3 to already be resident nearby. And the register-blocked micro-kernel in the next section is, in every real high-performance GEMM implementation, written using SIMD registers rather than scalar ones — the "registers" a production micro-kernel keeps its output tile in are vector registers, each holding several output elements at once.

This post treats SIMD at the level of "what it is and why it matters" and leaves the full mechanics — vector register width across ISAs, masking, horizontal reductions, and a complete SIMD-vectorized GEMM microkernel — to Lesson 4, Part 2, where it gets the dedicated depth it deserves.

6. Register Blocking: Computing a Micro-Tile Entirely in Registers

6.1 Why even L1 is too expensive to hit on every step

Registers are the fastest storage a CPU core has — reading or writing one costs essentially nothing extra, because it's wired directly into the ALU's operand paths with no addressing, no tag check, and no possibility of a miss. L1 cache, one step down, is fast by any absolute standard (typically 4-5 cycles of load-use latency on a modern core) but it is not free, and more importantly, a core only has a small number of load/store execution ports it can issue through per cycle. A kernel that reloads its accumulator from L1 on every single multiply-add is spending load-port and store-port bandwidth — a genuinely scarce resource — on values that never actually left the chip's fastest storage tier in any meaningful sense.

Register blocking is the fix: pick a small output sub-tile — an MR × NR block of the result matrix, small enough to fit entirely inside the register file — and accumulate the entire reduction over K for that sub-tile using only registers, touching memory only to stream in the A and B operands and, once at the very end, to write the finished tile back out.

6.2 A 4x4 register-blocked micro-kernel

Concretely: pick MR = NR = 4, meaning the kernel computes a 4x4 tile of C — 16 output values — for the entire K-length reduction before writing anything back. The 16 accumulators live in registers (in a real vectorized kernel these would be a handful of SIMD registers; written as scalars here for clarity, the structure is identical to what BLIS and OpenBLAS actually generate):

// 4x4 register-blocked micro-kernel.
// Computes a 4x4 tile of C by accumulating over the full K
// reduction entirely in registers, touching memory only to
// stream in operands and, once, to flush the finished tile.
void micro_kernel_4x4(int K,
                       const float* A_panel, int lda,   // 4 rows x K, packed
                       const float* B_panel, int ldb,   // K x 4 cols, packed
                       float C[4][4]) {
    float c[4][4] = {0};   // 16 accumulators — lives in registers for the whole call
 
    for (int k = 0; k < K; ++k) {
        float a[4];
        float b[4];
        for (int i = 0; i < 4; ++i) a[i] = A_panel[i * lda + k]; // 4 loads: one A column
        for (int j = 0; j < 4; ++j) b[j] = B_panel[k * ldb + j]; // 4 loads: one B row
 
        for (int i = 0; i < 4; ++i)
            for (int j = 0; j < 4; ++j)
                c[i][j] += a[i] * b[j];   // 16 FMAs, register-to-register only
    }
 
    for (int i = 0; i < 4; ++i)           // one-time flush: 16 stores, total, ever
        for (int j = 0; j < 4; ++j)
            C[i][j] = c[i][j];
}

This is exactly the rank-1 update (outer-product) formulation that BLIS's kernel-authoring documentation and the GotoBLAS lineage OpenBLAS descends from both build their micro-kernels around: at each step of k, load a length-MR slice of the A panel and a length-NR slice of the B panel, then form their outer product and accumulate it directly into the MR × NR block of C sitting in registers. BLIS's own terminology calls this MR and NR pair the micro-kernel's register blocksize, chosen once by whoever writes the kernel for a given target and then compiled into that target's configuration.

Drawn out, one step of that outer-product update looks like this — a length-4 column of A, a length-4 row of B, and every pairwise product landing in its own accumulator cell:

Repeat that outer-product update once per value of k, and the 16 c[i][j] cells accumulate the full reduction without ever leaving the register file until the loop over k is finished.

6.3 Counting what register blocking actually saves

Compare against the naive alternative: a literal triple loop over the same 4x4xK sub-block, with no register discipline at all, where C[i][j] is read from and written back to memory on every single (i, j, k) iteration:

Naive (no register blocking), per k-step over the 16 (i, j) pairs:
  loads:  16 x (1 for C[i][j] + 1 for A[i][k] + 1 for B[k][j]) = 48
  stores: 16 x (1 for C[i][j])                                 = 16
  total memory ops per k-step                                  = 64
 
Register-blocked, per k-step:
  loads:  4 (A column) + 4 (B row)                              = 8
  stores: 0 (accumulators stay in registers until the very end) = 0
  total memory ops per k-step                                   = 8
 
One-time flush at the end of the K reduction: 16 stores, total, for the entire call.

For a realistic reduction depth of K = 256 (a plausible inner dimension for a hidden layer or attention projection):

Naive:            64 x 256              = 16,384 memory operations
Register-blocked: (8 x 256) + 16 (flush) =  2,064 memory operations

That's roughly an 8x reduction in memory traffic to produce the identical 16 output values, purely from refusing to round-trip the accumulator through memory on every multiply-add. Nothing about the arithmetic changed — it's still 4 × 4 × 256 = 4,096 multiply-adds either way — only where the intermediate values live while that arithmetic happens.

Scale MR and NR up (real BLIS and OpenBLAS kernels commonly use blocks in the range of 4x8, 6x8, 8x8, or larger depending on the target's register file size and vector width) and the same argument scales with it: a bigger register-resident output tile amortizes each loaded A and B element over more FMAs before that element is discarded, which is exactly the "operations per byte moved" quantity Lesson 1's roofline model already taught this series to care about — register blocking is arithmetic-intensity engineering applied at the register-file scale instead of the DRAM-traffic scale.

7. Prefetching: Hiding Latency Before the Pipeline Needs To Ask

7.1 Why out-of-order execution alone isn't enough

Modern CPU cores hide a great deal of memory latency automatically, through out-of-order execution: while one instruction waits on a load that missed cache, the core can keep issuing and completing later, independent instructions, as long as it can find enough of them within its reorder buffer — typically a few hundred entries deep on a modern high-performance core. That window is not infinite, and a full round trip to DRAM on a cache miss commonly costs on the order of 200 to 400+ cycles, depending on the platform. If the load that missed sits on the critical dependency chain — the very next multiply-add genuinely needs that value before it can proceed — and there isn't several hundred cycles' worth of independent, already-decoded work available to fill the gap, the reorder buffer drains and the core stalls waiting, full stop.

Prefetching is the deliberate strategy of bringing data into cache before it's actually needed, so that by the time the instruction that consumes it executes, the data is already sitting in L1 or L2 rather than triggering a fresh DRAM round trip on the critical path.

7.2 Hardware prefetchers

Most cores include hardware stream prefetchers that watch the stream of recent memory accesses, recognize sequential or fixed-stride patterns, and speculatively issue fetches for addresses the pattern predicts will be needed soon — entirely automatically, with no software involvement. For a kernel whose access pattern is already simple and regular (which, not coincidentally, is exactly what the loop-ordering work in Section 2 is trying to produce), the hardware prefetcher does a genuinely good job on its own.

It has real limits, though: hardware prefetchers typically detect only a handful of concurrent streams, have a bounded lookahead distance, and are tuned conservatively to avoid polluting the cache with data that turns out not to be needed. A kernel's access pattern that jumps between several packed panels — the A micro-panel, the B micro-panel, and the output tile, all being touched in an interleaved sequence inside a blocked, register-tiled inner loop — can be regular enough for a human (or a kernel author who designed the packing layout) to predict perfectly, while still being irregular enough, or spread across enough simultaneous streams, that a general-purpose hardware prefetcher misses some of it.

7.3 Software prefetch

That gap is what explicit software prefetch instructions exist to close: __builtin_prefetch in GCC/Clang, _mm_prefetch intrinsics on x86, and PLD/PLDW on ARM all let a kernel author insert an instruction that says, in effect, "start fetching this address now; I won't actually use the value for another several iterations." Because the kernel author knows the exact packing layout and loop structure — this is, after all, hand-authored or auto-generated kernel code, not a black box — they can issue a prefetch for the A and B panel data needed N iterations ahead, timed so the DRAM or L2-to-L1 latency is fully absorbed by useful FMA work happening in the meantime.

Laid out on a timeline, the goal is for the fetch of tile k+1's data to run entirely underneath the compute on tile k, so the core never actually stalls waiting:

This is exactly why high-performance GEMM and convolution kernels are among the most common places to see explicit prefetch intrinsics in real production code: the access pattern is regular enough to predict exactly, and the reduction loop is long enough to have real latency to hide behind.

8. Alignment: Why the Low Bits of an Address Matter

8.1 The cost of a misaligned SIMD access

A SIMD load or store instruction moves a fixed-width chunk of memory — 16, 32, or 64 bytes, depending on the vector width — in one instruction. Alignment means the starting address of that chunk is itself a multiple of the chunk's width (a 32-byte AVX load reading from an address that's a multiple of 32, for instance).

When an access is misaligned and, worse, straddles a cache-line boundary, the hardware may need to fetch two separate cache lines and merge the relevant bytes out of each before the vector register is filled — an extra step with a real, measurable cost that a naturally-aligned access simply never pays. On some architectures and older microarchitectures, misaligned SIMD loads historically required entirely separate, slower instruction forms, or were disallowed outright and would fault. Modern x86 cores (post-Nehalem, roughly) have narrowed this gap considerably for the common case — an unaligned load that happens not to cross a cache-line boundary costs little to nothing extra — but the boundary-crossing case still costs real cycles, and on many embedded and DSP-class targets, including a fair number of NPU and microcontroller-class SIMD units, strict alignment requirements remain the norm rather than the exception, sometimes enforced by a hard fault rather than a graceful slowdown.

8.2 How compilers and allocators keep this from being a problem

Because getting this right by hand on every buffer would be tedious and error-prone, the standard toolchain enforces alignment automatically in a few places:

Explicit alignment annotations — alignas(32) in C++11 and later, or __attribute__((aligned(32))) — tell the compiler to guarantee a variable or buffer starts on a chosen boundary.
Aligned allocation functions — posix_memalign, aligned_alloc, or platform-specific equivalents — let heap-allocated buffers (the packed A/B panels a GEMM kernel streams from, for instance) start on a SIMD-friendly boundary rather than wherever the general-purpose allocator happens to place them.
Arena and tensor allocators inside inference runtimes commonly round every buffer's base address up to 32- or 64-byte boundaries as a matter of policy, precisely so that every kernel downstream can assume alignment without checking, and so the compiler is free to emit the faster aligned-load instruction form wherever it can statically prove that guarantee holds.
The compiler itself only emits the faster aligned-access instruction form when it can prove alignment from context (a type's alignas, a known-aligned allocation, or an intrinsic that takes alignment as a compile-time guarantee); otherwise it's forced to conservatively emit the general, unaligned-safe form, even on hardware where that form is only marginally slower, because correctness has to come first.

The upshot for a kernel author: alignment is mostly a discipline of setting up buffers correctly once, at allocation time, rather than something to think about on every individual load inside the inner loop — but getting that setup wrong silently costs cycles on every single SIMD access the kernel ever issues against that buffer, which makes it a cheap, high-leverage thing to get right.

9. Threading: Splitting the Same Kernel Across Cores

9.1 Which loop to parallelize

A kernel that's already well-ordered, blocked, register-tiled, and vectorized within a single core is still using only one core. Threading partitions the same work across multiple cores — but which loop gets split matters enormously, for a reason that falls directly out of the loop-ordering discussion in Section 2.

The M and N loops — the ones that index into distinct, non-overlapping regions of the output C — are the right ones to parallelize. Split the output matrix into disjoint row or column ranges, hand one range to each thread, and every thread writes to memory no other thread ever touches. Zero synchronization is needed until every thread is done, because there is no shared mutable state between them at all.

The K loop — the reduction dimension — is the wrong one to parallelize naively, because every thread working on a slice of K would be accumulating into the same output element C[i][j]. Splitting K across threads either forces unsafe concurrent writes to a shared accumulator, or requires each thread to keep a private partial-sum tile and then pay an explicit reduction step at the end to combine them — extra synchronization and extra memory traffic that the M/N split simply never needs. This is exactly why production GEMM and NPU compiler backends parallelize at the level of output tiles: each thread or core is handed one or more MR × NR (or larger, block-level) tiles of C to own completely, computing that tile's full K reduction locally using everything from Sections 2 through 8, and touching no memory any other thread is writing.

Even with a correct, disjoint partitioning of output tiles across threads, one more hazard is specific to threading and easy to introduce by accident: false sharing. Cache-coherence protocols (MESI and its variants) track ownership of data at cache-line granularity, not at the granularity of individual variables. If two threads write to two different variables that happen to land in the same 64-byte cache line — for instance, adjacent elements of an output array partitioned so finely that a single cache line straddles the boundary between two threads' assigned ranges — the coherence protocol has no way to know the two threads aren't actually contending over the same data.

Every write from one thread invalidates the other thread's cached copy of that line, forcing a fresh fetch, even though the two threads never touch a single byte in common. The result is a workload that looks embarrassingly parallel on paper and performs worse than single-threaded in practice, with no incorrect output and no obvious symptom beyond "adding threads made it slower" — which is exactly why it's one of the first suspects to check when a parallel kernel underperforms.

The fix is straightforward once the hazard is named: give each thread's output region enough alignment padding that its boundary lands on a cache-line multiple rather than splitting a line down the middle, or size per-thread partitions in units of whole cache lines (or whole MR × NR register-blocked tiles, which are typically already several cache lines wide) so that no two threads' write targets can ever alias the same line. This is a natural, low-cost consequence of the tile-based partitioning Section 9.1 already recommends for correctness reasons — output tiles sized to match the register-blocking granularity from Section 6 tend to be comfortably larger than a single cache line, which sidesteps false sharing as a side effect of getting the parallelization strategy right for other reasons.

10. Specialized Instructions: FMA, Dot-Product, and MAC Arrays

10.1 Fused multiply-add

The single most impactful specialized instruction in this whole discussion is the humble fused multiply-add (FMA): an instruction that computes a * b + c as one operation, with a single rounding step at the end, instead of a separate multiply followed by a separate add with two rounding steps. This matters twice over. It's faster — one instruction and one pipeline slot instead of two — which directly multiplies the effective throughput of every register-blocked micro-kernel from Section 6, since the inner accumulation loop is nothing but a long chain of exactly this operation. And it's more numerically accurate, because a single correctly-rounded fused result avoids the extra rounding error a separate multiply-then-add would introduce at the intermediate step. x86's FMA3 extension (instructions like vfmadd231ps) and ARM NEON's fused multiply-accumulate forms both expose this directly, and essentially every SIMD-capable compiler will select the fused form automatically once the target ISA supports it and fast-math or explicit FMA semantics are permitted.

10.2 Dot-product and matrix-multiply-accumulate instructions

Beyond a single scalar FMA, modern SIMD and NPU instruction sets increasingly expose instructions that do considerably more accumulation work per instruction issued. ARM's dot-product extensions (SDOT/UDOT) compute a full int8 dot product across a vector's worth of elements and accumulate the result into a wider accumulator in one instruction, rather than requiring separate multiply and horizontal-add steps. x86's AVX-512 VNNI extension (VPDPBUSD and related instructions) does the equivalent for int8 operands directly in the mainstream server and desktop ISA. ARM's newer matrix-multiply instructions (SMMLA/UMMLA) and the Scalable Matrix Extension push this further still, computing small matrix-multiply-accumulate tiles — conceptually very close to the register-blocked micro-kernel structure from Section 6 — as a hardware primitive rather than a software loop. Intel's AMX extension takes a similar tile-based approach on the server side.

The throughline connecting all of these back to Section 6 is direct: every one of them is hardware doing, in a single instruction, some amount of the exact outer-product-accumulate structure a hand-written register-blocked micro-kernel implements in software with a loop over FMAs. A dot-product instruction is a tiny rank-1 update executed atomically; a matrix-multiply-accumulate instruction is a small MR × NR register-blocked tile computed in one shot, no software loop required at all.

10.3 Where this lands on an NPU

NPU MAC (multiply-accumulate) arrays take the same idea to its natural endpoint. Where a CPU's vector FMA or dot-product instruction still fetches operands, issues an instruction, and writes a result through a conventional register file, an NPU's MAC array is a grid of physical accumulator cells wired to compute an entire outer-product-style update per cycle in hardware, in parallel, across the whole array. The register-blocked tile from Section 6 isn't held in a handful of vector registers being updated by a software loop anymore; it's held in silicon that performs that update as its native operation.

The kernel-authoring problem on that target stops being "write an inner loop of FMAs efficiently" and becomes "get operands to the edges of the array fast enough, and in the right packed layout, to keep every one of its MAC cells busy every cycle." That reframing is precisely the job of an NPU compiler's kernel-selection and instruction-lowering passes: recognizing a matmul or convolution in the graph IR and lowering it to whichever native tile-multiply-accumulate primitive the target's MAC array exposes, rather than emitting a generic scalar or even SIMD loop and hoping the hardware's own dataflow control can still keep the array fed.

11. Putting the Nine Dimensions Together

None of these nine dimensions is sufficient on its own, and most of them compound rather than substitute for each other. A rough order in which a kernel author typically applies them:

1. Loop order (§2)       — get the base access pattern right; everything else inherits it
2. Blocking (§3)         — fit the computation to a chosen memory tier
3. Cache reuse (§4)      — maximize reuse of whatever is currently resident
4. Register blocking (§6)— push the innermost accumulator into registers entirely
5. SIMD (§5, next post)  — vectorize that register-resident computation
6. Prefetch + alignment (§7, §8) — make the memory side deliver its full throughput
7. Threading (§9)        — parallelize across cores once a single core is efficient
8. Specialized ISA (§10) — replace hand-rolled loops with native fused/MMA instructions

This is, not coincidentally, close to the actual order BLIS's kernel-authoring framework asks a porting engineer to work through, and close to the layered blocking-then-microkernel-then-vectorization structure inside OpenBLAS's real GEMM implementations. None of the nine dimensions is exotic in isolation — loop reordering, keeping things in registers, aligning buffers — but a kernel that gets all nine right, stacked correctly, is routinely an order of magnitude or more faster than the mathematically-identical naive triple loop this post opened with, without changing a single FLOP of arithmetic anywhere in the computation.