A dense matrix multiply is the most-optimized kernel in numerical computing, and for good reason: it is the beating heart of linear algebra, and it is the one workload where the gap between "what you naively wrote" and "what the hardware can do" is widest. In Equana — a scientific-computing notebook that runs entirely in the browser — A * B on two Float64 matrices dispatches to a WebAssembly kernel. This is the story of how that kernel went from a 3.1 GFLOPS JavaScript triple loop to 410 GFLOPS across 32 threads, on a Ryzen 9950X, one optimization at a time.

First, the unit. Multiplying two N×N matrices costs N³ multiply-add pairs, so 2·N³ floating-point operations. Divide that by the wall-clock time and you get GFLOPS:

GFLOPS = (2 · N³) / (time_seconds · 10⁹)

At N=1000 that's 2 billion flops per multiply. If you finish in 649 ms you're at ~3.1 GFLOPS; if you finish in 15 ms you're at ~130. Same arithmetic, same answer — the only thing that changes down this ladder is how well the code feeds the CPU. All numbers below are Float64, measured on the same machine.

Rung 0 — naive JavaScript: 3.1 GFLOPS

The baseline is the triple loop everyone writes first:

for (i) for (j) { s = 0; for (k) s += A[i][k] * B[k][j]; C[i][j] = s }

At N=1000 this takes about 649 ms. The problem isn't the arithmetic — it's everything around it. Arrays of arrays mean every A[i][k] is a pointer chase into a separate heap object, values are boxed doubles, and the inner loop walks B down a column, striding across memory by a full row each step. That column walk blows the cache on every iteration. 3.1 GFLOPS is what "just write the loop" buys you.

Rung 1 — typed arrays and a flat layout

The first fix is representation, not algorithm. Store each matrix as one contiguous Float64Array in column-major order — the same layout BLAS and LAPACK have used since Fortran — and index it as A[i + k*lda]. No boxing, no pointer chasing, no per-row heap objects. The JIT can now keep the accumulator in a register and the loads become predictable. This alone moves the needle meaningfully, but the fundamental memory-traffic problem remains: we still stream whole matrices through the cache faster than the cache can hold them.

Rung 2 — 64×64 cache tiling

The real bottleneck is that at N=1000 a single matrix is 8 MB — far larger than L1 or L2. If you compute C one full row at a time you touch all of B before you come back to reuse any of it, and by then it has been evicted. Loop tiling fixes this: break the problem into 64×64 sub-blocks and finish all the arithmetic for one block of C while its slices of A and B are still hot in cache. A 64×64 Float64 tile is 32 KB, so a working set of a few tiles fits comfortably in L1/L2. The flop count is identical — you've just reordered the loops so each value loaded from memory gets used many times before eviction. This is the single highest-leverage change in the whole diary, because it turns a memory-bound loop into a compute-bound one.

Rung 3 — a packed 4×4 kernel: 8.8 GFLOPS

Tiling fixes cache reuse between tiles; register blocking fixes reuse inside a tile. Instead of computing one C[i][j] at a time, compute a 4×4 patch of C in the innermost loop, holding all 16 accumulators live. Each value loaded from A and B now feeds four multiply-adds instead of one, so the load-to-arithmetic ratio drops by 4×. Combined with packing the operands into small contiguous scratch buffers — so the inner kernel reads them with unit stride and no lda multiply — the pure-JavaScript path reaches 8.8 GFLOPS, nearly 3× the naive baseline, with zero WASM involved. This is roughly the ceiling for scalar JavaScript: the JIT is now limited by issuing one floating-point op at a time.

Rung 4 — WASM naive: 4.7 GFLOPS

To go further we cross into WebAssembly. Surprisingly, a naive WASM triple loop lands at only 4.7 GFLOPS — barely ahead of tiled JavaScript. That's the shape of the tradeoff: WASM gives you predictable, un-deoptimizable execution and direct control over memory, but a naive kernel doesn't automatically beat a well-tuned JIT. WASM's advantage isn't raw scalar speed — it's that it exposes two things the JavaScript engine won't hand you: explicit SIMD and real threads. The rest of the ladder is about spending those.

Rung 5 — WASM SIMD: 9.1 GFLOPS

