Optimize a Matrix Multiplication Kernel

Prove 5-30x speedups on a core AI workload using cache-aware C++.

Introduction

30 Second Summary

Every time you ask an AI model a question, billions of matrix multiplications happen under the hood. The speed of those calculations determines whether inference takes milliseconds or minutes.

In this project, you will write a naive matrix multiplication in C++, then apply cache-aware optimizations and compiler flags to achieve measurable 5-30x speedups. You will build an automated Python benchmarking tool to prove your results with hard numbers.

What You'll Build

A complete performance engineering mini-lab where you optimize the core computation behind AI workloads and measure every improvement with automated tooling.

By the end of this project, you'll have:

  • A working C++ matrix multiplication program with naive and cache-optimized variants.
  • A Python benchmarking harness that automates compilation, timing, and visualization of all optimization strategies.
  • A system architecture diagram mapping each optimization technique to its target in the CPU memory hierarchy.
  • Secret Mission: Implement cache blocking (tiling) to squeeze even more performance by fitting sub-matrices into L1/L2 cache.

Do I need to know C++ or Linux?

No prior C++ or Linux experience is required. You only need a basic understanding of loops and arrays in any programming language.

The project guides you through every command and line of code from scratch.

Set Up WSL, GCC, and Python

This project uses Linux-based tools to compile and benchmark C++ code. WSL (Windows Subsystem for Linux) gives you a full Ubuntu environment inside Windows, so you can run GCC, Python, and standard Linux utilities without leaving your machine.

In this step, you will install WSL 2 with Ubuntu, add the C++ compiler and Python plotting library, and create your project directory.

In this step, get ready to:
  • Install WSL 2 with Ubuntu on your Windows machine.
  • Install g++ and python3-matplotlib inside your Ubuntu terminal.
  • Create the project directory and verify all tools are working.
Install WSL 2 with Ubuntu

WSL 2 runs a real Linux kernel on your Windows machine. It lets you use Linux tools like g++ and python3 natively, which is exactly what performance engineering work requires.

  • Press the Windows key to open your search bar.
  • Type PowerShell, then right-click Windows PowerShell and select Run as administrator.
  • Check whether you already have Ubuntu on WSL by running this command:
wsl -l -v

✔️ I see Ubuntu in the list

You already have Ubuntu installed on WSL. Confirm the VERSION column shows 2. If it does, you can skip ahead to opening your Ubuntu terminal.

  • Type Ubuntu in PowerShell and press Enter to open your Ubuntu terminal.

ⓧ I see the WSL help text or an error

WSL is not installed yet (or has no distributions). The wsl --install command handles the full setup in one step.

  • Install WSL 2 with Ubuntu by running this command in the same administrator PowerShell:
wsl --install

What does this command do?

This enables the WSL feature, downloads the Linux kernel, and installs Ubuntu as your default distribution. New installations via wsl --install default to WSL 2.

  • Wait for the installation to complete. You will see progress messages in the terminal.
  • Restart your computer when prompted. This is required for WSL to finish setup.
  • After restarting, Ubuntu will launch automatically (or press the Windows key, type Ubuntu, and press Enter).
  • Ubuntu will ask you to create a Linux username. Type a short username (all lowercase, no spaces) and press Enter.
  • Type a password when prompted. The characters will not appear as you type. Press Enter, then retype the same password to confirm.

Password is invisible

Linux terminals never display password characters. This is normal. Just type your password and press Enter. You will need this password later when running sudo commands.

You should now see a Ubuntu terminal prompt that looks like yourname@MACHINE:~$.

Having trouble installing WSL?

Make sure you opened PowerShell as an administrator. The wsl --install command requires elevated permissions.

If you see the help text instead of an installation starting, WSL may already be partially installed. Try running wsl --install -d Ubuntu to install the Ubuntu distribution specifically.

Help me troubleshoot my WSL installation

Install the C++ compiler and Python plotting library

Ubuntu comes with Python 3 pre-installed, but you still need the C++ compiler and the matplotlib plotting library for your benchmarking harness.

  • In your Ubuntu terminal, install g++ and python3-matplotlib by running this command:
sudo apt update && sudo apt install g++ python3-matplotlib

What does this command do?

The first half (sudo apt update) refreshes Ubuntu's package list so it knows where to find the latest versions. The second half (sudo apt install) downloads and installs both packages.

The g++ package is the GNU C++ compiler. The python3-matplotlib package is the charting library your benchmarking script will use to generate PNG charts.

  • When prompted for your password, type the Linux password you created during Ubuntu setup and press Enter.
  • When prompted to confirm the installation (Do you want to continue? [Y/n]), type Y and press Enter.
  • Wait for the installation to finish. You will see package download progress in the terminal.

Seeing a 'Unable to locate package' error?

Make sure the sudo apt update part ran successfully before the install. If it failed, your package list may be out of date. Run sudo apt update on its own first, then retry the install.

Help me fix my apt install error

Create your project directory and verify tools

Your project needs a dedicated folder to keep all the C++ source files, the Python harness, and the benchmark results organized.

  • Create the project directory and move into it by running these commands:
mkdir ~/matmul-optimization && cd ~/matmul-optimization

What does this command do?

The mkdir command creates a new folder called matmul-optimization in your home directory (~). The cd command moves you into that folder so all subsequent commands run from there.

  • Verify that g++ is installed by running:
g++ --version
  • You should see version information printed to the terminal. The exact version number will depend on your Ubuntu release.
  • Verify that Python 3 is installed by running:
python3 --version
  • You should see a version number like Python 3.x.x printed to the terminal.
  • Confirm you are in the correct directory by running:
pwd
  • You should see /home/yourname/matmul-optimization (where yourname is the Linux username you created earlier).

g++ or python3 not found?

If either command says it is not found, the apt install step may not have completed successfully. Re-run the install command:

