Complete Mathematics of "Optimization 4: Register-Bound Accumulators"¶

Optimization 3 updated the four elements of \(C\) through memory on every iteration of \(p\):

C(0,c) += A(0,p) * B(p,c);   // read C from memory, add, write C back to memory

Optimization 4 captures the four accumulators in CPU registers for the entire \(p\)-loop, and only writes them back to \(C\) once, after the loop. This dramatically reduces memory traffic to the output matrix.

1. The Mathematical Operation (unchanged)¶

The arithmetic is identical to Optimization 3:

\[ C_{i,\,j+c} \;\mathrel{+}=\; \sum_{p=0}^{k-1} A_{i,p}\, B_{p,\,j+c}, \qquad c = 0,1,2,3. \]

The change is purely about where the partial sum lives during the loop.

2. Register Accumulation¶

2.1 Before (Optimization 3) — read-modify-write to memory¶

Let \(c_{pc}\) denote the value of \(C_{i,j+c}\) as it resides in main memory. Invoice iteration the kernel does:

\(\text{load } c_{0c}\) from memory into a register,
\(\text{compute } c_{0c} + A_{i,p}B_{p,j+c}\),
\(\text{store } c_{0c}\) back to memory.

Because \(A_{ip}B_{pj}\) must accumulate over all \(p\), the CPU performs a memory round-trip every iteration, for every one of the four outputs. For a full kernel this is \(4k\) loads plus \(4k\) stores of \(C\).

2.2 After (Optimization 4) — accumulate in registers¶

We introduce four real-valued accumulator variables held in registers:

\[ \gamma_0,\; \gamma_1,\; \gamma_2,\; \gamma_3. \]

They are initialized to zero, updated in the inner loop, and flushed at the end:

\[ \gamma_c^{(0)} = 0, \qquad \gamma_c^{(p+1)} = \gamma_c^{(p)} + A_{i,p}\, B_{p,j+c}, \]

\[ C_{i,j+c} \;\mathrel{+}=\; \gamma_c^{(k-1)}. \]

No access to the address of \(C(i,j+c)\) happens inside the \(p\)-loop.

3. Memory-Traffic Reduction for \(C\)¶

Let \(T_C(\text{optim 3})\) and \(T_C(\text{optim 4})\) be the number of memory transfers to \(C\) per \(1\times4\) tile.

Metric	Optim 3	Optim 4
Loads of\(C\) in loop	\(4k\)	\(0\)
Stores of\(C\) in loop	\(4k\)	\(0\)
Loads of\(C\) after loop	\(0\)	\(0\) (uses `+=`)
Stores of\(C\) after loop	\(0\)	\(4\)
Total \(C\) traffic	\(8k\)	\(4\)

Because \(k\) is large (here \(2000\)), this removes roughly \(8k \approx 1.6 \times 10^4\) memory operations per tile — a huge win for a memory-latency-limited inner loop.

4. Dependency-Chain / Latency Analysis¶

4.1 The recurrence¶

The accumulation is a chain recurrence:

\[ \gamma_c^{(p+1)} = \gamma_c^{(p)} + A_{i,p} B_{p,j+c}. \]

Each step depends on the previous step's value of \(\gamma_c\). A single such chain has a latency of one floating-point add per iteration (the multiply can happen off the critical path).

Let:

\(L_{\text{add}}\) = FP add latency (e.g., 4 cycles),
\(L_{\text{mem}}\) = memory latency if \(\gamma\) were in memory (e.g., 30–100 cycles).

Optimization 3's chain goes through memory, so its period is dominated by \(L_{\text{mem}}\). Optimization 4's chain stays in registers, so its period is \(L_{\text{add}}\) (a few cycles).

Since the four chains \(c=0..3\) are independent, the CPU can interleave them. With enough independent chains the FP units stay saturated:

\[ \text{throughput} \approx \frac{\text{FPUs}}{\max(L_{\text{add}}, \text{issue width})}. \]

5. The `register` Keyword¶

The code hints the compiler with:

register double c_00, c_01, c_02, c_03, a_0p;

Here \(a_{0p}\) also lives in a register and holds the single row element \(A(0,p)\) reused for all four products (just as in Optimization 3's fused kernel). The register keyword is only a hint; with -O2/-O3 the compiler will generally keep these hot variables in registers anyway. Its explicit use documents intent and, on some older compilers, changes register allocation for the better.

6. Summary of the Math¶

Quantity	Optim 3	Optim 4
FLOPs	\(2mnk\)	\(2mnk\)
\(C\) traffic per tile	\(8k\)	\(4\)
Accumulator location	memory	registers
Inner-loop latency bound	memory latency	FP-add latency

The central mathematical idea is deferring the final memory write: instead of writing an amortized result to memory every step, we carry partial sums in registers and commit them once. This converts the inner loop from a memory-latency-bound operation into a compute-bound one.