A pedagogical SGEMM (single-precision matrix multiply) kernel in C++, shipped as five progressively-optimized variants so you can watch each optimization earn its keep. Each variant adds one technique compared to the previous file, and each file lives on its own so the diff is the optimization.
Final speedup on Intel i9-13900K, FP32, N=2048:
naive -> 0.5 GFLOPS (1.0x baseline)
+ ikj reorder -> 30.6 GFLOPS (60x)
+ cache blocks -> 52.4 GFLOPS (103x)
+ AVX2 4x16 -> 104.0 GFLOPS (204x)
+ OpenMP 32T -> 909.6 GFLOPS (1783x)
(numpy / OpenBLAS reference on the same box: 492 GFLOPS)
The parallel variant beats numpy's bundled OpenBLAS on this hardware --
not because the kernel is better engineered (OpenBLAS is decades of
professional work; this is ~80 lines of intrinsics), but because OpenBLAS's
default thread heuristic doesn't fully exploit the 13900K's hybrid
8 P-core + 16 E-core layout, while a flat #pragma omp parallel for does.
See the honest analysis below.
| # | File | Optimization added | GFLOPS (N=2048) |
|---|---|---|---|
| 1 | src/gemm_naive.cpp |
textbook ijk loop | 0.5 |
| 2 | src/gemm_reorder.cpp |
ikj loop order (stride-1 on B and C) | 30.6 |
| 3 | src/gemm_blocked.cpp |
three-level cache blocking (MC, NC, KC) | 52.4 |
| 4 | src/gemm_avx2.cpp |
hand-written 4x16 AVX2+FMA micro-kernel | 104.0 |
| 5 | src/gemm_parallel.cpp |
OpenMP across the M dimension | 909.6 |
Each file is self-contained and small (the AVX2 micro-kernel is ~30 lines of intrinsics). You should be able to read them top-to-bottom in 20 minutes.
The heart of the project. Holds a 4-row × 16-column tile of C in 8 YMM
accumulators while the inner k-loop streams 16 floats of B and 4 scalar
broadcasts of A per iteration:
B (k-th row, 16 floats)
+------+------+
| b0 | b1 | (2 ymm vectors)
+------+------+
| |
broadcast a0 ----> | FMA | a0*b0, a0*b1 -> c00, c01
broadcast a1 ----> | FMA | a1*b0, a1*b1 -> c10, c11
broadcast a2 ----> | FMA | a2*b0, a2*b1 -> c20, c21
broadcast a3 ----> | FMA | a3*b0, a3*b1 -> c30, c31
Each k-iteration does 8 FMAs producing 16 outputs (a 4×16 update),
keeping both FMA pipes on the P-core busy. The 8 accumulators stay in
registers across the entire kc loop; only b0/b1 and the four ai
broadcasts are reloaded each iteration. Loads come from B (stride-1,
ideal for L1) and A (broadcast of one scalar, also cheap).
Peak FP32 SGEMM per Raptor Lake P-core @ 5.5 GHz:
2 FMA pipes × 8-wide AVX2 × 2 ops (mul+add) × 5.5 GHz = 176 GFLOPS / core
We measure 130 GFLOPS single-threaded (gemm_avx2 at N=1024) -- about
74% of single-core peak for a clean intrinsic kernel without
prefetching, packing, or hand-tuned register tiles. Not bad for an
afternoon of work.
Block sizes are tuned to the i9-13900K cache hierarchy:
| Cache | Approx. size | Block dimension | Hosts |
|---|---|---|---|
| L1d | 48 KB | MC * KC = 12K floats = 48 KB |
A_block |
| L2 | 2 MB | KC * NC = 60K floats = 240 KB |
B_block |
| L3 | 36 MB (shared) | (whole working set) | full A, B, C |
KC N
+--------------+ +---------------------+
| | | B_block (KC x NC) |
KC | | KC | |
| A_block | -----> | in L2 |
MC | (MC x KC) | +---------------------+
| in L1d |
+--------------+ +---------------------+
M | C_block (MC x NC) |
+---------------------+
We chose MC=64, KC=192, NC=320. These are good (within 2x of optimal
on this CPU) but not heroic -- a real BLIS-style autotuner would sweep them
per-machine.
# Build everything with -O3 -march=native -mavx2 -mfma -fopenmp
make # or: mingw32-make on MSYS2
make test # correctness vs naive (4/4 pass to ~1e-6 ULP)
make bench N=1024 # GFLOPS for all five variants
make bench N=2048
make bench-numpy N=2048 # numpy/OpenBLAS referenceSingle-file compile if you prefer:
g++ -std=c++20 -O3 -march=native -mavx2 -mfma -fopenmp \
src/gemm_*.cpp bench/bench.cpp -o benchmake test checks every optimized variant against gemm_naive on a
deliberately-irregular shape M=200 N=200 K=137 (none of these dimensions
are multiples of MC, NC, or KC -- so the scalar fallback paths in
gemm_avx2.cpp get exercised).
[ok] gemm_reorder max|diff| = 3.815e-06 (tol=1e-02)
[ok] gemm_blocked max|diff| = 3.815e-06 (tol=1e-02)
[ok] gemm_avx2 max|diff| = 3.815e-06 (tol=1e-02)
[ok] gemm_parallel max|diff| = 3.815e-06 (tol=1e-02)
The 3.8e-6 is FP32 ULP-level error from accumulating 137 products of
[-1, 1] values -- not algorithmic drift.
| variant | time (s) | GFLOPS | vs naive |
|---|---|---|---|
gemm_naive |
2.3621 | 0.91 | 1.0x |
gemm_reorder |
0.0542 | 39.62 | 43.6x |
gemm_blocked |
0.0362 | 59.40 | 65.3x |
gemm_avx2 |
0.0165 | 130.18 | 143.2x |
gemm_parallel |
0.0026 | 822.16 | 904.3x |
| numpy (OpenBLAS) | 0.0048 | 448.25 | - |
| variant | time (s) | GFLOPS | vs naive |
|---|---|---|---|
gemm_naive |
33.6783 | 0.51 | 1.0x |
gemm_reorder |
0.5606 | 30.65 | 60.1x |
gemm_blocked |
0.3280 | 52.39 | 102.7x |
gemm_avx2 |
0.1652 | 104.02 | 203.9x |
gemm_parallel |
0.0189 | 909.56 | 1783.0x |
| numpy (OpenBLAS) | 0.0350 | 492.07 | - |
It's honest to address this -- a 200-line teaching kernel should not outperform OpenBLAS. What's happening:
-
OpenBLAS thread heuristic. The build of OpenBLAS bundled with
numpy 1.26.4on Windows picks a thread count that maxes out around 8. It also doesn't have specific tuning for Raptor Lake's hybrid 8 P-core + 16 E-core layout. -
Our OpenMP scheme is trivial. A flat
#pragma omp parallel for schedule(static)over the M dimension uses all 32 logical threads. With independent output stripes there's no contention. Sometimes brute simplicity wins. -
Single-thread comparison is the honest one. Our
gemm_avx2does 104-130 GFLOPS single-threaded. A single-threaded OpenBLAS call would do somewhere in 150-200 GFLOPS thanks to packing, prefetching, and register-level tuning. OpenBLAS is still better on a per-thread basis. We just out-scale it with naive parallelism on this specific CPU.
If you build OpenBLAS from source with make USE_OPENMP=1 NUM_THREADS=32
on the same machine, the numbers reverse.
- v0.2.0: hand-written packing of A and B into contiguous panels (BLIS step that we currently skip). This is what closes the per-thread gap with OpenBLAS.
- v0.2.0 alt: DGEMM (FP64) and HGEMM (FP16/BF16) variants.
- v0.3.0: AVX-512 micro-kernel for CPUs that have it (server Xeons, Zen 4/5, M4 in Rosetta).
- v0.3.0 alt: a 6x16 micro-kernel (BLIS preferred shape on Skylake-era).
- Goto, K. & van de Geijn, R. (2008). Anatomy of High-Performance Matrix Multiplication. ACM TOMS 34(3). paper
- Smith, T. M. et al. (2014). Anatomy of High-Performance Many-Threaded Matrix Multiplication. (BLIS paper.) paper
- Intel Intrinsics Guide -- the canonical reference for
_mm256_*names: https://www.intel.com/content/www/us/en/docs/intrinsics-guide - Marat Dukhan's gemmlowp internals -- great practical commentary on micro-kernel design.