sudo apt install g++ python3-matplotlib

Help me troubleshoot my g++ or Python installation

Your development environment is ready. Next up, you will write your first matrix multiplication program in C++ and establish a performance baseline.

Write and Run the Naive Matrix Multiply

Your WSL environment is set up with g++ and Python ready to go. Now it's time to write the code that forms the foundation of your entire benchmarking lab.

Matrix multiplication is the core operation behind every neural network. Before you can optimize it, you need a baseline implementation to measure against.

In this step, get ready to:
  • Write a naive matrix multiplication program in C++.
  • Compile and run it to record your baseline performance.
Create the naive matrix multiplication source file

Your naive implementation uses the most intuitive loop order (ijk) for multiplying two matrices. It stores matrices as flat 1D arrays in row-major order, where each row sits contiguously in memory.

  • Create and open a new file by running this command:
nano matmul_naive.cpp

This opens the nano text editor with a new file called matmul_naive.cpp.

  • Copy the following code into the editor:
#include <iostream>
#include <chrono>
#include <cstdlib>

int main(int argc, char* argv[]) {
    // Default matrix size, overridable from command line
    int N = 512;
    if (argc > 1) N = std::atoi(argv[1]);

    // Allocate three NxN matrices as 1D arrays (row-major layout)
    double* A = new double[N * N];
    double* B = new double[N * N];
    double* C = new double[N * N];

    // Fill A and B with random values, zero out C
    srand(42);
    for (int i = 0; i < N * N; i++) {
        A[i] = (double)rand() / RAND_MAX;
        B[i] = (double)rand() / RAND_MAX;
        C[i] = 0.0;
    }

    // Start timing the multiplication
    auto start = std::chrono::high_resolution_clock::now();

    // Naive ijk loop order: C[i][j] += A[i][k] * B[k][j]
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                C[i * N + j] += A[i * N + k] * B[k * N + j];
            }
        }
    }

    // Stop timing
    auto end = std::chrono::high_resolution_clock::now();
    double elapsed = std::chrono::duration<double, std::milli>(end - start).count();

    // Print elapsed time to stdout (used by benchmarking harness later)
    std::cout << elapsed << std::endl;

    // Print checksum to stderr for verification
    double checksum = 0.0;
    for (int i = 0; i < N * N; i++) checksum += C[i];
    std::cerr << "Checksum: " << checksum << std::endl;

    // Free allocated memory
    delete[] A;
    delete[] B;
    delete[] C;
    return 0;
}

What does this code do?

This program multiplies two 512x512 matrices and times how long the multiplication takes.

  • N = 512 sets the matrix dimensions. You can override this from the command line later.
  • The three arrays (A, B, C) are stored as flat 1D arrays. A 2D position [i][j] maps to index i * N + j in the 1D array.
  • The ijk loop order means the innermost loop (k) strides through column positions of B. Each access to B[k * N + j] jumps N elements forward in memory, which is terrible for cache locality.
  • srand(42) seeds the random number generator with a fixed value so every run produces identical matrices. This makes your benchmarks reproducible.
  • The checksum printed to stderr lets you verify that different optimization variants produce the same result.
  • Save the file by pressing Ctrl+O, then press Enter to confirm the filename.
  • Exit nano by pressing Ctrl+X.
  • Confirm the file was saved by running:
ls

You should see matmul_naive.cpp listed in the output.

Don't see the file?

  • Make sure you pressed Ctrl+O (the letter O, not zero) and then Enter before exiting. Nano shows a confirmation at the bottom of the screen.
  • Verify you're in the right directory by running pwd. It should show /home/your-username/matmul-optimization.

Help me figure out why my file didn't save in nano.

Compile and run your baseline

Now you'll compile the source into an executable binary with all optimizations disabled. The -O0 flag tells GCC to skip all optimizations so you can see the raw algorithmic performance.

  • Compile the program with no optimization by running:
g++ matmul_naive.cpp -o matmul_naive -O0

What does this command do?

  • g++ invokes the GNU C++ compiler.
  • matmul_naive.cpp is your source file.
  • -o matmul_naive names the output binary.
  • -O0 disables optimization. GCC prioritizes fast compilation and straightforward debugging at this level.
  • Run the binary by executing:
./matmul_naive

The program will pause for several seconds while it multiplies the matrices. You should see a number printed to your terminal (the elapsed time in milliseconds) followed by a checksum line.

What should I see?

Your output will look something like this (your exact number will differ based on your hardware):

1842.37
Checksum: 33554816.123456

The first line is the elapsed time in milliseconds. On most machines, the naive ijk implementation takes several hundred milliseconds to a few seconds for a 512x512 matrix. This number is your baseline.

The checksum line goes to stderr and confirms the computation produced a result. You'll use this later to verify that optimized versions produce the same answer.

Seeing a compilation error?

  • Check for typos in the code. Common mistakes include missing semicolons, mismatched braces, or using std:chrono instead of std::chrono (double colon).
  • If you see command not found, verify g++ is installed by running g++ --version.
  • You can re-open the file to fix errors by running nano matmul_naive.cpp again.

Help me fix my g++ compilation error for matmul_naive.cpp.

✔️ Awesome, I've got everything!

Great. Double check you've saved your file and note down the elapsed time. That number is the one you'll beat in the next steps.

ⓧ I'd like to double check the full code

Here is the complete matmul_naive.cpp file. Compare your code character by character:

#include <iostream>
#include <chrono>
#include <cstdlib>

