Complete Mathematics of the Naive Matrix–Vector Multiplication (matvec)¶

1. Problem Setup¶

Given a matrix and two vectors:

\[ A \in \mathbb{R}^{m \times k}, \qquad x \in \mathbb{R}^{k}, \qquad y \in \mathbb{R}^{m}, \]

the operation computed is:

\[ \boxed{\; y \;\leftarrow\; A\,x \;+\; y \;} \]

elementwise:

\[ y_i \;\leftarrow\; \sum_{p=0}^{k-1} A_{i,p}\, x_p \;+\; y_i, \qquad i = 0,1,\dots,m-1. \]

For the main() routine in this project:

\[ m = k = 2000, \]

so \(A \in \mathbb{R}^{2000\times2000}\), \(x \in \mathbb{R}^{2000}\), \(y \in \mathbb{R}^{2000}\).

2. Row-Major vs Column-Major Storage (`gemv_row` vs `gemv_col`)¶

The naive version provides both storage conventions.

2.1 Row-major (`gemv_row`)¶

Element \((i,j)\) is at linear address:

\[ \operatorname{addr}(A_{ij}) = i\cdot\text{lda} + j, \]

matching the macro \(A_{\text{row}}(i,j) = a[i\,\text{lda} + j]\).

Inner loop over \(p\):

\[ y_i \mathrel{+}= A_{\text{row}}(i,p)\cdot x_p, \]

so the row of \(A\) is read contiguously (stride 1), but each row of \(A\) is used exactly once (never reused), so this is a pure streaming pass.

2.2 Column-major (`gemv_col`)¶

Element \((i,j)\) is at linear address:

\[ \operatorname{addr}(A_{ij}) = j\cdot\text{lda} + i, \]

matching \(A_{\text{col}}(i,j) = a[j\,\text{lda} + i]\).

Inner loop over \(p\):

\[ y_i \mathrel{+}= A_{\text{col}}(i,p)\cdot x_p, \]

where \(A_{\text{col}}(i,p) = a[p\,\text{lda} + i]\). Varying \(p\) means a stride of \(\text{lda}=m\) — non-contiguous, cache-unfriendly.

3. Memory Access Patterns (naive)¶

Variant	A access (inner)	x access (inner)	y access
`gemv_row`	contiguous (stride 1)	contiguous	single write per \(i\)
`gemv_col`	strided (stride \(m\))	contiguous	single write per \(i\)

In both cases the matrix \(A\) is read once (each element used exactly once for the dot product of its row). There is no reuse of any element of \(A\) across the whole operation — this is intrinsic to matvec and central to its character: unlike matmul, matvec has no data reuse on the matrix side, only on the vector \(x\).

4. Total Arithmetic Work¶

Each \(y_i\) needs \(k\) multiplications and \(k-1\) additions (plus the initial += y_i):

\[ \text{FLOPs per element} = 2k. \]

Total over all \(m\) output elements:

\[ \boxed{\text{FLOPs} = 2\cdot m\cdot k} \]

For \(m=k=2000\):

\[ \text{FLOPs} = 2\times 2000^2 = 8\times10^6 = 8\ \text{MFLOP}. \]

This matches the code's gflops = 2.0 * m*k * 1.0e-09.

Performance metric:

\[ \text{GFLOPS} = \frac{2\cdot m\cdot k \times 10^{-9}}{t_{\text{seconds}}}. \]

5. Roofline: matvec Is Memory-Bound by Construction¶

The total bytes moved (assuming perfect streaming) are:

\[ \text{Traffic} \gtrsim \underbrace{mk\cdot 8}_{\text{A}} \;+\; \underbrace{k\cdot 8}_{\text{x reused}} \;+\; \underbrace{m\cdot 8}_{\text{y}}. \]

The arithmetic intensity is:

\[ I = \frac{2mk\ \text{FLOP}}{(mk + k + m)\cdot 8\ \text{bytes}} \approx \frac{2}{8} = 0.25\ \text{FLOP/byte}. \]

Because each \(A\) element participates in exactly one multiply, the intensity is \(O(1)\) and tiny. matvec cannot be made compute-bound the way matmul can — there is simply not enough arithmetic per byte of the matrix. This is why the optimization ladder here focuses on:

Vectorizing the inner loop (compute several rows at once from the same \(x\) broadcast),
Improving cache/TLB behavior of the matrix reads,
Register blocking of the output accumulators.

rather than on cache blocking for reuse (which is the dominant technique in matmul).

6. Loop Structure¶

Both naive versions use the \(i\)-\(p\) loop order (outer over rows of \(y\), inner over \(p\)):

\[ y_i \mathrel{+}= \sum_{p=0}^{k-1} A_{ip}\,x_p, \qquad i=0,\dots,m-1. \]

Each \(y_i\) is a single dot product of row \(i\) of \(A\) with the full vector \(x\).

7. Summary¶

Aspect	`gemv_row`	`gemv_col`
Layout	row-major	column-major
Inner access (A)	stride 1	stride \(m\)
Inner access (x)	contiguous	contiguous
FLOPs	\(2mk\)	\(2mk\)
Arithmetic intensity	\(O(1)\)	\(O(1)\)
Bottleneck	memory bandwidth	memory bandwidth (worse: strided)

Bottom line: matvec does \(2mk\) FLOPs and is memory-bound by nature — each of the \(mk\) entries of \(A\) is used once. gemv_col is additionally harmed by strided access. The ladder's job is to make the streaming as efficient as possible and vectorize the arithmetic.