Complete Mathematics of the Given Matrix Multiplication Code¶

1. Problem Setup¶

Given three matrices:

\[ A \in \mathbb{R}^{m \times k}, \quad B \in \mathbb{R}^{k \times n}, \quad C \in \mathbb{R}^{m \times n} \]

The operation computed is:

\[ C \leftarrow A \cdot B + C \]

For the main() function specifically:

\[ m = n = k = 2000 \]

So \(A, B, C \in \mathbb{R}^{2000 \times 2000}\) (with \(A, B\) being \(2000 \times 2000\) as well since \(m=k=n=2000\)).

2. Element-wise Definition¶

Each element of the product is defined as:

\[ C_{ij} = \sum_{p=0}^{k-1} A_{ip} \cdot B_{pj} + C_{ij}^{(\text{old})} \]

for all \(0 \le i < m\) and \(0 \le j < n\).

This corresponds directly to the inner loop:

C_row(i,j) += A_row(i,p) * B_row(p,j);

3. Row-Major vs Column-Major Memory Layout¶

3.1 Row-Major Layout (used by `matmul_row`)¶

In row-major order, element \(A_{ij}\) is stored at linear address:

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

This matches the macro:

\[ A_{\text{row}}(i,j) = a[i \cdot \text{lda} + j] \]

Memory access pattern in the triple loop \((i, j, p)\):

Loop variable	Access	Stride
\(p\) (innermost)	\(A_{\text{row}}(i,p)\)	1 (contiguous)
\(p\) (innermost)	\(B_{\text{row}}(p,j)\)	\(\text{ldb} = k\) (non-contiguous)
\(j\)	\(C_{\text{row}}(i,j)\)	1 (contiguous)

So in matmul_row, A is accessed contiguously, but B is accessed with stride \(k\) — poor cache locality on B.

3.2 Column-Major Layout (used by `matmul_col`)¶

In column-major order (Fortran/BLAS style), element \(A_{ij}\) is stored at:

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

This matches the macro:

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

Memory access pattern in the triple loop \((i, j, p)\):

Loop variable	Access	Stride
\(p\) (innermost)	\(A_{\text{col}}(i,p)\)	\(\text{lda} = m\) (non-contiguous)
\(p\) (innermost)	\(B_{\text{col}}(p,j)\)	1 (contiguous)
\(j\)	\(C_{\text{col}}(i,j)\)	\(\text{ldc} = m\) (non-contiguous)

So in matmul_col, B is accessed contiguously, but A and C are strided.

4. Total Arithmetic Work¶

Each \(C_{ij}\) requires \(k\) multiplications and \(k\) additions:

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

Total over all \(m \times n\) elements:

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

For \(m = n = k = 2000\):

\[ \text{FLOPs} = 2 \times 2000^3 = 1.6 \times 10^{10} = 16 \text{ GFLOP} \]

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

Performance metric:

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

5. Memory Traffic (Cache Analysis)¶

For a naive triple loop, the theoretical minimum data movement is \(O(mn + mk + nk)\), but the naive implementation incurs much more.

5.1 Naive Row-Major Cost¶

For each of the \(m \cdot n\) output elements, the inner loop over \(p\) reads \(k\) elements of a row of A and \(k\) elements of a column of B.

Reads of A: \(m \cdot n \cdot k\) (each row of A re-read \(n\) times)
Reads of B: \(m \cdot n \cdot k\) (each column of B re-read \(m\) times, with cache misses due to stride \(k\))

Total memory traffic:

\[ \text{Traffic} = O(m \cdot n \cdot k) \]

This is the same order as the arithmetic work, giving an arithmetic intensity of:

\[ I = \frac{2mnk}{O(mnk)} = O(1) \text{ FLOP/byte} \]

This is memory-bound on modern hardware — the CPU spends most of its time waiting for cache misses rather than computing.

5.2 Cache Block Considerations¶

Let \(B\) be the cache block size (in doubles). For matmul_row, accessing \(B_{\text{row}}(p,j)\) with stride \(k\) pulls in a full cache line for each element used:

\[ \text{Effective traffic per element} \approx \frac{\text{cache line size}}{\text{sizeof(double)}} = \frac{64}{8} = 8\times \text{ over-fetch} \]

6. Loop Nest and Accumulator Form¶

The code's loop order is \((i, j, p)\):

\[ C_{ij} \mathrel{+}= \sum_{p=0}^{k-1} A_{ip} B_{pj} \]

This is the "dot product" formulation, where each output element is computed independently as a single dot product of a row of A with a column of B.

An equivalent but more cache-friendly axpy / rank-1 update form (loop order \(i, p, j\)) would be:

\[ C_{ij} \mathrel{+}= A_{ip} \cdot B_{pj} \quad \text{for all } j \]

which makes the innermost loop over \(j\) contiguous in row-major layout. The given code does not use this, which is why both versions perform poorly.

7. Summary Table¶

Aspect	`matmul_row`	`matmul_col`
Layout	Row-major	Column-major
Inner loop stride (A)	1 (good)	\(m\) (bad)
Inner loop stride (B)	\(k\) (bad)	1 (good)
Inner loop stride (C)	1 (good)	\(m\) (bad)
FLOPs	\(2mnk\)	\(2mnk\)
Arithmetic intensity	\(O(1)\)	\(O(1)\)
Bottleneck	Memory bandwidth / cache misses	Memory bandwidth / cache misses

Bottom line: both versions do the same \(2mnk\) arithmetic work, but the naive \((i,j,p)\) loop order makes them memory-bound, with the row-major version generally being faster on x86 because the contiguous A access plus contiguous C write partially offsets the strided B reads.

Want me to show the math for a blocked/tiled version that turns this into a compute-bound kernel?