matrix · high

Matrix Multiplication

Matrix multiplication computes each output element as the dot product of a row of A and a column of B, so multiplying an n x k by a k x m matrix is O(n*k*m) — O(n^3) for square matrices. The naive i-j-p loop order strides down columns of B and thrashes the cache at scale; reordering to i-k-j makes the inner loop unit-stride, enabling prefetch and SIMD for a large speedup with an identical result. Strassen's algorithm lowers the exponent to ~2.807 but is rarely worth it; production uses tuned BLAS libraries (OpenBLAS, MKL) or Eigen.

The code

// C = A * B   (A is n x k, B is k x m, C is n x m)
for (int i = 0; i < n; ++i)
    for (int j = 0; j < m; ++j) {
        int s = 0;
        for (int p = 0; p < k; ++p)
            s += A[i][p] * B[p][j]; // dot of row i of A and col j of B
        C[i][j] = s;
    }

// Cache-friendly i-k-j order: innermost loop walks C and B by ROW
for (int i = 0; i < n; ++i)
    for (int p = 0; p < k; ++p)
        for (int j = 0; j < m; ++j)
            C[i][j] += A[i][p] * B[p][j]; // unit-stride on C and B

What this lesson walks through

1. 01Definition — each C cell is a row · a column
2. 02Compute one cell — 1·5 + 2·7 = 19
3. 03O(n³) — and loop order decides cache behavior
4. 04Gotcha — inner dims must match; not commutative

Matrix multiply isn't element-wise. Every cell C[i][j] is the DOT PRODUCT of row i of A with column j of B: multiply pairs and sum them. Picture row i lying across, column j standing up, meeting at one output cell.