Matrix Multiplication in C, Java & Python

1. What matrix multiplication actually is

Matrix multiplication is the operation of combining two matrices to produce a third. If A is an m × k matrix and B is a k × n matrix, their product C = A × B is an m × n matrix. The one shape rule to remember: the inner dimensions must match (the number of columns in A must equal the number of rows in B); the outer dimensions become the shape of the result.

Formally, each entry of the result is:

C[i][j] = Σ (k=0 .. k-1)  A[i][k] · B[k][j]

In words: to compute the number at row i, column j of the output, walk down row i of A and down column j of B at the same time, multiply matching pairs, and add them all up. That sum is C[i][j].

A worked 2×2 example

Multiply two 2×2 matrices. Here k = 2, so each dot product has two terms:

Let's compute the top-left entry C[0][0]: it's the dot product of row 0 of A (1, 2) with column 0 of B (5, 7):

C[0][0] = 1·5 + 2·7 = 5 + 14 = 19

The top-right entry C[0][1] uses the same row of A but column 1 of B (6, 8):

C[0][1] = 1·6 + 2·8 = 6 + 16 = 22

The bottom row follows the same pattern using row 1 of A (3, 4):

C[1][0] = 3·5 + 4·7 = 15 + 28 = 43
C[1][1] = 3·6 + 4·8 = 18 + 32 = 50

2. The algorithm we're measuring

There are six ways to order the three loops (i-j-k, i-k-j, j-i-k, …). They all compute the same mathematical result, but they touch memory in very different patterns — and that turns out to matter more than the language. This worksheet deliberately uses the textbook i-j-k order with a per-element accumulator:

for i in 0 .. n-1:
    for j in 0 .. n-1:
        s = 0
        for k in 0 .. n-1:
            s = s + A[i][k] * B[k][j]   # walk A's row, B's column
        C[i][j] = s

This is the most literal translation of the definition C[i][j] = Σ A[i][k]·B[k][j]: one result entry at a time, accumulating its dot product. Every language below runs this exact loop structure, unchanged, so that the only thing varying between runs is the language runtime — not the algorithm.

3. Interactive: from dot product to parallel blocks

The three optimizations ahead are all about how the work is shaped, not the math. The math is identical in every tab below — same C = A·B, same n³ multiply-adds on a 6×6 matrix. Step through each to feel the difference between computing one cell at a time, many cells at once across cores, and whole tiles that stay in cache.

4. Choosing the data size (and staying under 5 minutes)

The constraint: the slowest run must finish in under 5 minutes, single-threaded. Python's pure-interpreter triple loop is the binding constraint — C and Java finish in seconds. Because the cost grows as n³, the safest way to pick n is to measure at small sizes and extrapolate, rather than guess and risk a 6-minute Python run.

Calibration runs (Python, pure CPython)

Two things to notice. First, between n=200 and n=300 the scaling is almost perfectly cubic (3.49× observed vs 3.38× predicted) — the data still fits in cache. Second, by n=800 it's running super-linearly (101× vs the 64× that pure n³ would predict): the 15 MB working set no longer fits in this chip's L3 cache, so each access starts paying real memory-latency cost. That's a crucial detail: extrapolating from a cache-resident size understates the large-size cost.

So I extrapolate from the n=800 point — which already lives in the cache-miss regime, making its per-op cost representative of larger sizes:

time(n) ≈ 45.86 · (n / 800)³

Setting that equal to a 3.5-minute budget (210 s, leaving margin below 300 s) solves to n ≈ 1329. Rounding to a clean benchmark size gives n = 1300:

predicted Python time(1300) = 45.86 · (1300/800)³ ≈ 197 s  (~3.3 min)

5. The exact code

All three programs are structurally identical: deterministic fill, time only the multiply, print a checksum so cross-language correctness is verifiable. C and Java use flat 1-D arrays in row-major order; Python uses a list-of-lists (also row-major). The fill uses small deterministic values A[i][j] = (i+j)·0.001, B[i][j] = (i·j+1)·0.001 so the checksum is reproducible.