int main(int argc, char* argv[]) {
    int N = 512;
    if (argc > 1) N = std::atoi(argv[1]);

    double* A = new double[N * N];
    double* B = new double[N * N];
    double* C = new double[N * N];

    srand(42);
    for (int i = 0; i < N * N; i++) {
        A[i] = (double)rand() / RAND_MAX;
        B[i] = (double)rand() / RAND_MAX;
        C[i] = 0.0;
    }

    auto start = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                C[i * N + j] += A[i * N + k] * B[k * N + j];
            }
        }
    }

    auto end = std::chrono::high_resolution_clock::now();
    double elapsed = std::chrono::duration<double, std::milli>(end - start).count();

    std::cout << elapsed << std::endl;

    double checksum = 0.0;
    for (int i = 0; i < N * N; i++) checksum += C[i];
    std::cerr << "Checksum: " << checksum << std::endl;

    delete[] A;
    delete[] B;
    delete[] C;
    return 0;
}

You now have a concrete number to beat. In the next step, you'll document the system architecture to understand exactly where in the memory hierarchy this code is bottlenecked.

Create the System Architecture Diagram

You have a working naive matrix multiply with a baseline timing number. Before you start optimizing, you need a map of the system you are working with.

In this step, you will create three architecture diagrams that show how your benchmarking pipeline works, how your CPU's memory hierarchy is structured, and which optimization technique targets which level of that hierarchy. These diagrams make the "why" behind each optimization concrete and visual.

In this step, get ready to:
  • Create a benchmarking pipeline diagram in the Mermaid Live Editor.
  • Create a CPU memory hierarchy diagram with latencies and capacities.
  • Create an optimization targets diagram and save all sources as architecture.md.
Create the benchmarking pipeline diagram

The Mermaid Live Editor is a free web tool that renders diagrams from plain text. You write code that describes nodes and arrows, and the editor renders a visual flowchart instantly.

  • Open your web browser (on Windows, outside WSL) and navigate to https://mermaid.live.
  • You will see a code editor on the left and a diagram preview on the right. Select all the default code in the left panel and delete it.
  • Paste the following Mermaid code into the editor:
flowchart LR
    A[benchmark.py\nPython Harness] -->|compiles with| B[g++\nGCC Compiler]
    B -->|produces| C[C++ Binaries\nnaive + optimized]
    C -->|runs and outputs| D[Timing Data\nstdout: elapsed ms]
    D -->|collected by| A
    A -->|generates| E[Results\nspeedup chart + table]

What does this diagram show?

This is the pipeline you will build in Step 5. The Python harness (benchmark.py) compiles your C++ source files using g++, runs the resulting binaries, collects timing from stdout, and generates a chart comparing all variants.

The \n in the node labels creates line breaks inside each box. The -->|text| syntax adds labels to the arrows between nodes.

The preview panel on the right should now show a left-to-right flowchart with five connected boxes: benchmark.py, g++, C++ Binaries, Timing Data, and Results.

Don't see a rendered diagram?

Make sure you deleted all the default code before pasting. Any leftover text can cause a syntax error. Also check that the first line says flowchart LR with no extra spaces before it.

Why isn't my Mermaid diagram rendering in the live editor?

Create the CPU memory hierarchy diagram

Your next diagram visualizes why cache locality matters so much. A CPU does not access RAM directly for every operation. It uses a hierarchy of progressively larger but slower caches.

  • Select all the code in the Mermaid editor's left panel and delete it.
  • Paste the following code to create a top-to-bottom hierarchy diagram:
flowchart TB
    R[CPU Registers\n~1 cycle latency\n~1 KB capacity]
    L1[L1 Cache\n~4 cycles latency\n32-64 KB capacity]
    L2[L2 Cache\n~12 cycles latency\n256 KB to 1 MB capacity]
    L3[L3 Cache\n~30 cycles latency\n8-32 MB capacity]
    RAM[Main Memory RAM\n~100+ cycles latency\n8-64 GB capacity]

    R --- L1 --- L2 --- L3 --- RAM

What does this diagram show?

Each tier is labeled with its approximate latency (in CPU cycles) and capacity. Registers are fastest but tiny. RAM is huge but over 100x slower than registers.

The --- connector draws lines without arrowheads, showing the tiers as a connected chain rather than a directed flow. TB means "top to bottom" layout.

You should see a vertical chain of five boxes, from CPU Registers at the top down to Main Memory RAM at the bottom.

Why does this matter for matrix multiplication?

Your 512x512 matrices use double-precision floats. Each matrix is 512 * 512 * 8 bytes = 2 MB. Three matrices total 6 MB. That exceeds L1 cache (32-64 KB) and L2 cache (256 KB - 1 MB) on most CPUs.

When your loop accesses memory in a pattern that jumps around (like the naive ijk order accessing columns of B), the CPU constantly fetches data from slow RAM instead of fast cache. Loop reordering fixes this by accessing memory sequentially.

Create the optimization targets diagram and save all sources

Your final diagram maps each optimization technique you will apply in this project to the specific level of the memory hierarchy it targets.

  • Select all the code in the Mermaid editor and delete it.
  • Paste the following code to create a diagram with two subgraphs:
flowchart LR
    subgraph Optimizations
        O1[Loop Reordering\nijk to ikj]
        O2[Compiler Flags\n-O2 -O3 -march=native]
        O3[Cache Blocking\ntiling into sub-matrices]
    end

    subgraph Targets
        T1[Reduce L1 and L2\ncache misses]
        T2[Better register use\nand instruction pipeline]
        T3[Maximize L2 and L3\ndata reuse]
    end

    O1 --> T1
    O2 --> T2
    O3 --> T3

What does this diagram show?

The left subgraph lists the three optimization strategies: loop reordering (Step 4), compiler flags (Step 5), and cache blocking (Secret Mission). The right subgraph shows what each targets in the hardware.

Loop reordering reduces L1/L2 cache misses by accessing memory sequentially. Compiler flags like -O3 and -march=native let the compiler optimize register usage and instruction scheduling. Cache blocking fits sub-matrices into L2/L3 cache for maximum reuse.

You should see two grouped boxes (Optimizations on the left, Targets on the right) with arrows connecting each optimization to its target.

  • Export each diagram as a PNG by clicking the **Actions** menu above the preview panel, then selecting **Download PNG**. Repeat this for all three diagrams (paste each one back in to export it).

