Complete Mathematics of "Optimization 11: SSE Vectorized 4×4 Kernel"¶

Optimization 11 switches from scalar floating-point arithmetic to SIMD vector registers using x86-64 SSE intrinsics. Each 128-bit register (__m128d) holds two doubles. This lets the kernel process two rows of the 4×4 tile per vector operation, doubling the arithmetic throughput per instruction while keeping exactly the same number of mathematical operations.

1. The Mathematical Operation (unchanged)¶

Still the block product:

\[ \mathbf{C} \;\mathrel{+}=\; \mathbf{A}_{4\times k}\,\mathbf{B}_{k\times4} \;=\; \sum_{p=0}^{k-1} \mathbf{a}_p\, \mathbf{b}_p^{T}. \]

But now each element-pair is computed together in one 128-bit vector lane.

2. The Vector Data Layout¶

A v2df_t union packs two doubles into one __m128d:

\[ \mathrm{vreg} = \big( d_0,\; d_1 \big), \qquad d_0,d_1 \in \mathbb{R}. \]

2.1 Vectorizing the rows of A¶

Column-major storage means rows \(i\) and \(i+1\) of a column \(p\) are adjacent in memory, so one aligned load grabs two rows at once:

\[ \mathbf{a}_{0,1}(p) = \big( A(i,p),\; A(i+1,p) \big) = \texttt{\_mm\_load\_pd(\&A(0,p))}, \]

\[ \mathbf{a}_{2,3}(p) = \big( A(i+2,p),\; A(i+3,p) \big) = \texttt{\_mm\_load\_pd(\&A(2,p))}. \]

2.2 Duplicating the scalars of B¶

_mm_loaddup_pd loads a scalar and broadcasts it into both lanes:

\[ \mathbf{b}_{c}(p) = \big( B(p,j+c),\; B(p,j+c) \big) = \texttt{\_mm\_loaddup\_pd(b\_pc\_ptr++)}. \]

This is the SIMD form of the "reuse \(B\) across two rows" idea: the same column value multiplies both row-accumulators.

2.3 Vector accumulators¶

The sixteen scalar accumulators become eight vector accumulators, each holding two rows of one column:

\[ \mathbf{\Gamma}_{0c} = (\gamma_{0c},\, \gamma_{1c}),\qquad \mathbf{\Gamma}_{2c} = (\gamma_{2c},\, \gamma_{3c}),\qquad c=0,1,2,3. \]

3. The Vectorized Rank-1 Update¶

Each iteration of \(p\) performs the vector rank-1 update:

\[ \mathbf{\Gamma}_{0c} \;\mathrel{+}=\; \mathbf{a}_{0,1}(p)\cdot \mathbf{b}_c(p), \qquad c = 0,\dots,3; \]

\[ \mathbf{\Gamma}_{2c} \;\mathrel{+}=\; \mathbf{a}_{2,3}(p)\cdot \mathbf{b}_c(p), \qquad c = 0,\dots,3. \]

For example, the pair of updates for \((r,c)=(0,0)\) and \((1,0)\) collapses into one vector multiply-add:

\[ \begin{bmatrix} \gamma_{00} \\ \gamma_{10} \end{bmatrix} \;\mathrel{+}=\; \begin{bmatrix} A(i,p) \\ A(i+1,p) \end{bmatrix}\cdot B(p,j). \]

4. FLOP and Instruction Accounting¶

Each vector operation handles 2 scalar FLOPs per element × 2 elements = 4 FLOPs per vector FMA.

Scalar FLOPs per kernel: \(2\cdot16\cdot k = 32k\).
Vector FMAs per \(p\)-step: 8 (two row-groups × four columns).
Vector FMA count per kernel: \(8k\).

Since \(\frac{32k}{8k}=4\), each vector FMA does 4 scalar-equivalent FLOPs, exactly matching the two-element lanes.

Throughout: total work remains

\[ \text{FLOPs} = \frac{mn}{16}\cdot 32k = 2mnk. \]

5. Register Pressure Improves¶

The 16 scalar accumulators needed 16 registers at the limit of the budget. Vectorizing folds pairs into 8 vector registers, freeing registers for operands and pointers:

State	Scalar (Optim 8)	Vector (Optim 11)
Accumulators	16	8
A operands	4	2
B operands	4	4
Pointers	4	4

6. Alignment Requirement¶

_mm_load_pd requires 16-byte alignment. For a column-major matrix with leading dimension 2000 (even), the two double-groups (&A(0,p), &A(2,p)) land on 16-byte-aligned addresses only if the base allocation and offset align. In practice the union type ensures the accumulator registers align; A-loads with even lda keep the offset p*lda even (8-byte multiples give 16-byte alignment when lda is even). This is why the code uses the aligned load.

7. Summary¶

Metric	Optim 10 (scalar)	Optim 11 (SSE)
FLOPs	\(2mnk\)	\(2mnk\)
Vector width	1 double	2 doubles
Accumulators	16 scalar	8 vector
FMAs per kernel	\(16k\)	\(8k\)
FLOPs per FMA	2	4

The mathematics is identical; the change is that every pair of output rows is now computed in a single 128-bit vector lane, halving the instruction count and (on machines with full SIMD FMA throughput) roughly doubling achievable GFLOPS, all without altering a single computed value.