From 0.25 to 3,749 GFLOPS —
Optimising GEMM Kernels Across the Memory Hierarchy

What was done

Implemented a naive scalar CPU matrix multiplication, then optimized the loop ordering to respect CPU cache lines.

The Key Technical Decision

Naïve i–j–k order causes a cache miss on matrix B on every inner iteration because B[k][p] strides across rows. Reordering the loops to i–k–j fixes this: the inner j-loop accesses B[k][j] sequentially, allowing it to stay within the L1 cache line.

// Naïve i-j-k (Bad cache locality on B)
for (int i = 0; i < m; i++)
  for (int j = 0; j < p; j++)
    for (int k = 0; k < n; k++)
      C[i*p+j] += A[i*n+k] * B[k*p+j];

// Cache-aware i-k-j (Good cache locality)
for (int i = 0; i < m; i++)
  for (int k = 0; k < n; k++)
    for (int j = 0; j < p; j++)
      C[i*p+j] += A[i*n+k] * B[k*p+j];

What was done

Introduced multi-level loop tiling (blocking) to keep active working sets entirely within the L1 cache, combined with Advanced Vector Extensions (AVX2) for Single Instruction Multiple Data (SIMD) execution.

The Key Technical Decision

Tiles were sized explicitly to fit in L1 data cache. For float32, T=64 was chosen: 3 matrices × (64×64) × 4 bytes = 48 KB, precisely matching the Intel P-core L1 cache size. The tile size was empirically swept over T ∈ {32, 64, 128, 256}. Under AVX2, the _mm256_fmadd_ps intrinsic executes 8 float Fused Multiply-Adds (FMAs) per cycle, with a horizontal reduction performed only once per tile row.

// AVX2 SIMD Core Computation
__m256 a_val = _mm256_set1_ps(A[i*n+k]);
for (int j = 0; j < T; j+=8) {
  __m256 b_vec = _mm256_loadu_ps(&B[k*p+j]);
  __m256 c_vec = _mm256_loadu_ps(&C[i*p+j]);
  c_vec = _mm256_fmadd_ps(a_val, b_vec, c_vec);
  _mm256_storeu_ps(&C[i*p+j], c_vec);
}

What was done

Distributed the blocked AVX2 workload across multiple CPU cores using OpenMP.

The Key Technical Decision

The outer loop was parallelised using #pragma omp parallel for with static scheduling. Static scheduling was chosen over dynamic because all tile rows have identical computational cost; dynamic scheduling overhead would be pure waste. Due to the Intel i5-12th Gen P-core/E-core asymmetry, work was pinned to P-cores only.

False sharing investigation: Early runs showed sub-linear scaling. Performance profiling revealed cache-line invalidation traffic from adjacent C tiles being updated by different cores. Aligning the tiles to 64-byte boundaries resolved the false sharing hazard.

// Parallelizing the outer tile loops
#pragma omp parallel for schedule(static)
for (int ii = 0; ii < m; ii += T) {
  for (int kk = 0; kk < n; kk += T) {
    for (int jj = 0; jj < p; jj += T) {
      compute_tile_avx2(ii, kk, jj, T);
    }
  }
}

What was done

Transitioned execution to the NVIDIA GPU. First, implemented shared-memory tiling, then escalated to utilize hardware Tensor Cores via the WMMA (Warp Matrix Multiply Accumulate) API.

The Key Technical Decision

Shared-memory tiling: Thread blocks load 16×16 tiles of matrices A and B from global DRAM into fast SRAM (shared memory) to compute partial products. This reduces global memory accesses from O(n³) to O(n³/T), delivering a 23× speedup over the naïve global-memory CUDA kernel.

WMMA: Replaced thread-level scalar math with warp-level 16×16×16 FMA fragments via the nvcuda::wmma API. D = A·B + C is executed in a single warp-synchronous operation at FP16 precision. Required defining correct fragment types (matrix_a, matrix_b, accumulator), managing shared-memory double-buffering to overlap compute with data loading, and rigorous synchronization barriers. The final throughput sits ~6% below theoretical peak due to shared-memory bank conflicts and boundary tile handling.

// WMMA Tensor Core Fragment Math
using namespace nvcuda;
wmma::fragment16, 16, 16, half, wmma::row_major> a_frag;
wmma::fragment16, 16, 16, half, wmma::col_major> b_frag;
wmma::fragment16, 16, 16, float> c_frag;

wmma::fill_fragment(c_frag, 0.0f);

// Load from shared memory
wmma::load_matrix_sync(a_frag, shared_A, 16);
wmma::load_matrix_sync(b_frag, shared_B, 16);

// Hardware FMA
wmma::mma_sync(c_frag, a_frag, b_frag, c_frag);

// Store back to global memory
wmma::store_matrix_sync(&C[idx], c_frag, p, wmma::mem_row_major);