Complete Mathematics of "Optimization 7: The Fused 4×4 Micro-Kernel"¶
Optimization 7 fuses the sixteen scalar dot products of Optimization 6 into a single innermost loop over \(p\), keeping all sixteen \(C\) accumulators alive across that one loop. This is the classic GotoBLAS-style 4×4 micro-kernel.
1. The Mathematical Operation¶
Same block product as Optimization 6:
elementwise, for \(r, c \in \{0,1,2,3\}\):
2. The Key Structural Difference: One Loop, Sixteen Updates¶
2.1 Optimization 6 — sixteen loops¶
Each AddDot runs a separate loop of length \(k\), so sixteen full traversals:
for p (len k): C(i,j) += A(i,p)B(p,j)
for p (len k): C(i,j+1) += A(i,p)B(p,j+1)
... x16 total ...
for p (len k): C(i+3,j+3) += A(i+3,p)B(p,j+3)
The A and B panels are traversed almost entirely redundantly.
2.2 Optimization 7 — one loop¶
A single pass over \(p\) performs all sixteen updates:
for p = 0..k-1:
// fetch 4 A elements (column p, rows i..i+3)
a0 = A(i ,p); a1 = A(i+1,p); a2 = A(i+2,p); a3 = A(i+3,p);
// fetch 4 B elements (row p, cols j..j+3)
b0 = B(p,j); b1 = B(p,j+1); b2 = B(p,j+2); b3 = B(p,j+3);
// 16 rank-1 updates, all independent
c00 += a0*b0; c01 += a0*b1; c02 += a0*b2; c03 += a0*b3;
c10 += a1*b0; c11 += a1*b1; c12 += a1*b2; c13 += a1*b3;
c20 += a2*b0; c21 += a2*b1; c22 += a2*b2; c23 += a2*b3;
c30 += a3*b0; c31 += a3*b1; c32 += a3*b2; c33 += a3*b3;
In matrix language this is the rank-1 update at each \(p\):
where \(\mathbf{a}_p \in \mathbb{R}^{4}\) is the \(p\)-th column slice of \(\mathbf{A}\) and \(\mathbf{b}_p \in \mathbb{R}^{4}\) is the \(p\)-th row slice of \(\mathbf{B}\) (in column-major terms, the four values \(A(i:i+3,\,p)\) and \(B(p,\,j:j+3)\)).
3. Data-Fetch Savings¶
For one tile, Optimization 7 loads each of the \(4k\) A elements and \(4k\) B elements exactly once per \(p\)-pass (they are held in the eight working registers \(a_0..a_3, b_0..b_3\)). Compare:
| Resource per tile | Optim 6 | Optim 7 |
|---|---|---|
| A loads | \(16k\) (redundant) | \(4k\) |
| B loads | \(16k\) (redundant) | \(4k\) |
| Loop trips | \(16k\) | \(k\) |
| \(C\) traffic | \(16\cdot(2k)\) memory | \(16\) register accumulators, flushed once |
The inner loop is now compute-bound rather than load/loop-bound.
4. Instruction-Level Parallelism (ILP)¶
The sixteen accumulator chains \(c_{rc}\) are mutually independent — each is updated once per trip with fresh \(a_r,b_c\). This gives:
- dependency depth per chain: \(L_{\text{add}}\) (one FP add per step);
- sixteen interleavable chains: far more than the 8–10 needed to hide FP-add latency on a modern superscalar core with several FP pipes.
The multiply \(a_r \cdot b_c\) and the accumulate into \(c_{rc}\) form a fused multiply–add (FMA), giving \(16\) FMAs per trip for \(8\) loads of \(A\)+\(B\) — an arithmetic intensity inside the kernel of \(2\) FLOP per loaded element, at the register level.
5. Residual Limitation: Strided A¶
Inside the \(p\) loop, the column-major address of \(A(i+r,p)\) is:
Varying \(p\) means a stride of \(\text{lda}\), and the four rows \(r=0..3\) are at base offsets \(0,1,2,3\) within each column. So the \(a_0..a_3\) loads are strided by \(\text{lda}\), which is cache-unfriendly. Later optimizations fix this by packing A (Optimization 13) into a contiguous buffer. B's four values \(B(p,j:j+3)\) are contiguous (stride 1) in column-major layout.
6. Summary¶
| Metric | Optim 6 | Optim 7 |
|---|---|---|
| Kernel shape | \(4\times4\) | \(4\times4\) |
| Loops per tile | 16 | 1 |
| A+B loads per tile | \(32k\) | \(8k\) |
| Independent accumulators | 0 in loop | 16 |
| FLOPs | \(2mnk\) | \(2mnk\) |
| In-kernel arithmetic intensity | ~1 FLOP/elem | 2 FLOP/elem |
| A access | strided | strided (fixed later) |
The essential math is the rank-1 update accumulation
which keeps sixteen independent accumulator registers busy and makes the inner loop compute-bound, setting up the fully register-tiled kernel of Optimization 8.