WebAssembly SIMD adds a 128-bit vector type, v128. For Float64 that's exactly two lanes — an f64x2. Rewriting the inner loop to load, multiply, and accumulate two doubles at a time with wasm_f64x2_mul and wasm_f64x2_add roughly doubles throughput to 9.1 GFLOPS. The shipped kernel builds its micro-tile out of these lanes: it holds a column of A as one f64x2 and broadcasts each B element across a vector with wasm_f64x2_splat, so a single vector multiply-add updates two rows of C at once. 128 bits is the widest lane WASM SIMD gives you — there is no f64x4 in the browser — and that ceiling shapes every decision that follows.

Rung 6 — full 4×4 SIMD: 35 GFLOPS

SIMD and register blocking multiply together. Take the 4×4 register-blocked kernel and make every accumulator a v128: now the innermost loop keeps a grid of vector accumulators live, each fed by a shared A vector and a splatted B scalar, and every instruction does two flops. With the accumulators resident in WASM's virtual registers and the operands packed for unit-stride access, the kernel jumps to 35 GFLOPS — a 3.8× leap over plain SIMD, and more than 10× the naive baseline. At N=1000 the SIMD kernel finishes a multiply in about 57 ms, down from 649. This is the point where the browser stops feeling like a toy for numerics.

Rung 7 — 6×8 / 8×8 register blocking

The remaining single-thread headroom is in the tile shape. A wider micro-kernel — 6×8 or 8×8 accumulators — amortizes each operand load across even more multiply-adds and gives the scheduler more independent accumulator chains to hide instruction latency. The art is not going too wide: WASM exposes a stack-machine bytecode, and the engine's register allocator has to fit your live accumulators into the CPU's actual physical registers. Ask for more live values than there are registers and they spill to memory, and your beautiful kernel starts thrashing the stack. The sweet spot is the largest tile whose accumulators still live entirely in registers — found empirically, not by formula.

Rung 8 — the SSE-style kernel: 37.9 GFLOPS single-thread

Tuning the tile, the packing, and the loop order together produces the production single-thread kernel: an SSE-style micro-kernel built on 128-bit f64x2 lanes, hand-shaped for WASM's constraints. At N=4000 it sustains 37.9 GFLOPS on one thread. That number matters beyond the ladder: OpenBLAS — the reference hand-tuned BLAS — compiled straight to WASM manages about 30.7 GFLOPS single-thread on the same problem. The kernel written for WebAssembly beats the world-class native library ported to WebAssembly by roughly 23%. That result gets its own post; the short version is that OpenBLAS's micro-kernels assume a native register file and hardware prefetch that don't survive the compile to a stack machine.

Rung 9 — pthreads: 130 GFLOPS at N=1000, 410 at N=4000

The last order of magnitude is parallelism. Matrix multiply is embarrassingly parallel across the columns of C, and WebAssembly threads — WASM's pthreads implementation over a shared SharedArrayBuffer — let the kernel split the output block across workers with no data copying, because every thread already sees the same linear memory. Partition the columns, run one instance of the single-thread kernel per worker, and throughput scales close to linearly with core count. At N=1000 the multi-threaded kernel finishes in about 15 ms — 130 GFLOPS. At N=4000, where there's enough work to keep all 32 threads saturated and amortize the spawn cost, it reaches 410 GFLOPS. For reference, a native SSE OpenBLAS build on the same machine sits around 400. The browser is not in a different league anymore; on this workload it's in the same seat.

The whole ladder

Read top to bottom, the pattern is consistent: cache tiling turns a memory-bound loop compute-bound, register blocking cuts the load-to-flop ratio, SIMD does two flops per instruction, and threads buy the last order of magnitude. Each technique is old — GEMM tuning is a decades-deep field — but the interesting part is how cleanly they all transfer to a 128-bit stack machine running inside a browser tab. Nothing here required leaving the sandbox.

Every rung on this ladder ships in Equana today; open a notebook at equana.dev and A * B runs the kernel from rung 9. If you want the architecture that makes zero-copy threading possible, read how Equana shares one WASM memory via SharedArrayBuffer; and for why the hand-written kernel outruns OpenBLAS-in-WASM, here's the full teardown.