Matrix Multiplication Optimization in C

From a naive triple loop to a cache-blocked, vectorized implementation, one change at a time, with measurements at each step.

Course project for Introduction to Parallel and Distributed Computing, IME-USP (2024).

The goal was not to beat BLAS. It was to understand, with numbers, how much of matrix multiplication performance comes from memory access patterns before any parallelism is involved. The answer: almost all of it.

Step	Build	Technique	Time (s)	Speedup vs. baseline
1	exe_1	naive i-j-k, no flags		1.0x
2	exe_2	row-major access
3	exe_3	+ -O3
4	exe_4	+ loop reordering (i-k-j)
5	exe_5	+ auto-vectorization (-march=native)
6	exe_6	+ cache blocking (tiling)
7	exe_7	+ recursive block decomposition

Final result: 99.3% less execution time than the baseline on the same input size and machine.

---

The optimization path

Each step keeps the previous ones.

1. Row-major access. C stores arrays row by row. Traversing B column-wise touches a new cache line on every step. Reorganizing the access so inner loops walk rows gives spatial locality for free.

2. Compiler flags (-O3). Lets the compiler unroll, inline and schedule instructions. Cheap win, but it does not fix a bad access pattern.

3. Loop reordering (i-k-j). With i-k-j, the innermost loop updates a contiguous row of C using a contiguous row of B and a single scalar A[i][k]. The hardware prefetcher can keep up and cache misses drop sharply. This is usually the single largest gain.

4. Vectorization (SIMD). With contiguous inner loops, -march=native -ftree-vectorize lets GCC emit vector instructions (AVX2 or AVX-512 depending on the CPU) that process several elements per instruction.

5. Cache blocking (tiling). For large matrices even a row does not fit in L1. Splitting the matrices into sub-blocks sized to the cache keeps each block resident while it is reused, instead of streaming it from memory repeatedly.

6. Recursive decomposition. Dividing the matrices into quadrants recursively gives cache-oblivious behavior: the recursion naturally reaches a size that fits whichever cache level is available.

---

Profiling

Each version was profiled to confirm where the time went, not just that it dropped:

• gprof for the call graph and hot functions.
• perf stat for cycles, instructions per cycle (IPC), L1 and LLC cache misses.
• perf record / perf report for the hot loop.

# example
perf stat -e cycles,instructions,cache-misses,L1-dcache-load-misses ./bin/exe_6

The interesting comparison is between steps 3 and 4: instruction count barely changes, but cache misses collapse and IPC roughly doubles.

---

Build and run

# all versions
make

# single version
make exe_5

# run
./bin/exe_5 <matrix_size>

Or by hand:

gcc -o bin/exe_1 src/mat_mul_1.c
gcc -o bin/exe_2 src/mat_mul_2.c
gcc -O3 -o bin/exe_3 src/mat_mul_2.c
gcc -O3 -o bin/exe_4 src/mat_mul_3.c
gcc -O3 -march=native -ftree-vectorize -o bin/exe_5 src/mat_mul_3.c
gcc -O3 -march=native -ftree-vectorize -o bin/exe_6 src/mat_mul_4.c
gcc -O3 -march=native -ftree-vectorize -o bin/exe_7 src/mat_mul_5.c

---

Repository structure

matrix_multiply_optimizer/
├── src/            # one file per optimization stage (mat_mul_1.c ... mat_mul_5.c)
├── bin/            # compiled executables
├── benchmark/      # timing scripts and raw results
├── benchmark.png   # summary chart
└── README.md

---

What carried over into production work

The lesson generalizes beyond matrices: in data systems, layout and access pattern decide performance long before you add threads or nodes. The same reasoning drives partition design in Cassandra, file layout in Iceberg and how you order columns in a Parquet scan.

---

License

GPL-3.0. See LICENSE.