Now save all three Mermaid sources together in a single markdown file in your project directory.

  • In your Ubuntu terminal, create architecture.md by running:
nano architecture.md
  • Paste the following content into nano. This combines all three diagram sources into one reference file:
# System Architecture: Matrix Multiplication Optimization

## Benchmarking Pipeline

```mermaid
flowchart LR
    A[benchmark.py\nPython Harness] -->|compiles with| B[g++\nGCC Compiler]
    B -->|produces| C[C++ Binaries\nnaive + optimized]
    C -->|runs and outputs| D[Timing Data\nstdout: elapsed ms]
    D -->|collected by| A
    A -->|generates| E[Results\nspeedup chart + table]
```

## CPU Memory Hierarchy

```mermaid
flowchart TB
    R[CPU Registers\n~1 cycle latency\n~1 KB capacity]
    L1[L1 Cache\n~4 cycles latency\n32-64 KB capacity]
    L2[L2 Cache\n~12 cycles latency\n256 KB to 1 MB capacity]
    L3[L3 Cache\n~30 cycles latency\n8-32 MB capacity]
    RAM[Main Memory RAM\n~100+ cycles latency\n8-64 GB capacity]

    R --- L1 --- L2 --- L3 --- RAM
```

## Optimization Targets

```mermaid
flowchart LR
    subgraph Optimizations
        O1[Loop Reordering\nijk to ikj]
        O2[Compiler Flags\n-O2 -O3 -march=native]
        O3[Cache Blocking\ntiling into sub-matrices]
    end

    subgraph Targets
        T1[Reduce L1 and L2\ncache misses]
        T2[Better register use\nand instruction pipeline]
        T3[Maximize L2 and L3\ndata reuse]
    end

    O1 --> T1
    O2 --> T2
    O3 --> T3
```

What is this file for?

This markdown file serves as living documentation for your project. Any tool that supports Mermaid rendering (GitHub, GitLab, VS Code) will display these as visual diagrams automatically. The file also acts as a quick reference for how your pipeline fits together.

  • Save and exit nano by pressing **Ctrl+O**, then **Enter** to confirm the filename, then **Ctrl+X** to exit.
  • Verify the file was saved by running:
ls

You should see architecture.md alongside matmul_naive.cpp and matmul_naive in the output.

Don't see architecture.md in the output?

Make sure you pressed **Ctrl+O** (the letter O, not zero) to write the file before exiting nano. If you pressed **Ctrl+X** without saving, nano will ask if you want to save. Press **Y** then **Enter**.

Also confirm you are in the ~/matmul-optimization directory by running pwd.

Help me troubleshoot saving a file with nano in WSL.

✔️ Awesome, I've got everything!

Great. Double check you have saved architecture.md and exported your PNG diagrams from the Mermaid Live Editor.

ⓧ I'd like to double check the full code

Here is the complete architecture.md file. Compare it against what you have saved:

# System Architecture: Matrix Multiplication Optimization

## Benchmarking Pipeline

```mermaid
flowchart LR
    A[benchmark.py\nPython Harness] -->|compiles with| B[g++\nGCC Compiler]
    B -->|produces| C[C++ Binaries\nnaive + optimized]
    C -->|runs and outputs| D[Timing Data\nstdout: elapsed ms]
    D -->|collected by| A
    A -->|generates| E[Results\nspeedup chart + table]
```

## CPU Memory Hierarchy

```mermaid
flowchart TB
    R[CPU Registers\n~1 cycle latency\n~1 KB capacity]
    L1[L1 Cache\n~4 cycles latency\n32-64 KB capacity]
    L2[L2 Cache\n~12 cycles latency\n256 KB to 1 MB capacity]
    L3[L3 Cache\n~30 cycles latency\n8-32 MB capacity]
    RAM[Main Memory RAM\n~100+ cycles latency\n8-64 GB capacity]

    R --- L1 --- L2 --- L3 --- RAM
```

## Optimization Targets

```mermaid
flowchart LR
    subgraph Optimizations
        O1[Loop Reordering\nijk to ikj]
        O2[Compiler Flags\n-O2 -O3 -march=native]
        O3[Cache Blocking\ntiling into sub-matrices]
    end

    subgraph Targets
        T1[Reduce L1 and L2\ncache misses]
        T2[Better register use\nand instruction pipeline]
        T3[Maximize L2 and L3\ndata reuse]
    end

    O1 --> T1
    O2 --> T2
    O3 --> T3
```

Your architecture documentation is in place. You now have a visual map of the entire system: the benchmarking pipeline, the hardware constraints, and exactly where each optimization fits. Next up, you will write the cache-optimized matrix multiply and see these diagrams come to life as measurable speedups.

Write the Cache-Optimized Matrix Multiply

In the previous step, you mapped out the CPU memory hierarchy and identified where each optimization targets. Now it's time to put that theory into practice.

Your naive multiplication uses ijk loop order, which causes the inner loop to jump across rows of matrix B. Every jump means a cache miss. In this step, you will reorder the loops to ikj, which makes both B and C arrays flow sequentially through memory. This single change eliminates most L1 cache misses and delivers a dramatic speedup.

In this step, get ready to:
  • Create a cache-optimized matrix multiplication with ikj loop order.
  • Compile and run the optimized version.
  • Compare timing against your naive baseline to measure the speedup.
Create the optimized source file

The key insight is about row-major memory layout. C++ stores 2D arrays row by row in a flat memory block. When your inner loop walks along a row (incrementing j), it reads consecutive memory addresses. The CPU prefetcher loves this pattern and loads entire cache lines ahead of time.

Why does loop order matter so much?

In your naive ijk code, the inner k-loop accesses B[k * N + j] where k changes each iteration. That jumps 512 doubles (4 KB) between reads, missing the cache every single time.

