kernels(cpu)(matmul): F32 outer-product tiled GEMM kernel for prefill

[jamesburton]

Closes #312.

ROADMAP reference

This PR addresses the F32 portion of ROADMAP Phase 3 Step 26 — Outer-product tiled matmul kernels (currently unticked).

Step 26: Replace inner-product GEMM with outer-product formulation: unroll M×N output tile, share one activation load across all weight-row dot products. AVX2: 4×3 tile (12 ymm accumulators)… Blocked on RyuJIT register pressure — 4×3 tile needs 23 YMM registers (only 16 available), causing spills that negate the weight-reuse benefit.

Why F32 (and not Q8_0)

PR #61 landed the Q8_0 / Q5_0 / K-quant outer-product kernels but reverted GEMM dispatch — the 4×3 Q8_0 microkernel was ~14% slower than 3× VecDot4RowsR4 on AVX2 because RyuJIT spilled 7+ YMM registers (23 needed: 12 acc + 6 token + 1 ones + 4 data).

That register pressure is quantization-specific. The F32 path has no ones mask, no Q8_0 scale extraction, and no Half→float conversion. A 4×3 F32 tile budgets to exactly 16 YMM (12 accumulators + 3 token vectors + 1 reloaded weight vector), which fits the AVX2 register file naturally.

Kernel design

OuterProductGemm.OuterProductGemmF32 in a new file src/DotLLM.Cpu/Kernels/OuterProductGemm.cs:

• 4 weight rows × 3 input tokens register tile, applying dotLLM's C[N,M] = B[N,K] × A[M,K]^T convention.
• Vectorises along K (the contraction axis, 8 floats/lane via Vector256<float>). The maintainer's literal ROADMAP description envisions broadcast-style outer-product à la BLIS, which requires per-tile A-packing/transpose. Vectorising along K instead achieves the same register-reuse property (each load is reused across the orthogonal-axis FMAs in the tile) while avoiding a separate pack pass — and still ships within the literal "12 YMM accumulators" budget.
• At each K-step: 3 token vectors are held in registers and reused across 4 FMAs per row; the weight vector is reloaded for each row and reused 3× across the token FMAs.
• Horizontal-reduces 12 accumulators into the C tile at the end of the K loop.
• Scalar K-tail (k % 8), row tail (m % 4), and token tail (n % 3) handle non-tile-aligned shapes.
• AVX2/FMA detection at the public entry point; falls back to a scalar reference (OuterProductGemmF32Scalar) otherwise.

[SkipLocalsInit] + [MethodImpl(AggressiveOptimization)] on the AVX2 microkernel; [AggressiveInlining] on the helpers per project convention.

Tests

tests/DotLLM.Tests.Unit/Cpu/Kernels/OuterProductGemmTests.cs — 25 new test cases extending the existing Q8_0 outer-product test class.

• Scalar path: bit-exact equality with MatMul.GemmF32Scalar (Assert.Equal(reference, value) — no tolerance). The scalar microkernel accumulates the same terms in the same ascending-index order as the production scalar reference, and .NET does not auto-contract mul + add to FMA, so true bit equality is required and achieved.
• Vector path: tolerance ≤ 4e-6 · √K vs MatMul.GemmF32. The √K scaling comes from the standard error of an unbiased K-term sum for unit-magnitude operands with float32 ULP ≈ 1.2e-7. Mirrors the tolerance convention used in the existing Q8_0 outer-product tests (1e-2f/1e-3f).
• Coverage: tile-aligned shapes (M=4, N=3, K∈{8,64}; M=8, N=6; M=12, N=9; up to M=128, N=32, K=1024), all-tail-combination shapes (row tail only; token tail only; both; both plus K tail), and edge cases (M=1, N=1, K=1, pure-inner-product M=1 K=4096).
• All 53 outer-product tests pass (existing 28 Q8_0 + new 25 F32). No existing tests modified.

Benchmark

benchmarks/DotLLM.Benchmarks/OuterProductGemmF32Benchmark.cs — OuterProductGemmF32 vs production MatMul.GemmF32 at three prefill profiles (M ∈ {128, 512, 2048}, K=4096, N=32).

Intel Core Ultra 7 155H (AVX2, no AVX-512), single-threaded, BenchmarkDotNet warmupCount=5, iterationCount=15:

M	Baseline GemmF32	OuterProductGemmF32	Speedup
128	2.507 ms ± 0.189	1.238 ms ± 0.107	2.02×
512	10.216 ms ± 1.889	5.588 ms ± 0.580	1.83×
2048	47.070 ms ± 11.420	32.714 ms ± 2.369	1.44×

GFLOPS at M=512 jumps from ~13 to ~24; at M=2048 from ~11 to ~16. Speedup narrows at large M as memory bandwidth (rather than register-reuse) becomes the dominant cost — consistent with outer-product GEMM theory.

This is a single-threaded kernel-vs-kernel comparison. The pooled GemmF32(..., ComputeThreadPool) overload — which is what production prefill actually dispatches — is not touched in this PR; caller wiring is the natural follow-up.

Scope (explicit)

In scope:

• New OuterProductGemm kernel (scalar + AVX2).
• Parity tests + benchmark.

Out of scope (deliberate follow-ups):

