optim1
Mathematics of the "Optimized" Column-Major Matrix Multiplication¶
This code implements the same mathematical operation as the previous version, but with a different loop order and column-major memory layout. While the code claims "Optimization 1," let's analyze what actually changes mathematically and why the performance might differ.
1. Core Mathematical Operation¶
The fundamental operation remains identical:
Where:
- \(A \in \mathbb{R}^{m \times k}\)
- \(B \in \mathbb{R}^{k \times n}\)
- \(C \in \mathbb{R}^{m \times n}\)
Each element is computed as:
For \(m = n = k = 2000\), this is still 16 GFLOPs total work.
2. Loop Order Change: \((j, i, p)\) vs \((i, j, p)\)¶
Previous Code (Naive Row-Major)¶
for(int i = 0; i < m; i++) // outer: rows of C
for(int j = 0; j < n; j++) // middle: cols of C
for(int p = 0; p < k; p++) // inner: dot product
C_row(i,j) += A_row(i,p) * B_row(p,j);
Mathematical order: Compute \(C_{00}, C_{01}, \dots, C_{0,n-1}, C_{10}, \dots\)
This Code (Optimization 1)¶
for(int j = 0; j < n; j++) // outer: columns of C
for(int i = 0; i < m; i++) // middle: rows of C
AddDot(k, &A(i,0), lda, &B(0,j), &C(i,j));
Mathematical order: Compute \(C_{00}, C_{10}, \dots, C_{m-1,0}, C_{01}, \dots\)
The total set of computed values is identical, but the access pattern changes.
3. Memory Layout: Column-Major (Fortran-style)¶
This code uses column-major layout, where element \((i,j)\) is stored at:
This is explicit in the macros:
Memory Access Pattern in Inner Loop¶
The AddDot function computes:
Translated to our matrices:
X(p)accessesA(i, p)→ address:p * lda + i→ stride = lda = my[p]accessesB(p, j)→ address:j * ldb + p→ stride = 1 (contiguous!)gammaaccessesC(i, j)→ single location, reused
Stride Analysis¶
| Matrix | Access Pattern | Stride | Cache Behavior |
|---|---|---|---|
| \(A_{ip}\) | Fixed\(i\), varying \(p\) | \(m = 2000\) | Bad (strided read) |
| \(B_{pj}\) | Fixed\(j\), varying \(p\) | \(1\) | Good (contiguous read) |
| \(C_{ij}\) | Fixed\((i,j)\) | - | Good (single write) |
Key insight: In this layout, B is accessed contiguously, but A is strided.
4. Why This Might Be "Optimized"¶
4.1 Reduced Redundancy in C Access¶
In the previous code, C_row(i,j) was accessed in the inner loop, but since \(j\) is the middle loop variable, the same \(C_{ij}\) is not reused across iterations of \(p\).
In this code, C(i,j) is computed once per \((i,j)\) pair, and the entire dot product is performed in a single call to AddDot. This is semantically identical but conceptually cleaner.
4.2 Better Contiguity for B¶
Since \(B\) is stored column-major and we iterate over rows of \(B\) (fixed column \(j\), varying row \(p\)), the access pattern B(p,j) becomes contiguous in memory. This can significantly reduce cache misses compared to the row-major version where \(B\) was strided.
4.3 Single Accumulator Pattern¶
The AddDot function uses a single accumulator (*gamma) that is updated in the inner loop:
Then:
This is the same as the previous code's C_row(i,j) += ... but encapsulated in a function call.
5. Performance Implications¶
5.1 Arithmetic Intensity¶
Still \(O(1)\) FLOPs/byte because:
- Total FLOPs: \(2mnk\)
- Total memory traffic: \(O(mnk)\) due to strided access on A
5.2 Cache Behavior Comparison¶
| Metric | Previous (Row-Major) | This (Column-Major) |
|---|---|---|
| A access | Contiguous in inner loop | Strided (stride = 2000) |
| B access | Strided (stride = 2000) | Contiguous |
| C access | Contiguous | Single write per (i,j) |
| Expected performance | Slower (B strided) | Faster (B contiguous) |
Why? Modern CPUs prefetch contiguous memory much better than strided access. Even though A is now strided, the contiguous B access may dominate performance because B is read \(m\) times per column (once for each row of the result).
5.3 Function Call Overhead¶
The AddDot function adds a function call per \((i,j)\) pair (\(m \times n = 4 \times 10^6\) calls). This overhead is negligible compared to the \(k=2000\) inner loop iterations, but it's not zero.
6. Mathematical Equivalence¶
Both versions compute the exact same result:
The only differences are:
- Loop order: \((i,j,p)\) vs \((j,i,p)\)
- Memory layout: Row-major vs Column-major
- Encapsulation: Direct accumulation vs function call
7. Summary Table¶
| Aspect | Previous (Row-Major) | This (Column-Major) |
|---|---|---|
| Loop order | \((i, j, p)\) | \((j, i, p)\) |
| Memory layout | Row-major | Column-major |
| A access (inner) | Contiguous | Strided (stride = 2000) |
| B access (inner) | Strided (stride = 2000) | Contiguous |
| C access | Contiguous writes | Single write per (i,j) |
| FLOPs | \(2mnk\) | \(2mnk\) |
| Expected speedup | Baseline | Faster (due to B contiguity) |
8. Why It's Not Truly "Optimized"¶
While this version is better than the naive row-major version due to contiguous B access, it's still not optimal:
- Still \(O(mn)\) cache misses on A due to strided access
- No loop blocking/tiling to improve cache reuse
- No vectorization hints (SIMD)
- No parallelization (multi-threading)
A truly optimized version would use blocking/tiled multiplication to keep data in cache, achieving \(O(mn + k(m+n))\) memory traffic instead of \(O(mnk)\).
9. The Real Optimization: Blocked Matrix Multiplication¶
For completeness, here's the mathematical form of a blocked version:
Where \(B\) is the block size (e.g., 64). This keeps blocks of A, B, and C in cache, reducing memory traffic by a factor of \(\approx k/B\).
Conclusion: This "Optimization 1" is a layout and loop-order improvement that makes B access contiguous, which is generally faster on modern hardware. However, it's still a naive \(O(mnk)\) algorithm and doesn't achieve the performance of true blocked/tiled matrix multiplication.
Want me to show the math for a blocked version that achieves near-peak performance?