In ikj order, the inner j-loop accesses B[k * N + j] where j changes. That walks sequentially through memory, hitting the cache on almost every read. The CPU fetches a 64-byte cache line and uses all 8 doubles in it before moving on.

  • Create and open matmul_optimized.cpp by running:
nano matmul_optimized.cpp
  • Type the following code into the file:
#include <iostream>
#include <chrono>
#include <cstdlib>

int main(int argc, char* argv[]) {
    int N = 512;
    if (argc > 1) N = std::atoi(argv[1]);

    double* A = new double[N * N];
    double* B = new double[N * N];
    double* C = new double[N * N];

    srand(42);
    for (int i = 0; i < N * N; i++) {
        A[i] = (double)rand() / RAND_MAX;
        B[i] = (double)rand() / RAND_MAX;
        C[i] = 0.0;
    }

    auto start = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < N; i++) {
        for (int k = 0; k < N; k++) {
            double a_ik = A[i * N + k];
            for (int j = 0; j < N; j++) {
                C[i * N + j] += a_ik * B[k * N + j];
            }
        }
    }

    auto end = std::chrono::high_resolution_clock::now();
    double elapsed = std::chrono::duration<double, std::milli>(end - start).count();

    std::cout << elapsed << std::endl;

    double checksum = 0.0;
    for (int i = 0; i < N * N; i++) checksum += C[i];
    std::cerr << "Checksum: " << checksum << std::endl;

    delete[] A;
    delete[] B;
    delete[] C;
    return 0;
}

What changed from the naive version?

The only difference is the loop section. Everything else (includes, array setup, timing, output, cleanup) is identical to matmul_naive.cpp. Here is what changed:

  • The middle loop now iterates over k instead of j. The inner loop iterates over j instead of k.
  • double a_ik = A[i * N + k] caches the value of A at position (i, k) in a local variable. Since i and k are fixed in the inner loop, this avoids re-reading A from memory on every iteration.
  • The inner j-loop now accesses both B[k * N + j] and C[i * N + j] sequentially. Both arrays flow through cache lines in order.
  • Save the file by pressing Ctrl+O, then Enter to confirm the filename.
  • Exit nano by pressing Ctrl+X.
  • Confirm the file exists by running:
ls

You should see matmul_optimized.cpp listed alongside matmul_naive.cpp and architecture.md.

Don't see matmul_optimized.cpp?

Make sure you pressed Ctrl+O (the letter O, not zero) to write the file before exiting nano. If the file is missing, run nano matmul_optimized.cpp again and re-enter the code.

Also check you are in the ~/matmul-optimization directory. Run pwd to confirm your current location.

Still stuck? Help me figure out why my file wasn't saved in nano on WSL Ubuntu.

Compile and run the optimized version

You will compile with the same -O0 flag (no compiler optimization) so the speedup you see comes purely from your loop reordering, not from the compiler doing extra work.

  • Compile the optimized version by running:
g++ matmul_optimized.cpp -o matmul_optimized -O0
  • Run the optimized binary by running:
./matmul_optimized

You should see a time in milliseconds printed to your terminal. This number should be significantly lower than your naive baseline from Step 2.

What should I see?

If your naive version took around 800-2000 ms, the optimized version should finish in roughly 80-300 ms. That is a 5-12x speedup from changing nothing but the loop order. The exact numbers depend on your CPU and its cache sizes.

The checksum printed to stderr should match the naive version exactly. Both programs compute the same matrix product. If the checksums differ, there is a typo in the loop logic.

  • Run the naive version again to get a fresh comparison:
./matmul_naive

Note both timing values. Divide the naive time by the optimized time to get your speedup ratio.

Seeing a compilation error?

The most common issue is a typo in the loop variables. Double check that your loops are ordered as i (outer), k (middle), j (inner). Make sure the line double a_ik = A[i * N + k]; sits inside the k-loop but before the j-loop.

If you see a "No such file or directory" error, confirm you are in ~/matmul-optimization by running pwd.

Need more help? Help me debug my matmul_optimized.cpp compilation error.

✔️ Awesome, I've got everything!

Great. Double check you have saved matmul_optimized.cpp and that both binaries run without errors.

ⓧ I'd like to double check the full code

Here is the complete matmul_optimized.cpp file:

#include <iostream>
#include <chrono>
#include <cstdlib>

int main(int argc, char* argv[]) {
    int N = 512;
    if (argc > 1) N = std::atoi(argv[1]);

    double* A = new double[N * N];
    double* B = new double[N * N];
    double* C = new double[N * N];

    srand(42);
    for (int i = 0; i < N * N; i++) {
        A[i] = (double)rand() / RAND_MAX;
        B[i] = (double)rand() / RAND_MAX;
        C[i] = 0.0;
    }

    auto start = std::chrono::high_resolution_clock::now();

    for (int i = 0; i < N; i++) {
        for (int k = 0; k < N; k++) {
            double a_ik = A[i * N + k];
            for (int j = 0; j < N; j++) {
                C[i * N + j] += a_ik * B[k * N + j];
            }
        }
    }

    auto end = std::chrono::high_resolution_clock::now();
    double elapsed = std::chrono::duration<double, std::milli>(end - start).count();

    std::cout << elapsed << std::endl;

    double checksum = 0.0;
    for (int i = 0; i < N * N; i++) checksum += C[i];
    std::cerr << "Checksum: " << checksum << std::endl;

    delete[] A;
    delete[] B;
    delete[] C;
    return 0;
}

You just proved that understanding your hardware's memory hierarchy can deliver a 5-12x speedup without changing the algorithm itself. Next up, you will build a Python benchmarking harness that automates compilation, runs all variants with different compiler flags, and charts the results.

Build and Run the Python Benchmarking Harness