C — matmul.c

C · GCC 16/*
 * matmul.c - Single-threaded matrix multiplication (textbook i-j-k).
 *
 * Build: gcc -O3 -march=native matmul.c -o matmul
 * Run:   ./matmul [N]        (default N = 1000)
 *
 * Computes C = A * B for two N x N matrices of doubles using the
 * naive O(N^3) triple-nested loop with a per-element accumulator.
 * Only the multiply is timed; allocation/init/checksum are excluded.
 */

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main(int argc, char **argv) {
    int n = (argc > 1) ? atoi(argv[1]) : 1000;

    double *A = malloc((size_t)n * n * sizeof(double));
    double *B = malloc((size_t)n * n * sizeof(double));
    double *C = calloc((size_t)n * n, sizeof(double));
    if (!A || !B || !C) {
        fprintf(stderr, "allocation failed for n=%d\n", n);
        return 1;
    }

    /* Deterministic fill, identical across all three languages. */
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            A[(size_t)i * n + j] = (double)(i + j) * 0.001;
            B[(size_t)i * n + j] = (double)(i * j + 1) * 0.001;
        }
    }

    struct timespec t0, t1;
    clock_gettime(CLOCK_MONOTONIC, &t0);

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            double s = 0.0;
            for (int k = 0; k < n; k++) {
                s += A[(size_t)i * n + k] * B[(size_t)k * n + j];
            }
            C[(size_t)i * n + j] = s;
        }
    }

    clock_gettime(CLOCK_MONOTONIC, &t1);
    double secs = (t1.tv_sec - t0.tv_sec) +
                  (t1.tv_nsec - t0.tv_nsec) * 1e-9;

    /* Checksum: sum of all C elements, in row-major order (matches Java/Python). */
    double sum = 0.0;
    for (size_t i = 0; i < (size_t)n * n; i++) {
        sum += C[i];
    }

    printf("C       n=%d  time=%.4fs  checksum=%.6f\n", n, secs, sum);

    free(A);
    free(B);
    free(C);
    return 0;
}

Java — MatMul.java

Java · OpenJDK 21/*
 * MatMul.java - Single-threaded matrix multiplication (textbook i-j-k).
 *
 * Build: javac MatMul.java
 * Run:   java MatMul [N]      (default N = 1000)
 *
 * Computes C = A * B for two N x N matrices of doubles using the
 * naive O(N^3) triple-nested loop with a per-element accumulator.
 * Only the multiply is timed; allocation/init/checksum are excluded.
 * Uses 1-D arrays (row-major) so memory layout matches the C version.
 */

public class MatMul {
    public static void main(String[] args) {
        int n = args.length > 0 ? Integer.parseInt(args[0]) : 1000;

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

        // Deterministic fill, identical across all three languages.
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                A[i * n + j] = (i + j) * 0.001;
                B[i * n + j] = (i * j + 1) * 0.001;
            }
        }

        long t0 = System.nanoTime();

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

        long t1 = System.nanoTime();
        double secs = (t1 - t0) / 1e9;

        // Checksum: sum of all C elements, in row-major order (matches C/Python).
        double sum = 0.0;
        for (double v : C) {
            sum += v;
        }

        System.out.printf("Java    n=%d  time=%.4fs  checksum=%.6f%n",
                          n, secs, sum);
    }
}

Python — matmul.py

Python · CPython 3.14#!/usr/bin/env python3
"""
matmul.py - Single-threaded matrix multiplication (textbook i-j-k).

Run:  python3 matmul.py [N]    (default N = 1000)

Computes C = A * B for two N x N matrices of doubles using the naive
O(N^3) triple-nested loop with a per-element accumulator. Pure CPython,
no NumPy. Only the multiply is timed; allocation/init/checksum excluded.
Uses nested lists (list-of-lists) in row-major order.
"""

import sys
import time


