Expand description
Register-tiled, panel-packed f32 GEMM — the BLAS-shaped dense route for hosts with no BLAS.
§Why this exists
On macOS the f32 dense path issues the reference’s own Accelerate SGEMM and is exact against
the oracle. Off that platform — Linux, and above all wasm, where no BLAS exists at all —
the same call degrades to a dot-product loop: for every output element, walk k and reduce.
That formulation re-reads the entire activation row once per output column and gets no reuse
out of the weights, which is why the browser codec measured 89.1 s of a 97.3 s frame (92%).
This is the standard answer, and it is what every serious GEMM does: hold an MR x NR tile of
the output in registers, stream one packed k-panel of the weights past it, and pay for each
loaded weight MR times instead of once.
§Why it is plain scalar Rust with no intrinsics
Doctrine #3: hand-rolled wide SIMD over scalar inner loops measured ~5x SLOWER than LLVM
autovectorization in the sibling repos. The inner loop below is a fixed-size [[f32; NR]; MR]
accumulator updated by a broadcast scalar — precisely the shape LLVM turns into NR/4 v128
multiply-adds per row with no help. The structure is the lever; the instruction selection is
the compiler’s job.
Ported from franken_numpy/crates/fnp-linalg/src/lib.rs (packed_gemm_serial_tiled), f64 to
f32, with the packing adapted to this project’s [n, k] weight layout.
§Exactness
Bit-identical to the scalar reference, and that is a design constraint rather than a happy
accident. Each output element accumulates over ascending k into its own slot, one f32 add
at a time — the same values in the same order as crate::f32ref’s scalar dot product. No
partial-sum splitting, no reassociation, no fused multiply-add. Blocking and packing change
only WHICH element is computed WHEN, never how any single element is summed.
That matters here more than speed: the current wasm path uses eight independent partial chains
(a different, non-reference reduction order), so adopting this kernel moves the codec CLOSER to
the reference while making it faster. packed_matches_scalar_bit_for_bit pins the claim.
Functions§
- linear_
packed out[m, n] = x[m, k] @ weight[n, k]^T + bias[n].