You now have two C++ programs that multiply matrices at very different speeds. You proved the speedup manually by compiling and running each one yourself.

Manual testing does not scale. To measure performance rigorously, you need multiple compilation strategies, multiple runs per strategy, and statistical analysis of the results. In this step, you will build a Python script that automates the entire pipeline: compiling variants, running benchmarks, computing statistics, and generating a chart.

In this step, get ready to:
  • Write a Python benchmarking harness that compiles and runs multiple C++ variants.
  • Add chart generation using matplotlib to visualize speedups.
  • Execute the harness and observe 15-30x combined speedup.
Write the benchmarking harness core functions

The harness needs three utility functions: one to compile a C++ source file with specified flags, one to run a binary multiple times and collect timing, and one to format the results.

  • Create benchmark.py by running:
nano benchmark.py
  • Type the following imports and utility functions at the top of the file:
import subprocess
import sys
import statistics


def compile_program(source, output, flags=None):
    """Compile a C++ source file with specified flags."""
    cmd = ["g++", source, "-o", output]
    if flags:
        cmd.extend(flags)
    result = subprocess.run(cmd, capture_output=True, text=True)
    if result.returncode != 0:
        print(f"Compilation failed for {source}:")
        print(result.stderr)
        sys.exit(1)
    print(f"  Compiled: {output} (flags: {' '.join(flags) if flags else 'none'})")


def run_benchmark(executable, matrix_size, runs=5):
    """Run a benchmark executable multiple times and return timing results."""
    times = []
    for _ in range(runs):
        result = subprocess.run(
            [f"./{executable}", str(matrix_size)],
            capture_output=True, text=True
        )
        if result.returncode != 0:
            print(f"Error running {executable}: {result.stderr}")
            sys.exit(1)
        elapsed = float(result.stdout.strip())
        times.append(elapsed)
    return times


def print_results(label, times):
    """Print formatted benchmark results."""
    avg = statistics.mean(times)
    std = statistics.stdev(times) if len(times) > 1 else 0.0
    print(f"  {label:.<45} {avg:>8.1f} ms (std: {std:.1f} ms)")
    return avg

What do these functions do?

  • compile_program() builds a C++ source file into a binary by calling g++ as a subprocess. It exits immediately if compilation fails.
  • run_benchmark() runs a compiled binary multiple times, reads the elapsed milliseconds from stdout each time, and returns the list of timings.
  • print_results() computes the mean and standard deviation of a timing list, prints a formatted line, and returns the average.

You won't see a change yet. The next section adds the chart generation and the main function that ties everything together.

Add chart generation and the main function

The matplotlib library can render charts to a PNG file without needing a display server. This is critical inside WSL where no GUI is available.

  • Continue typing in the same file, directly below the print_results function. Add the chart generation function:

def generate_chart(results, matrix_size):
    """Generate a bar chart of benchmark results."""
    try:
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt
    except ImportError:
        print("\n  matplotlib not available, skipping chart.")
        return

    labels = list(results.keys())
    values = list(results.values())
    baseline = values[0]
    speedups = [baseline / v for v in values]

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

    colors = ["#e74c3c", "#3498db", "#f39c12", "#2ecc71", "#9b59b6"]

    bars1 = ax1.bar(range(len(labels)), values, color=colors)
    ax1.set_ylabel("Time (ms)")
    ax1.set_title(f"Matrix Multiplication Time (N={matrix_size})")
    ax1.set_xticks(range(len(labels)))
    ax1.set_xticklabels(labels, rotation=25, ha="right", fontsize=8)
    for bar, val in zip(bars1, values):
        ax1.text(
            bar.get_x() + bar.get_width() / 2, bar.get_height() + 10,
            f"{val:.0f}", ha="center", va="bottom", fontsize=9
        )

    bars2 = ax2.bar(range(len(labels)), speedups, color=colors)
    ax2.set_ylabel("Speedup (x)")
    ax2.set_title(f"Speedup vs Naive (N={matrix_size})")
    ax2.set_xticks(range(len(labels)))
    ax2.set_xticklabels(labels, rotation=25, ha="right", fontsize=8)
    ax2.axhline(y=1, color="gray", linestyle="--", alpha=0.5)
    for bar, val in zip(bars2, speedups):
        ax2.text(
            bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.1,
            f"{val:.1f}x", ha="center", va="bottom", fontsize=9
        )

    plt.tight_layout()
    plt.savefig("benchmark_results.png", dpi=150)
    print(f"\n  Chart saved to: benchmark_results.png")

What does generate_chart do?

matplotlib.use("Agg") selects a non-interactive backend that renders directly to a file. This avoids errors about missing display servers in WSL.

The function creates two side-by-side panels. The left panel shows raw execution time in milliseconds. The right panel shows speedup relative to the naive baseline. Both are saved as benchmark_results.png.

  • Now add the main() function directly below generate_chart. This orchestrates the full benchmark:

