Complete Mathematics of "Optimization 3: The Fused `AddDot1x4` 4-Output Kernel"¶

Optimization 3 introduces the canonical matvec micro-kernel: AddDot1x4, which computes four output elements \(y_i, y_{i+1}, y_{i+2}, y_{i+3}\) simultaneously, in a single loop over \(p\). This is the matvec analogue of matmul's fused 4×4 kernel, but with a shape tailored to the vector product.

1. The Mathematical Operation¶

For a group of four consecutive rows \(i, i+1, i+2, i+3\), the kernel computes:

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

which in suppressed indices reads:

\[ \begin{bmatrix} y_i \\ y_{i+1} \\ y_{i+2} \\ y_{i+3} \end{bmatrix} \;\mathrel{+}=\; \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}. \]

This is the rank of the vector product: at each \(p\), the scalar \(x_p\) is multiplied into a four-row column slice of \(A\), and the results accumulate into four distinct outputs.

2. The Kernel Code¶

void AddDot1x4(int k, double *a, int lda, double *x,
               double *y0, double *y1, double *y2, double *y3){
    for (p = 0; p < k; p++){
        double x_p = x[p];
        y0[0] += A(0,p) * x_p;   // row i
        y1[0] += A(1,p) * x_p;   // row i+1
        y2[0] += A(2,p) * x_p;   // row i+2
        y3[0] += A(3,p) * x_p;   // row i+3
    }
}

called from matvec as:

for (i = 0; i < m; i += 4)
    AddDot1x4(k, &A(i,0), lda, x, &y[i], &y[i+1], &y[i+2], &y[i+3]);

3. Key Structural Insight: "Broadcast x, Walk 4 Contiguous Rows"¶

In column-major layout the entries of one column \(p\) across rows \(i..i+3\) are contiguous in memory:

\[ A(i,p),\;A(i+1,p),\;A(i+2,p),\;A(i+3,p) \quad = \quad a[p\,\text{lda}+i],\; a[p\,\text{lda}+i+1],\; a[p\,\text{lda}+i+2],\; a[p\,\text{lda}+i+3]. \]

So at each \(p\):

load the scalar \(x_p\) once — it is broadcast and reused by all four rows;
read the contiguous 4-run \(A(i:i+3,\ p)\).

The innermost loop is a perfect SIMD shape: one scalar load (_dup/broadcast x) times a contiguous 4-element vector load (the A column slice), producing a 4-wide vector accumulation. This is why the natural matvec kernel is "4 outputs per pass," not a square tile.

4. ILP and Independent Chains¶

The four accumulators y0..y3 are distinct output locations, hence four independent dependency chains. Each chain updates once per \(p\)-step:

\[ y_{r}^{(p+1)} \;\leftarrow\; y_{r}^{(p)} + A_{i+r,p}\,x_p. \]

With four independent chains the FP units can overlap the multiply-add latency of one row with the arithmetic of the others — a substantial ILP win over the single-accumulator Optimization 2.

5. Data-Reuse Accounting per Kernel Call¶

For one \(4\times k\) panel of \(A\) and the vector \(x\):

Quantity	Accesses
\(x_p\) loaded	\(k\) (each read once, reused 4×)
\(A\) elements read	\(4k\) each once
outputs produced	4

The vector \(x\) is fully reused across the four outputs of the group, but each \(A\) element is still read exactly once — consistent with matvec's memory-bound nature.

6. Memory Access vs Previous¶

Access	Optim 2 (1 output)	Optim 3 (4 outputs)
A read per trip	4 strided (row)	4 contiguous (column slice)
x	4 contiguous	1 broadcast
A access pattern	strided `lda`	contiguous column slice

A notable improvement: the kernel now reads \(A\) in contiguous 4-runs (full cache-line-friendly) rather than chasing the full row at stride lda. Even though successive \(p\) still jump by lda, each step pulls a 4-double contiguous run instead of a single isolated element, improving cache-line utilization.

7. FLOP Accounting¶

Per kernel call: \(4\) rows \(\times\) (\(k\) mult + \(k\) add) \(= 8k\) FLOPs. Number of calls \(= \frac{m}{4}\).

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

8. Summary¶

Metric	Optim 2	Optim 3
Outputs per kernel	1	4
Loop trips per group	\(k\) (×1)	\(k\) (for 4 rows)
Independent chains	1	4
A read	strided row	contiguous 4-run column slice
x reuse	none across outputs	broadcast across 4 rows
FLOPs	\(2mk\)	\(2mk\)

Optimization 3 is the defining turn in the matvec ladder: it establishes the 4-output kernel that (a) reads \(A\) in contiguous column slices and (b) broadcasts \(x\) across four independent accumulators — exactly the structure the compiler can auto-vectorize. Subsequent steps hold these four accumulators in registers and optimize the addressing of \(A\).

Complete Mathematics of "Optimization 3: The Fused AddDot1x4 4-Output Kernel"¶