• Wiring TransformerModel (or any caller) to the new kernel. Caller switching belongs in a separate PR after the benchmark is reproduced on a second machine and validated against the threaded path.
• Q8_0 / Q5_0 / K-quant outer-product re-enablement. Those kernels exist in MatMul.cs (PR Step 26: Outer-product tiled matmul — kernels + investigation (blocked on AVX2) #61) but remain blocked on AVX2 register pressure as documented in the roadmap. AVX-512 (32 ZMM) or a native C microkernel is a separate investigation.

Verification

• dotnet build src/DotLLM.Cpu/DotLLM.Cpu.csproj -p:EnableSourceControlManagerQueries=false -c Release — green.
• dotnet build tests/DotLLM.Tests.Unit/DotLLM.Tests.Unit.csproj -p:EnableSourceControlManagerQueries=false -c Release — green.
• dotnet test tests/DotLLM.Tests.Unit/ -c Release --filter "FullyQualifiedName~OuterProductGemmF32" — 25/25 pass.
• dotnet test tests/DotLLM.Tests.Unit/ -c Release --filter "FullyQualifiedName~OuterProduct" — 53/53 pass (existing 28 + new 25).
• dotnet build benchmarks/DotLLM.Benchmarks/DotLLM.Benchmarks.csproj -p:EnableSourceControlManagerQueries=false -c Release — green.
• Benchmark --filter "*OuterProductGemmF32Benchmark*" produces the table above.

Inspiration: llamafile / tinyBLAS.

[Copilot]

Pull request overview

Note

Copilot was unable to run its full agentic suite in this review.

Adds a new F32 outer-product (4×3) tiled GEMM kernel for prefill workloads, along with parity tests and a benchmark to compare against the existing GEMM implementation.

Changes:

• Introduces OuterProductGemm.OuterProductGemmF32 (AVX2/FMA) plus a scalar reference implementation.
• Adds unit tests validating scalar bit-exact parity and vector-mode tolerance-based parity vs existing GEMM paths.
• Adds a BenchmarkDotNet benchmark to compare the new kernel against MatMul.GemmF32 on prefill-shaped sizes.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 4 comments.

File	Description
tests/DotLLM.Tests.Unit/Cpu/Kernels/OuterProductGemmTests.cs	Adds parity, tail-handling, and edge-case tests for the new F32 outer-product GEMM kernels.
src/DotLLM.Cpu/Kernels/OuterProductGemm.cs	Implements scalar + AVX2/FMA outer-product GEMM microkernels (4×3 tile) and helpers.
benchmarks/DotLLM.Benchmarks/OuterProductGemmF32Benchmark.cs	Adds a benchmark comparing the new kernel vs the existing production GEMM for prefill-shaped workloads.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[jamesburton]

Thanks — went through all four review comments. Pushed 7c34e46 addressing them.

1. HorizontalSum uses Sse3.HorizontalAdd under an Avx2/Fma-only guard — fixed (with a caveat on the premise).

The premise is not quite right: .NET's x86 ISA classes form a hierarchy (Avx2 : Avx : Sse42 : Sse41 : Ssse3 : Sse3 : …), and the JIT sets support flags hierarchically, so Avx2.IsSupported == true implies Sse3.IsSupported == true. There was no reachable runtime throw. (Note the same would have applied to the Sse.Add on the line above, which the comment didn't flag.)

That said, the suggestion is still an improvement, so I took it — the reduction is now expressed on the cross-platform surface and depends on no x86 ISA class at all:

return Vector128.Sum(v.GetLower() + v.GetUpper());

This also matches the repo convention of preferring Vector128<T>/Vector256<T> over platform-specific intrinsics unless the platform-specific form is measurably faster. Since this changes the summation order of the final 4-wide reduction, I re-ran the full parity suite against both the scalar reference and MatMul.GemmF32 — all shapes still within tolerance.

2 + 3. Tolerance commentary is internally inconsistent — valid, fixed.

Both comments describe the same defect (the header block; the in-body comment at the second location was already correct). The header claimed a relative tolerance growing linearly in K, while the code uses an absolute tolerance 4e-6f * MathF.Sqrt(k). The code is the intended behaviour — √K is the right growth rate for the error of an unbiased K-term float32 sum — so I corrected the comment rather than the tolerance, and made it explicit that the scalar path is compared bit-exactly while only the vector path uses the tolerance.

4. Assertion-free return on non-AVX2 hardware reports as passed — valid, fixed.

Switched to the Xunit.SkippableFact mechanism already used elsewhere in this test project (e.g. the CUDA tests' Skip.IfNot(CudaDevice.IsAvailable(), …)), so a non-AVX2 agent now reports skipped:

[SkippableTheory]
public void OuterProductGemmF32_Avx2_MatchesReference(int m, int n, int k)
{
    Skip.IfNot(Avx2FmaSupported, "AVX2/FMA not supported on this machine");
    …
}

The guard moved out of the shared RunF32ParityCase helper so it can never silently no-op. The two mixed theories (_HandlesAllTails, _EdgeCases) exercise both paths and their scalar half always asserts, so those stay plain [Theory] with only the vector half gated — skipping them entirely on a non-AVX2 agent would lose real scalar coverage.

Verification: dotnet test --filter "FullyQualifiedName~OuterProductGemm" → 53 passed, 0 failed, 0 skipped (this machine has AVX2/FMA).