def main():
    matrix_size = 512
    runs = 5

    print(f"\n{'='*60}")
    print(f"  Matrix Multiplication Benchmark")
    print(f"  Matrix size: {matrix_size}x{matrix_size} | Runs per variant: {runs}")
    print(f"{'='*60}")

    print(f"\n  Compiling variants...\n")
    compile_program("matmul_naive.cpp", "matmul_naive", ["-O0"])
    compile_program("matmul_optimized.cpp", "matmul_optimized", ["-O0"])
    compile_program("matmul_naive.cpp", "matmul_naive_O2", ["-O2"])
    compile_program("matmul_optimized.cpp", "matmul_optimized_O2", ["-O2"])
    compile_program(
        "matmul_optimized.cpp", "matmul_optimized_O3",
        ["-O3", "-march=native"]
    )

    print(f"\n  Running benchmarks...\n")
    results = {}

    times = run_benchmark("matmul_naive", matrix_size, runs)
    results["Naive (ijk, -O0)"] = print_results("Naive (ijk, -O0)", times)

    times = run_benchmark("matmul_optimized", matrix_size, runs)
    results["Optimized (ikj, -O0)"] = print_results("Optimized (ikj, -O0)", times)

    times = run_benchmark("matmul_naive_O2", matrix_size, runs)
    results["Naive (ijk, -O2)"] = print_results("Naive (ijk, -O2)", times)

    times = run_benchmark("matmul_optimized_O2", matrix_size, runs)
    results["Optimized (ikj, -O2)"] = print_results("Optimized (ikj, -O2)", times)

    times = run_benchmark("matmul_optimized_O3", matrix_size, runs)
    results["Optimized (ikj, -O3 native)"] = print_results(
        "Optimized (ikj, -O3 native)", times
    )

    baseline = results["Naive (ijk, -O0)"]
    print(f"\n{'='*60}")
    print(f"  Speedup Summary (vs Naive ijk -O0 baseline)")
    print(f"{'='*60}\n")
    for label, avg in results.items():
        speedup = baseline / avg
        bar = "#" * min(int(speedup * 3), 50)
        print(f"  {label:.<45} {speedup:>5.1f}x  {bar}")

    print(f"\n{'='*60}")

    generate_chart(results, matrix_size)


if __name__ == "__main__":
    main()

What does main() do?

The function compiles 5 variants of the matrix multiply: naive with -O0, optimized with -O0, naive with -O2, optimized with -O2, and optimized with -O3 -march=native.

Each variant runs 5 times to smooth out measurement noise. The script then prints a comparison table showing average time and standard deviation, followed by a speedup summary with a visual bar chart in the terminal.

  • Save the file by pressing Ctrl+O, then Enter to confirm the filename.
  • Exit nano by pressing Ctrl+X.

Seeing a syntax error when you save?

Nano does not check Python syntax on save. You will only see errors when you run the script in the next substep. If you suspect a typo, reopen the file with nano benchmark.py and compare against the full code in the double-check tab below.

help me find a syntax error in my benchmark.py file

✔️ Awesome, I've got everything!

Double check you have saved the file and exited nano.

ⓧ I'd like to double check the full code

import subprocess
import sys
import statistics


def compile_program(source, output, flags=None):
    """Compile a C++ source file with specified flags."""
    cmd = ["g++", source, "-o", output]
    if flags:
        cmd.extend(flags)
    result = subprocess.run(cmd, capture_output=True, text=True)
    if result.returncode != 0:
        print(f"Compilation failed for {source}:")
        print(result.stderr)
        sys.exit(1)
    print(f"  Compiled: {output} (flags: {' '.join(flags) if flags else 'none'})")


def run_benchmark(executable, matrix_size, runs=5):
    """Run a benchmark executable multiple times and return timing results."""
    times = []
    for _ in range(runs):
        result = subprocess.run(
            [f"./{executable}", str(matrix_size)],
            capture_output=True, text=True
        )
        if result.returncode != 0:
            print(f"Error running {executable}: {result.stderr}")
            sys.exit(1)
        elapsed = float(result.stdout.strip())
        times.append(elapsed)
    return times


def print_results(label, times):
    """Print formatted benchmark results."""
    avg = statistics.mean(times)
    std = statistics.stdev(times) if len(times) > 1 else 0.0
    print(f"  {label:.<45} {avg:>8.1f} ms (std: {std:.1f} ms)")
    return avg


def generate_chart(results, matrix_size):
    """Generate a bar chart of benchmark results."""
    try:
        import matplotlib
        matplotlib.use("Agg")
        import matplotlib.pyplot as plt
    except ImportError:
        print("\n  matplotlib not available, skipping chart.")
        return

    labels = list(results.keys())
    values = list(results.values())
    baseline = values[0]
    speedups = [baseline / v for v in values]

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

    colors = ["#e74c3c", "#3498db", "#f39c12", "#2ecc71", "#9b59b6"]

    bars1 = ax1.bar(range(len(labels)), values, color=colors)
    ax1.set_ylabel("Time (ms)")
    ax1.set_title(f"Matrix Multiplication Time (N={matrix_size})")
    ax1.set_xticks(range(len(labels)))
    ax1.set_xticklabels(labels, rotation=25, ha="right", fontsize=8)
    for bar, val in zip(bars1, values):
        ax1.text(
            bar.get_x() + bar.get_width() / 2, bar.get_height() + 10,
            f"{val:.0f}", ha="center", va="bottom", fontsize=9
        )

    bars2 = ax2.bar(range(len(labels)), speedups, color=colors)
    ax2.set_ylabel("Speedup (x)")
    ax2.set_title(f"Speedup vs Naive (N={matrix_size})")
    ax2.set_xticks(range(len(labels)))
    ax2.set_xticklabels(labels, rotation=25, ha="right", fontsize=8)
    ax2.axhline(y=1, color="gray", linestyle="--", alpha=0.5)
    for bar, val in zip(bars2, speedups):
        ax2.text(
            bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.1,
            f"{val:.1f}x", ha="center", va="bottom", fontsize=9
        )

    plt.tight_layout()
    plt.savefig("benchmark_results.png", dpi=150)
    print(f"\n  Chart saved to: benchmark_results.png")


