Complete Mathematics of "Optimization 14: Fully Packed, Pointer-Driven `4×1` Kernel"¶

Optimization 14 is the culmination of the matvec ladder. It packs each 4×k A-panel into a contiguous buffer and makes the AddDot4x1 kernel consume that buffer with pure += 4 pointer advances — so the leading dimension collapses to 4 and no stride-lda arithmetic remains anywhere in the hot loop. Combined with the mc×kc blocking (Optimization 12), this is the fully mature matvec pipeline.

1. The Mathematical Operation (unchanged)¶

For the 4×k panel of rows \(i..i+3\):

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

equivalently,

\[ \begin{bmatrix} y_i \\ y_{i+1} \\ y_{i+2} \\ y_{i+3} \end{bmatrix} \;\mathrel{+}=\; \mathbf{A}_{4\times k}\; x \;=\; \sum_{p=0}^{k-1} x_p\, \begin{bmatrix} A_{i,p} \\ A_{i+1,p} \\ A_{i+2,p} \\ A_{i+3,p} \end{bmatrix}. \]

2. Packing: The `4`-Interleaved Layout¶

PackMatrixA writes the panel columns interleaved so each column's 4 rows are contiguous:

\[ \texttt{packedA}[\,4p + r\,] = A(i+r,\;p), \qquad 0\le p<k,\;0\le r<4. \]

The packed element for column \(p\), row \(r\) is at offset \(4p+r\).

3. The Pointer-Driven Kernel¶

double *a_0_ptr = a;        // row 0 of col 0
double *a_1_ptr = a + 1;    // row 1 of col 0
double *a_2_ptr = a + 2;    // row 2 of col 0
double *a_3_ptr = a + 3;    // row 3 of col 0

for (p = 0; p < k; p++){
    x_p = x[p];
    c_0 += *a_0_ptr * x_p;
    c_1 += *a_1_ptr * x_p;
    c_2 += *a_2_ptr * x_p;
    c_3 += *a_3_ptr * x_p;

    a_0_ptr += 4;   // next column: same 4-interleaved stride
    a_1_ptr += 4;
    a_2_ptr += 4;
    a_3_ptr += 4;
}

3.1 Why `+= 4`, not `+= lda`¶

In the packed layout, consecutive columns differ by 4 doubles. So the four pointers advance by a constant 4 each step instead of by lda. The kernel is called with the packed buffer and lda = 4:

AddDot4x1(k, &packedA[i*k], 4, x, &y[i]);

The lda parameter is effectively the interleave factor (4), so all address arithmetic becomes a uniform += 4.

4. Address Arithmetic: Totally Uniform¶

Scheme	Step to next column	Per-step cost
Strided (Optim 9–13)	all pointers`+= lda` (16-KB jump)	large jump, scattered loads
Packed (Optim 14)	all pointers`+= 4`	tiny constant, fully sequential

Every load is now sequential (dense packed buffer), giving:

Perfect hardware prefetch,
Minimal TLB pressure,
Full cache-line utilization,
And the compiler can still auto-vectorize the four independent accumulator lanes (as in Optimization 11).

5. The Full Pipeline¶

for (p = 0; p < k; p += kc) {        // reduce by kc-panels
    pb = min(k-p, kc);
    for (i = 0; i < m; i += mc) {    // rows by mc-blocks
        ib = min(m-i, mc);
        InnerKernel(ib, pb, &A(i,p), lda, &x[p], &y[i]);
    }
}

// InnerKernel:
double *packedA = calloc(m*k, sizeof(double));
for (i = 0; i < m; i += 4){
    PackMatrixA(k, &A(i,0), lda, &packedA[i*k]);
    AddDot4x1(k, &packedA[i*k], 4, x, &y[i]);
}

Blocking confines the working set (Optimization 12); packing makes reads sequential (Optimization 13); the pointer kernel consumes them with += 4 (this level).

6. Mathematical Equivalence of Packing¶

Packing is a bijection (each element copied exactly once, in a known order), so:

\[ \sum_p x_p\, \texttt{packedA}[4p : 4p+3]^{T} \;=\; \sum_p x_p \,[A_{i,p},A_{i+1,p},A_{i+2,p},A_{i+3,p}]^{T}. \]

The computed \(y\) is identical to every earlier level.

7. FLOP Accounting¶

Per kernel call: \(4\cdot2k = 8k\) FLOPs; over \(m/4\) calls:

\[ \text{FLOPs} = \frac{m}{4}\cdot 8k = 2mk. \quad\checkmark \]

8. Summary — The Complete Matvec Ladder¶

Level	Core idea	A access	x access
optim0	naive row/col	strided/row	contiguous
optim1	`AddDot` extraction	strided (lda)	contiguous
optim2	4× reduction unroll	strided	contiguous
optim3	fused 4-output kernel	contiguous 4-run	broadcast
optim4	register accumulators	contiguous 4-run	broadcast
optim5	x-pointer + reduce unroll	strided	contiguous
optim6	`AddDot4x1` (via AddDot)	strided	contiguous
optim7	fused`AddDot4x1`	contiguous 4-run	broadcast
optim8	register-bound 4 accums	contiguous 4-run	broadcast
optim9	4 A pointers (+= lda)	strided columns	broadcast
optim10	statement reorder (+OpenMP hooks)	strided columns	broadcast
optim11	portable auto-vectorization	strided columns	broadcast
optim12	cache blocking (mc×kc)	panels	kc-slices
optim13	pack A	packed, kernel still strided	broadcast
optim14	pack A + raw-pointer kernel (`+= 4`)	packed sequential	broadcast

9. Conclusion¶

Optimization 14 assembles every technique in the ladder into its final form for

\[ y \;\leftarrow\; A\,x + y: \]

the fused, register-bound 4×1 micro-kernel (4 independent accumulator lanes, x broadcast),
cache blocking (mc×kc) for cache/TLB locality,
packing of each A-panel into a dense 4-interleaved buffer,
and a pure pointer-driven kernel whose only per-step work is four += 4 increments and four sequential loads.

Everything upstream of the arithmetic is set up to keep the memory subsystem streaming densely and the FP/SIMD pipes fully occupied — the practical limit for a memory-bound matrix–vector product.

Complete Mathematics of "Optimization 14: Fully Packed, Pointer-Driven 4×1 Kernel"¶