def main() -> None:
    n = int(sys.argv[1]) if len(sys.argv) > 1 else 1000

    # Deterministic fill, identical across all three languages.
    A = [[(i + j) * 0.001 for j in range(n)] for i in range(n)]
    B = [[(i * j + 1) * 0.001 for j in range(n)] for i in range(n)]
    C = [[0.0] * n for _ in range(n)]

    t0 = time.perf_counter()

    for i in range(n):
        for j in range(n):
            s = 0.0
            for k in range(n):
                s += A[i][k] * B[k][j]
            C[i][j] = s

    t1 = time.perf_counter()
    secs = t1 - t0

    # Checksum: sum of all C elements, in row-major order (matches C/Java).
    total = 0.0
    for row in C:
        for v in row:
            total += v

    print(f"Python  n={n}  time={secs:.4f}s  checksum={total:.6f}")


if __name__ == "__main__":
    main()

6. Results

Same algorithm, same n = 1300, same checksum. Three languages, three radically different runtimes.

Throughput = 2n³ / time. Full checksum for all three: 1404884765266.165039 (Java shows …5000 at the 6th decimal — see the FMA note below). All three finish well under the 5-minute ceiling; Python is the only one anywhere near it.

Visual: linear scale (the gap is the point)

C and Java are the barely-visible slivers at the left edge — under 2.3% of Python's bar. The interpreted language is doing the same arithmetic ~56× slower.

Visual: zoomed to C vs Java

Among the compiled/JIT languages, C is ~24% faster than Java on this workload — a much smaller gap than the Python chasm.

7. Why the numbers look like this

Why Python is ~56× slower than C

CPython doesn't compile to machine code ahead of time. Each iteration of the inner loop s += A[i][k] * B[k][j] is executed by the interpreter: dispatch on the bytecode, look up A[i] (a list, then a list-item lookup), look up [k], look up B[k], look up [j], multiply, add, and store back — all through type-checked, reference-counted generic machinery. That's on the order of a few dozen nanometers of overhead per operation, where C does a single mulsd/addsd in a nanosecond. There is no vectorization and no escape from the type dispatch. Net effect: ~20 MFLOPS in Python versus ~1140 MFLOPS in C — a 56× gap that is almost entirely the cost of interpretation, not the cost of the math.

Why C and Java are within 24% of each other

Java's HotSpot JIT watches the inner loop, discovers it's hot, and compiles it to native x86-64 — at which point it's doing essentially the same mulsd/addsd sequence as the C binary. The ~24% gap is the JIT's residual overhead: array bounds checks that C omits, a slightly less aggressive loop optimizer than GCC 16 at -O3 -march=native, and the fact that the JIT must speculate and guard its optimizations. Notably, neither language manages to vectorize this loop, because the i-j-k access pattern reads B down a column (stride-n) — there's no contiguous run for SIMD to grab. So both compiled languages are stuck doing scalar floating point, which is exactly why they end up so close.

The real headline: everyone is at ~1–2% of the hardware's capability

A note on the checksum mismatch at the 6th decimal

C and Python report …65165039; Java reports …65165000. They agree to ~14 significant digits. The difference is fused multiply-add (FMA) contraction: HotSpot is allowed to rewrite s += a * b as a single FMA instruction fmadd, which rounds once instead of twice (multiply, then add). That changes the last bit of each accumulator, and the effect accumulates across n additions. C and Python perform separate multiply-then-add rounding. Neither is "wrong" — this is a classic floating-point subtlety worth knowing about when you compare numeric output across languages.

8. Environment & build flags

All builds use only each language's standard library — no external packages, no pip, no Maven. C uses clock_gettime(CLOCK_MONOTONIC); Java uses System.nanoTime(); Python uses time.perf_counter(). Only the multiply loop is timed; allocation, initialization, and checksum are excluded in every language.

9. Where the next 50× comes from

The baseline is set. Here's the optimization ladder the rest of the course climbs, each rung recoverable without changing the mathematics:

By the end of that ladder, the same mathematical operation that takes Python 215 seconds here can finish in well under a second. That ratio — the distance between a careless implementation and a careful one — is what performance engineering is about.