def main():
    matrix_size = 512
    runs = 5

    print(f"\n{'='*60}")
    print(f"  Matrix Multiplication Benchmark")
    print(f"  Matrix size: {matrix_size}x{matrix_size} | Runs per variant: {runs}")
    print(f"{'='*60}")

    print(f"\n  Compiling variants...\n")
    compile_program("matmul_naive.cpp", "matmul_naive", ["-O0"])
    compile_program("matmul_optimized.cpp", "matmul_optimized", ["-O0"])
    compile_program("matmul_naive.cpp", "matmul_naive_O2", ["-O2"])
    compile_program("matmul_optimized.cpp", "matmul_optimized_O2", ["-O2"])
    compile_program(
        "matmul_optimized.cpp", "matmul_optimized_O3",
        ["-O3", "-march=native"]
    )

    print(f"\n  Running benchmarks...\n")
    results = {}

    times = run_benchmark("matmul_naive", matrix_size, runs)
    results["Naive (ijk, -O0)"] = print_results("Naive (ijk, -O0)", times)

    times = run_benchmark("matmul_optimized", matrix_size, runs)
    results["Optimized (ikj, -O0)"] = print_results("Optimized (ikj, -O0)", times)

    times = run_benchmark("matmul_naive_O2", matrix_size, runs)
    results["Naive (ijk, -O2)"] = print_results("Naive (ijk, -O2)", times)

    times = run_benchmark("matmul_optimized_O2", matrix_size, runs)
    results["Optimized (ikj, -O2)"] = print_results("Optimized (ikj, -O2)", times)

    times = run_benchmark("matmul_optimized_O3", matrix_size, runs)
    results["Optimized (ikj, -O3 native)"] = print_results(
        "Optimized (ikj, -O3 native)", times
    )

    baseline = results["Naive (ijk, -O0)"]
    print(f"\n{'='*60}")
    print(f"  Speedup Summary (vs Naive ijk -O0 baseline)")
    print(f"{'='*60}\n")
    for label, avg in results.items():
        speedup = baseline / avg
        bar = "#" * min(int(speedup * 3), 50)
        print(f"  {label:.<45} {speedup:>5.1f}x  {bar}")

    print(f"\n{'='*60}")

    generate_chart(results, matrix_size)


if __name__ == "__main__":
    main()
Run the full benchmark

The harness will compile all 5 variants, run each one 5 times, and display the results. This takes a few minutes because it runs 25 total matrix multiplications.

  • Run the benchmarking harness by executing:
python3 benchmark.py

What should I see?

The script prints a compilation summary, then a timing table for each variant showing average milliseconds and standard deviation. At the bottom, a speedup summary shows how many times faster each variant is compared to the naive -O0 baseline.

You should see the optimized -O3 -march=native variant achieving 15-30x speedup over the naive baseline. The final line confirms the chart was saved.

  • Confirm the chart file was created by running:
ls benchmark_results.png

You should see benchmark_results.png listed in the output.

Script fails with a compilation or runtime error?

Make sure both matmul_naive.cpp and matmul_optimized.cpp exist in the same directory as benchmark.py. Run ls *.cpp to confirm.

If you see a Python syntax error with a line number, reopen benchmark.py in nano and check that line against the full code in the double-check tab above.

help me debug my benchmark.py script

Why automate benchmarking?

Running each variant by hand introduces human error. You might forget to recompile after a flag change, or test one variant when the machine is under heavier load. The harness eliminates these variables by running everything in sequence with consistent methodology.

Multiple runs per variant plus standard deviation gives you confidence that the results are real performance differences, not noise.

Secret mission

Implement Cache Blocking (Tiling)

You have already achieved impressive speedups with loop reordering and compiler flags. But there is still performance left on the table. Can you break the matrix into cache-sized tiles to squeeze out even more speed?

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely inside WSL on your local machine. There are no cloud resources or ongoing costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

  • WSL 2 Ubuntu distribution (contains all project files and installed packages).
  • Project directory ~/matmul-optimization/ (C++ source files, Python harness, compiled binaries, benchmark chart).
  • Packages installed via apt: g++ and python3-matplotlib.

Keep everything running

No action needed. Choose this if you are still actively building or plan to extend the project with the Secret Mission.

  • Your WSL Ubuntu distribution stays installed and ready to use at any time.
  • All project files, compiled binaries, and the benchmark chart remain in ~/matmul-optimization/.
  • WSL uses minimal resources when not running. It does not consume CPU or significant memory when idle.

Pause - I'll come back to this later

Shut down WSL to free up memory while keeping all your files intact for later.

  • Open PowerShell and shut down all WSL instances by running:
  • This stops the WSL virtual machine and frees up memory. Your files, packages, and configurations remain saved inside the Ubuntu distribution.
wsl --shutdown
  • To resume later, open your Ubuntu terminal from the Start menu and everything will be exactly as you left it.

Delete - I don't want to use this again

Remove all project resources and start fresh. This permanently deletes your files and the entire Ubuntu distribution.

  • Open your Ubuntu terminal and delete the project directory by running:
rm -rf ~/matmul-optimization
  • Remove the installed packages by running:
sudo apt remove g++ python3-matplotlib -y
  • If you want to remove the entire WSL Ubuntu distribution (this deletes all data, settings, and software inside it), open an administrator PowerShell and run:
wsl --unregister Ubuntu

This is permanent

Running wsl --unregister Ubuntu permanently deletes everything inside that distribution. Make sure you have copied out any files you want to keep before running this command.

  • If you no longer need WSL at all, you can uninstall it through Windows Settings.
  • Open Settings.
  • Select Apps, then Installed apps.
  • Search for Windows Subsystem for Linux and click Uninstall.

That's a Wrap!

That's a Wrap!

Nice work! You just built a complete performance engineering lab from scratch, proving with hard numbers that cache-aware programming can speed up a core AI computation by 15-30x.

Here's what you learned:

  • Written a C++ matrix multiplication program and demonstrated measurable 5-30x performance improvement through cache-aware loop reordering (ijk to ikj).
  • Built a Python benchmarking harness that automates compilation, execution, timing, and visualization of multiple optimization strategies using subprocess and matplotlib.
  • Created system architecture diagrams mapping each optimization technique to its target in the CPU memory hierarchy, from registers through L1/L2/L3 cache to RAM.
  • Secret Mission: Implemented cache blocking (tiling) to squeeze out additional performance by fitting sub-matrices into L1/L2 cache.

Ready to quiz yourself?