Optimize Matrix Multiplication Cache Hits in C

Optimizing matrix multiplication in C requires aligning computational access patterns with CPU cache architectures. Because standard matrix multiplication algorithms frequently access memory with large strides, they cause severe cache thrashing and memory latency bottlenecks. By understanding C's row-major storage format and employing techniques such as loop reordering, matrix transposition, cache blocking (tiling), and memory alignment, developers can maximize spatial and temporal locality, drastically increasing L1/L2 cache hit rates and computational throughput.

Understanding Row-Major Order in C

C stores multidimensional arrays in row-major order. In a contiguous two-dimensional array allocated as matrix[ROWS][COLS], elements of the same row are placed adjacent to each other in memory:

&matrix[i][j] + 1 == &matrix[i][j + 1]

When a CPU accesses a specific memory address, it fetches an entire cache line (typically 64 bytes, or 8 double-precision floats) from main memory into the cache. Accessing memory sequentially along rows yields spatial locality, meaning subsequent operations access data already loaded into the cache line (a cache hit). Accessing memory down columns introduces a stride of COLS * sizeof(element) bytes, discarding unused elements in the cache line and generating a cache miss on almost every access.

The Pitfall of Naive Multiplication (IJK Order)

Consider the textbook triple-loop implementation:

for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
        for (int k = 0; k < N; k++) {
            C[i][j] += A[i][k] * B[k][j];
        }
    }
}
  • Matrix A (A[i][k]): Advances row-wise along k (stride 1). Good spatial locality.
  • Matrix C (C[i][j]): Stays stationary in the innermost loop. Excellent temporal locality.
  • Matrix B (B[k][j]): Advances column-wise along k (stride \(N\)). Poor spatial locality.

For large matrices where a column exceeds the cache capacity, every iteration of k causes a cache miss for B[k][j], stalling the CPU.

Technique 1: Loop Reordering (IKJ Order)

Simply changing the nesting order of the loops eliminates large-stride memory access in the innermost loop without allocating extra memory:

for (int i = 0; i < N; i++) {
    for (int k = 0; k < N; k++) {
        double r = A[i][k]; // Loaded once into a register
        for (int j = 0; j < N; j++) {
            C[i][j] += r * B[k][j];
        }
    }
}

In the IKJ variation:

  • A[i][k] is constant relative to the innermost loop index j and stays in a CPU register.
  • B[k][j] and C[i][j] both advance sequentially along their rows (stride 1) as j increments.
  • This maximizes spatial locality for all memory operations, significantly reducing cache misses.

Technique 2: Explicit Matrix Transposition

If the loop structure cannot be altered, transposing the second matrix before multiplication converts column reads into row reads:

// Transpose B into B_T
for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
        B_T[j][i] = B[i][j];
    }
}

// Compute multiplication using row-major reads
for (int i = 0; i < N; i++) {
    for (int j = 0; j < N; j++) {
        double sum = 0.0;
        for (int k = 0; k < N; k++) {
            sum += A[i][k] * B_T[j][k]; // Both stride 1
        }
        C[i][j] = sum;
    }
}

Although transposing incurs an \(O(N^2)\) overhead and requires extra memory, it pays off for large matrices during the \(O(N^3)\) multiplication phase because both inputs are read with stride-1 access.

Technique 3: Cache Blocking (Tiling)

When matrices exceed the total capacity of the L1, L2, or L3 caches, elements loaded into the cache during early iterations are evicted before they can be reused. Cache blocking decomposes the \(N \times N\) matrices into smaller submatrices (\(B_{size} \times B_{size}\)) that fit entirely into the target cache level, maximizing temporal locality:

#define BLOCK_SIZE 32 // Chosen to fit in the L1/L2 cache

for (int sj = 0; sj < N; sj += BLOCK_SIZE) {
    for (int si = 0; si < N; si += BLOCK_SIZE) {
        for (int sk = 0; sk < N; sk += BLOCK_SIZE) {
            // Multiply sub-blocks
            for (int i = si; i < si + BLOCK_SIZE; i++) {
                for (int k = sk; k < sk + BLOCK_SIZE; k++) {
                    double r = A[i][k];
                    for (int j = sj; j < sj + BLOCK_SIZE; j++) {
                        C[i][j] += r * B[k][j];
                    }
                }
            }
        }
    }
}

Choosing BLOCK_SIZE depends on cache capacity: \[\text{Memory per block} = 3 \times (\text{BLOCK\_SIZE}^2) \times \text{sizeof(element)} \le \text{Cache Size}\]

Technique 4: Data Alignment and Padding

  1. Cache Line Alignment: Align allocated memory on 64-byte boundaries using posix_memalign() or aligned_alloc() to ensure rows start at the beginning of a cache line. This prevents unaligned access penalties and facilitates SIMD vectorization:

    double *A = aligned_alloc(64, N * N * sizeof(double));

  2. Array Padding: When matrix dimensions are powers of two (e.g., \(1024 \times 1024\) or \(2048 \times 2048\)), different rows map to the exact same cache sets in set-associative caches. This causes conflict misses, where addresses prematurely evict one another even if total cache capacity remains. Adding padding to the row stride (e.g., using a leading dimension of \(1024 + 8\)) distributes memory addresses across different cache sets, eliminating conflict misses.