Complete Mathematics of "Optimization 11: Portable, Auto-Vectorizable 4-Output Kernel"

Optimization 11 is the vectorization step — but with a critical difference from the matmul ladder: it contains no architecture-specific SIMD intrinsics. Instead, it writes the AddDot4x1 kernel in a form the compiler can auto-vectorize into whatever SIMD ISA the target supports — NEON on Apple Silicon/ARM, SSE/AVX on x86. This keeps the code portable across macOS, Linux, and Windows while still capturing vector throughput.

---

1. The Mathematical Operation (unchanged)

For \(r=0,1,2,3\):

\[ y_{i+r} \;\mathrel{+}=\; \sum_{p=0}^{k-1} A_{i+r,p}\,x_p. \]

The kernel is textually identical to Optimization 9 (pointer-based, register accumulators) plus a documented, compiler-friendly structure:

register double c_0, c_1, c_2, c_3, x_p;
double *a_0_ptr, *a_1_ptr, *a_2_ptr, *a_3_ptr;

for (p = 0; p < k; p++){
    x_p = x[p];
    c_0 += *a_0_ptr * x_p;   // y[i]   stream
    c_1 += *a_1_ptr * x_p;   // y[i+1] stream
    c_2 += *a_2_ptr * x_p;   // y[i+2] stream
    c_3 += *a_3_ptr * x_p;   // y[i+3] stream
    a_0_ptr += lda; a_1_ptr += lda; a_2_ptr += lda; a_3_ptr += lda;
}

---

2. Why This Structure Auto-Vectorizes

The inner loop is a textbook SIMD pattern:

1. A broadcast scalar x_p is multiplied into four independent accumulators c_0..c_3.
2. The compiler recognizes the four independent accumulation chains as four vector lanes of a packed register.

Conceptually the compiler transforms the scalar loop into a vector operation:

\[ \begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \end{bmatrix} \;\mathrel{+}=\; \begin{bmatrix} *a_0Ptr \\ *a_1Ptr \\ *a_2Ptr \\ *a_3Ptr \end{bmatrix} \cdot\, (\text{broadcast } x_p). \]

On AVX2/NEON with 256/128-bit registers this maps to one packed vector FMA per \(p\)-step (4 doubles in 256-bit AVX, or pairs of SSE/NEON ops for 128-bit). The four c_r lanes are exactly the four-element vector width needed.

---

3. Portability: The Crucial Choice

The matmul ladder's Optimization 11 uses hand-written SSE intrinsics (__m128d, _mm_load_pd, _mm_loaddup_pd) which are x86-64 only and would fail to compile on Apple Silicon (ARM) or other non-x86 platforms.

This matvec Optimization 11 deliberately avoids intrinsics:

Feature	Intrinsic SSE (matmul 11)	Auto-vectorized (matvec 11)
SIMD ISA	SSE/AVX only (x86)	NEON, SSE, AVX, etc.
Apple Silicon / ARM	❌ won't compile	✅ auto(NEON)
Linux x86	✅	✅ auto(SSE/AVX)
Windows x86	✅	✅ auto(SSE/AVX)
Code portability	hardware-locked	fully portable

The compiler (with -O2/-O3/-ftree-vectorize on GCC, or -O3/-O2 + vectorizer on Clang, or /O2 on MSVC) turns the four independent chains into packed instructions. This realizes "portable vectorization" — same source, best SIMD on each platform.

---

4. Why matvec Especially Benefits

Unlike matmul, matvec reads each \(A\) element only once — there's no matrix–matrix reuse to block for. The performance ceiling is set by vectorized streaming: packing as many \(A\) elements per load as possible and keeping the FMAs maximally parallel. Four rows per kernel call give a 4-wide natural vector, matching typical SIMD width without any manual unrolling tricks.

---

5. FLOP Accounting (unchanged)

\[ \text{FLOPs} = \frac{m}{4}\cdot (4\cdot 2k) = 2mk. \]

Per \(p\)-step the (conceptual) vector FMA does 4 scalar FLOPs \(\times\) 2 (mult+add) \(= 8\) scalar-equivalent FLOPs.

---

6. Summary

Metric	Optim 9 (scalar)	Optim 11 (auto-vectorized)
SIMD intrinsics	none	none (portable)
ISA support	scalar only	NEON/SSE/AVX (by compiler)
Vector width	1 double	up to 4 doubles
Accumulators	4 scalar	4 vector lanes
FLOPs	\(2mk\)	\(2mk\)
Portability	✅	✅

Optimization 11 reaches vector throughput without sacrificing portability — the defining design goal of this matvec ladder. It relies on the compiler to vectorize the four independent accumulator chains, so the same source runs fast on ARM, Intel, and AMD across all three major OSes.