05 — matvec: the 4-Output Kernel¶

Matrix-vector product y = A·x + y has a different optimal micro-kernel shape than matmul. Understanding why is the key insight of this family of levels.

The operation¶

A is m×k column-major, x has length k, y has length m:

\[ y_i \;\mathrel{+}=\; \sum_{p=0}^{k-1} A_{ip}\,x_p \]

FLOPs = \(2mk\) — and matvec is memory-bound by construction, because each element of A is used exactly once.

Why the kernel is `4×1`, not `4×4`¶

In matmul, a square tile lets you reuse both operands. In matvec, the operand x is a scalar per step that you want to broadcast across several outputs, while A is streamed. So the natural kernel is the 4×1 (four outputs):

for (p = 0; p < k; p++){
    double x_p = x[p];                 // broadcast scalar
    c_0 += A(0,p) * x_p;               // y_i
    c_1 += A(1,p) * x_p;               // y_{i+1}
    c_2 += A(2,p) * x_p;               // y_{i+2}
    c_3 += A(3,p) * x_p;               // y_{i+3}
}

At each step it:

loads one scalar x_p and reuses it for four rows,
reads a contiguous 4-run of A (a column slice, ideal for SIMD),
drives four independent accumulators (great ILP).

This is exactly the pattern a compiler auto-vectorizes into NEON/SSE/AVX.

The matvec ladder in brief¶

Level	Idea
optim0–1	naive row/col; `AddDot` extraction
optim2	4× unroll of the reduction loop
optim3	fused 4-output kernel (broadcast `x`)
optim4–5	register accumulators; pointer-advance `x` + unroll
optim6–8	the named `AddDot4x1`; register-bound
optim9–11	A pointers; reordering; portable vectorization
optim12–14	cache blocking; pack `A`; fully-packed pointer kernel

Portable vectorization

Unlike the matmul ladder's SSE intrinsics, matvec's optim11 is portable: it deliberately uses no intrinsics and lets the compiler auto-vectorize. The same source runs on macOS (ARM), Linux, and Windows.

Watch the shape¶

Run the ladder and notice:

optim2 (unrolling p on one output) barely helps — it keeps a single accumulator chain.
optim3 (four outputs) is the big jump — four independent chains.
Later levels mostly improve memory density (packing, blocking) rather than adding more arithmetic.

for i in 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14; do
  printf "matvec%-2s " $i; ./build/matvec$i
done