From 049e1222c60e0039b1e60317b67490e7342c5a48 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 20:38:05 +0000 Subject: [PATCH 1/7] =?UTF-8?q?blackboard:=202026-07-06=20review=20pass=20?= =?UTF-8?q?=E2=80=94=20health=20green,=2010=20findings=20on=20subpel=5Ftap?= =?UTF-8?q?=5Ftile=20(#235)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Review-only session. fmt/clippy/2175 lib tests green; example runs on AMX (TDPBF16PS tier, asserts pass). Findings recorded: unguarded Gotcha-14 asserts, throughput leg measures pack/alloc overhead (PackedBf16B is the fix), false positive-operand comment, SplitMix64/gemm_f64 reuse (gemm_f64 measured bit-exact vs naive reference, 0/256 lanes differ), misc cleanups. Refuted: .gitattributes-accidental, needless_range_loop, checksum, transpose reuse. --- .claude/blackboard.md | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 61d0e1fc..e194da6d 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3,6 +3,45 @@ > **Read this first.** The "Polyglot Notebook" architecture below is a > separate/older program, not the current epoch. +## 2026-07-06 — Review pass: health check green + 10 findings on subpel_tap_tile (#235) + +Review-only session (no code changes). **Health check:** `cargo fmt --check` +clean, `cargo clippy -p ndarray --lib -- -D warnings` clean, clippy on the new +example clean, **2175/2175 lib tests pass** (30 ignored). The stale top-of- +CLAUDE.md "build fails (exit 101)" note again does NOT reproduce. The example +`subpel_tap_tile` runs green end-to-end on this EMR host (tier = AMX TDPBF16PS, +rel err 0.157% / 0.215%, asserts pass). + +**Findings on `examples/subpel_tap_tile.rs` (PR #235), most severe first:** +1. Lines 208-209 hard-assert rel err < 0.05 on AMX-dispatched results with no + contention guard — flakes on oversubscribed VMs per Gotcha 14 (precedent: + `#[ignore]` gating in bf16_tile_gemm.rs:572). +2. Throughput leg (3) times per-call allocs + f32→bf16 of BOTH operands + the + kernel's per-call VNNI pack of the CONSTANT H — printed 1.00 M/s measures + wrapper overhead, not the tile op; use `PackedBf16B` + + `bf16_tile_gemm_16x16_packed` (the API built for exactly this). +3. Line 167 comment "(fits u8 and i8; positive operand)" is false (x reaches + ≈ −19.6 for r ≥ 11) — porting hazard toward the u8×i8 int8 path. +4. Reuse: `mix()` is the 5th example-local splitmix64 copy (public bit-identical + `hpc::cam_index::SplitMix64` exists); `direct_matmul` duplicates + `backend::native::gemm_f64` — **[MEASURED this session] bit-exact, 0/256 + lanes differ** vs the naive loop on the probe's exact operands (operator + correction: earlier "independence" refutation was wrong — the BF16 tile + kernel shares zero code with gemm_f64). +5. Cleanups: Vec C + copy (kernel takes &mut [f32] — stack array works), + to_bf16 double-alloc, b_pad identical-index loop = prefix copy, direct_hv + duplicated FIR pass, padded-16→32 scaffold now copy-pasted across 2 probes + (suggest one public padded helper on ndarray::simd before probe #3). + +**Refuted during verify:** .gitattributes deletion is the deliberate, +documented revert PR #236 (union merge mangles [[example]] blocks — no residue +found); clippy needless_range_loop claim (empirically clean under -D warnings); +f32 checksum absorption (print-only anti-DCE); transpose_matrix reuse (net +loss — Vec round-trips vs 8-line stack helper). + +[LOOSE END] Findings 1-3 are worth a small follow-up PR (contention guard or +warning, packed-B throughput leg, comment fix); none applied here (review-only). + ## 2026-07-02 (later) — bf16 tile GEMM: VDPBF16PS middle tier + PackedBf16B (loose end closed) Closed the [LOOSE END] from the 1BRC entry below. `hpc/bf16_tile_gemm.rs` From c90651367eb31a31a4d0fad86689c529f9e71207 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 20:42:54 +0000 Subject: [PATCH 2/7] =?UTF-8?q?blackboard:=20addendum=20=E2=80=94=20own=20?= =?UTF-8?q?gemm=5Ff64=5Ftiled=20also=20bit-exact;=20public=20f64=20GEMM=20?= =?UTF-8?q?surface=20is=20external-backed?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Operator challenge answered with measurement: naive reference vs the crate's own scalar gemm_f64_tiled (native.rs:473) is bit-exact (0/256 lanes, probe operands and random f64 K=64), same as vs matrixmultiply. Structural gap surfaced: gemm_f64_tiled is dead code in private mod scalar while every public f64 GEMM path delegates to the external matrixmultiply crate. --- .claude/blackboard.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index e194da6d..fb5de613 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -42,6 +42,21 @@ loss — Vec round-trips vs 8-line stack helper). [LOOSE END] Findings 1-3 are worth a small follow-up PR (contention guard or warning, packed-B throughput leg, comment fix); none applied here (review-only). +**[ADDENDUM, same day — operator challenge "our gemm, not the stoneage +external":** the first bit-exactness probe compared naive-loop vs +`backend::native::gemm_f64`, which delegates to the EXTERNAL +`matrixmultiply::dgemm` (native.rs:249; registry dep, Cargo.toml:154). Re-ran +against the crate's OWN scalar `gemm_f64_tiled` (native.rs:473, verbatim): +**bit-exact too — 0/256 lanes differ** on the probe's operands AND on random +f64 at K=64. All three (naive / own tiled / matrixmultiply) agree bit-for-bit +on these shapes. **Structural finding the challenge surfaced:** the crate's +entire PUBLIC f64 GEMM surface (`gemm_f64`, `BlasLevel3::blas_gemm` via +`backend_gemm`) is external-backed, while the own-Rust `gemm_f64_tiled` sits +dead in the private `mod scalar` with zero callers. If the policy is +own-reverse-engineered-Rust-only, the right fix for finding #4 is to surface +`gemm_f64_tiled` (or route the scalar tier of `backend_gemm` through it) and +point probes at THAT — files under the UNUSED_INVENTORY dead-code thread.]** + ## 2026-07-02 (later) — bf16 tile GEMM: VDPBF16PS middle tier + PackedBf16B (loose end closed) Closed the [LOOSE END] from the 1BRC entry below. `hpc/bf16_tile_gemm.rs` From 5bbcfc64f5396223e8944c85526977754bc02a29 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 21:12:38 +0000 Subject: [PATCH 3/7] simd: surface crate-native gemm_f64_tiled on the canonical ndarray::simd surface MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Graduates the tiled f64 GEMM from dead code (backend/native.rs private mod scalar, zero callers) to src/simd_ops.rs, re-exported unconditionally via ndarray::simd. Innermost j-loop vectorized on the dispatched F64x8 — one source serving AVX-512/AVX2/NEON/WASM-SIMD128/scalar. Bit-exactness contract (documented): per-element c = c + (alpha*a)*b in ascending-k order, mul and add unfused on every backend — bit-identical across backends, and at alpha=1 beta=0 bit-identical to the naive triple-loop reference. This makes it the in-crate ground-truth GEMM for probe certification, replacing reliance on the external matrixmultiply engine behind backend::native::gemm_f64. W1a: 7 parity tests (fixed-seed corpus, 13 shapes incl. multi-tile 70^3, strided leading dims + sentinel padding, alpha/beta semantics, beta=0 over NaN, k=0, m/n=0, denormals) — all to_bits equality. Zero unsafe. Bench vs naive scalar (single thread, EMR): 128^3 2.23x, 256^3 2.54x, 512^3 6.74x, bit-equal at every size. Consumer sites: direct_matmul in examples/subpel_tap_tile.rs and examples/gridlake_field_tile.rs rewired to it (bit-identical swap — identical printed output pre/post). Dead native.rs copy removed. Gates: clippy -D warnings clean (lib + examples), 2182/2182 lib tests, doctest green, no_std build green. --- .claude/blackboard.md | 56 ++++++ examples/gridlake_field_tile.rs | 16 +- examples/subpel_tap_tile.rs | 15 +- src/backend/native.rs | 50 +---- src/simd.rs | 8 + src/simd_ops.rs | 335 ++++++++++++++++++++++++++++++++ 6 files changed, 415 insertions(+), 65 deletions(-) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index fb5de613..ec4212f9 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3,6 +3,62 @@ > **Read this first.** The "Polyglot Notebook" architecture below is a > separate/older program, not the current epoch. +## 2026-07-06 (later) — `ndarray::simd::gemm_f64_tiled` surfaced (operator directive) + +The crate-native tiled f64 GEMM graduated from dead code +(`backend/native.rs` private `mod scalar`, zero callers) to the canonical +simd surface: `src/simd_ops.rs::gemm_f64_tiled`, re-exported +**unconditionally** in `src/simd.rs` (alloc-free; `pub mod simd` itself is +std-gated in lib.rs — see reviewer note below if that changed). + +- **Bit-exactness contract (documented on the fn):** every C[i,j] gets + `c = c + (α·A[i,p])·B[p,j]` ascending-p, mul and add UNFUSED — the + `*`/`+` operators on F64x8 lower to plain mul/add intrinsics on ALL + five backends (AVX-512 `_mm512_mul/add_pd`, AVX2 per-half, NEON + `vmulq/vaddq_f64`, WASM `f64x2_mul/add`, scalar) and Rust never + FP-contracts explicit intrinsics → bit-identical across backends; at + α=1 β=0 bit-identical to the naive triple-loop reference. +- Innermost j-loop vectorized on dispatched `F64x8` (one source, every + backend, per the simd_ops polyfill model); TILE=64 blocking preserved + verbatim from the original. +- [MEASURED] vs naive scalar triple loop, single thread, this EMR VM: + 128³ 2.23×, 256³ 2.54×, 512³ 6.74× (4.6 GFLOP/s), **bit-equal: true** + at every size (W1a bench criterion; well above the 0.5× reject line). +- W1a compliance: parity tests = 7 new tests in + `simd_ops::gemm_f64_tiled_tests` (fixed-seed splitmix64 corpus, 13 + shape sweep incl. multi-tile 70³, strided lda/ldb/ldc with + sentinel-padding assert, α/β semantics, β=0-over-NaN, k=0, m/n=0, + denormals/−0.0) — all bit-equality (`to_bits`), not tolerance. Zero + `unsafe`. No new feature detection. Consumer sites named: the + `direct_matmul` f64 ground-truth in `examples/subpel_tap_tile.rs` + + `examples/gridlake_field_tile.rs`, both REWIRED to it (bit-identical + swap — subpel prints identical numbers pre/post). +- Free-fn shape note: matches the existing simd-surface GEMM family + (`bf16_tile_gemm_16x16`, `matmul_i8_to_i32`) — operator-directed, not + a speculative W1a-queue addition. +- Dead code removed: `native.rs` `gemm_f64_tiled` deleted (pointer + comment left); `gemm_f32_tiled` stays dead in `mod scalar` pending the + same treatment (UNUSED_INVENTORY thread). +- Gates: clippy `-D warnings` clean (lib + both examples), fmt clean, + **2182/2182 lib tests**, doctest green, `--no-default-features` build + green. + +[NOTE] `--tests` clippy surfaces 3 PRE-EXISTING test-code lints +(property_mask.rs:426 unusual_byte_groupings, bitwise.rs:637 identity_op, +palette_codec.rs:806 needless_range_loop) — not touched here; the house +gate (`cargo clippy -- -D warnings`, no --tests) is clean. + +[OBSERVED] `amx_available()` flipped true→false on this VM between runs +two hours apart (subpel ran AMX TDPBF16PS earlier, F32x16 polyfill now; +identical BF16-class errors either way). Host/hypervisor XTILEDATA +enablement drift — Gotcha 9/14 adjacent; worth a note in AMX_GOTCHAS if +observed again. + +[LOOSE END] subpel_tap_tile findings 1-3 from the same-day review still +queued (Gotcha-14 assert guard, PackedBf16B throughput leg, positive- +operand comment fix). Adversarial 3-angle verify workflow on this diff +was in flight at commit time; findings (if any) land as follow-up. + ## 2026-07-06 — Review pass: health check green + 10 findings on subpel_tap_tile (#235) Review-only session (no code changes). **Health check:** `cargo fmt --check` diff --git a/examples/gridlake_field_tile.rs b/examples/gridlake_field_tile.rs index b980dc45..7c57a81d 100644 --- a/examples/gridlake_field_tile.rs +++ b/examples/gridlake_field_tile.rs @@ -81,18 +81,14 @@ mod amx { out } - /// Direct `f64` reference GEMM `C = A·B` (16×16). + /// Direct `f64` reference GEMM `C = A·B` (16×16) via the crate-native + /// bit-exact tiled GEMM (`ndarray::simd::gemm_f64_tiled`) — bit-identical + /// to the naive triple loop at α=1 β=0, entirely in-crate. fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] { + let af: Vec = a.iter().map(|&v| v as f64).collect(); + let bf: Vec = b.iter().map(|&v| v as f64).collect(); let mut c = [0.0f64; M * M]; - for i in 0..M { - for j in 0..M { - let mut s = 0.0f64; - for k in 0..M { - s += a[i * M + k] as f64 * b[k * M + j] as f64; - } - c[i * M + j] = s; - } - } + ndarray::simd::gemm_f64_tiled(M, M, M, 1.0, &af, M, &bf, M, 0.0, &mut c, M); c } diff --git a/examples/subpel_tap_tile.rs b/examples/subpel_tap_tile.rs index a11e070c..b1edcf37 100644 --- a/examples/subpel_tap_tile.rs +++ b/examples/subpel_tap_tile.rs @@ -93,17 +93,14 @@ mod amx { out } + /// Direct `f64` ground truth via the crate-native bit-exact tiled GEMM + /// (`ndarray::simd::gemm_f64_tiled`) — bit-identical to the naive + /// triple loop at α=1 β=0, and entirely in-crate (no external matmul). fn direct_matmul(a: &[f32; M * M], b: &[f32; M * M]) -> [f64; M * M] { + let af: Vec = a.iter().map(|&v| v as f64).collect(); + let bf: Vec = b.iter().map(|&v| v as f64).collect(); let mut c = [0.0f64; M * M]; - for i in 0..M { - for j in 0..M { - let mut s = 0.0f64; - for k in 0..M { - s += a[i * M + k] as f64 * b[k * M + j] as f64; - } - c[i * M + j] = s; - } - } + ndarray::simd::gemm_f64_tiled(M, M, M, 1.0, &af, M, &bf, M, 0.0, &mut c, M); c } diff --git a/src/backend/native.rs b/src/backend/native.rs index 69b25e77..836297da 100644 --- a/src/backend/native.rs +++ b/src/backend/native.rs @@ -468,52 +468,10 @@ mod scalar { } } - /// Tiled GEMM (f64, scalar reference) - #[allow(dead_code)] - pub fn gemm_f64_tiled( - m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, - c: &mut [f64], ldc: usize, - ) { - const TILE: usize = 64; - - if beta == 0.0 { - for i in 0..m { - for j in 0..n { - c[i * ldc + j] = 0.0; - } - } - } else if beta != 1.0 { - for i in 0..m { - for j in 0..n { - c[i * ldc + j] *= beta; - } - } - } - - let mut kk = 0; - while kk < k { - let kb = TILE.min(k - kk); - let mut ii = 0; - while ii < m { - let ib = TILE.min(m - ii); - let mut jj = 0; - while jj < n { - let jb = TILE.min(n - jj); - for i in 0..ib { - for p in 0..kb { - let a_val = alpha * a[(ii + i) * lda + (kk + p)]; - for j in 0..jb { - c[(ii + i) * ldc + (jj + j)] += a_val * b[(kk + p) * ldb + (jj + j)]; - } - } - } - jj += jb; - } - ii += ib; - } - kk += kb; - } - } + // The f64 sibling (`gemm_f64_tiled`) graduated to `src/simd_ops.rs`, + // public as `ndarray::simd::gemm_f64_tiled` — same tiled algorithm, + // innermost loop F64x8-vectorized, with a documented bit-exactness + // contract (unfused mul+add; bit-identical to this scalar shape). pub fn gemv_f32(m: usize, n: usize, alpha: f32, a: &[f32], lda: usize, x: &[f32], beta: f32, y: &mut [f32]) { for i in 0..m { diff --git a/src/simd.rs b/src/simd.rs index 4bec0d94..552ce33d 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -566,6 +566,14 @@ pub use crate::hpc::fingerprint::{ // the reason the JIT-native option was deemed unnecessary. See the // "Foundation primitives — do not remove" notice in `src/simd_ops.rs`. pub use crate::simd_ops::{array_chunks, array_chunks_checked, array_windows, array_windows_checked}; +// Crate-native tiled f64 GEMM (`C := α·A·B + β·C`) with a bit-exactness +// contract: unfused mul+add in ascending-k order per element → bit-identical +// on every backend (AVX-512/AVX2/NEON/WASM/scalar) and, at α=1 β=0, +// bit-identical to the naive triple-loop reference. This is the in-crate +// ground-truth GEMM for probes/certification — `backend::native::gemm_f64` +// delegates to the external `matrixmultiply` crate instead. Alloc-free, so +// available without `std` (hence outside the std-gated kernel list below). +pub use crate::simd_ops::gemm_f64_tiled; pub use crate::simd_soa::MultiLaneColumn; pub use crate::hpc::quantized::{ diff --git a/src/simd_ops.rs b/src/simd_ops.rs index 455f915e..d1cff84b 100644 --- a/src/simd_ops.rs +++ b/src/simd_ops.rs @@ -880,3 +880,338 @@ mod add_mul_tests { } } } + +// ═══════════════════════════════════════════════════════════════════ +// Tiled f64 GEMM — the crate-native, bit-exact reference GEMM +// ═══════════════════════════════════════════════════════════════════ + +/// Tiled GEMM: `C := alpha·A·B + beta·C` (f64, row-major, leading dims). +/// +/// The crate-native f64 GEMM: a TILE=64-blocked loop nest with the +/// innermost j-run vectorized on the dispatched [`F64x8`] (AVX-512 / +/// AVX2 / NEON / WASM-SIMD128 / scalar — one source, every backend). +/// Unlike `backend::native::gemm_f64`, which delegates to the external +/// `matrixmultiply` crate, this kernel is entirely in-crate and carries +/// a bit-exactness contract, so probes can certify other kernels +/// against it as ground truth. +/// +/// # Bit-exactness contract +/// +/// Every `C[i,j]` receives the exact IEEE op sequence +/// `c = c + (alpha·A[i,p])·B[p,j]` in ascending-`p` order, with the +/// multiply and add **unfused** (the `*`/`+` operators on the SIMD +/// types lower to plain `mul`/`add` on every backend — never +/// FMA-contracted). The result is therefore **bit-identical across all +/// backends**, and for `alpha == 1.0, beta == 0.0` (the certification +/// shape) bit-identical to the naive triple-loop reference +/// `s = Σ_p A[i,p]·B[p,j]; C[i,j] = s`. +/// +/// Edge semantics: `m == 0` or `n == 0` is a no-op (`c` untouched); +/// `k == 0` applies only the beta pass; `beta == 0.0` overwrites `c` +/// without reading it (BLAS convention — pre-existing NaN/Inf in `c` +/// do not propagate); `beta == 1.0` accumulates into `c` as-is; +/// non-finite values in `a`/`b` propagate per IEEE. +/// +/// # Panics +/// +/// Panics if `lda < k`, `ldb < n`, `ldc < n`, or if `a`, `b`, or `c` +/// is shorter than `(rows − 1)·ld + cols` for its shape. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::gemm_f64_tiled; +/// +/// let a = [1.0, 2.0, 3.0, 4.0]; // 2×2 row-major +/// let b = [5.0, 6.0, 7.0, 8.0]; // 2×2 row-major +/// let mut c = [0.0f64; 4]; +/// gemm_f64_tiled(2, 2, 2, 1.0, &a, 2, &b, 2, 0.0, &mut c, 2); +/// assert_eq!(c, [19.0, 22.0, 43.0, 50.0]); +/// ``` +// BLAS gemm argument list kept verbatim (drop-in for `backend::native::gemm_f64`). +#[allow(clippy::too_many_arguments)] +pub fn gemm_f64_tiled( + m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, c: &mut [f64], + ldc: usize, +) { + if m == 0 || n == 0 { + return; + } + assert!(lda >= k, "lda ({lda}) must be >= k ({k})"); + assert!(ldb >= n, "ldb ({ldb}) must be >= n ({n})"); + assert!(ldc >= n, "ldc ({ldc}) must be >= n ({n})"); + if k > 0 { + assert!(a.len() >= (m - 1) * lda + k, "a too short for m×k with lda"); + assert!(b.len() >= (k - 1) * ldb + n, "b too short for k×n with ldb"); + } + assert!(c.len() >= (m - 1) * ldc + n, "c too short for m×n with ldc"); + + const TILE: usize = 64; + const LANES: usize = 8; // F64x8 + + // Beta pass — one scalar op per element (`= 0.0` / `*= beta`), so the + // pass is trivially bit-identical however it is traversed. + if beta == 0.0 { + for i in 0..m { + for x in &mut c[i * ldc..i * ldc + n] { + *x = 0.0; + } + } + } else if beta != 1.0 { + for i in 0..m { + for x in &mut c[i * ldc..i * ldc + n] { + *x *= beta; + } + } + } + + let mut kk = 0; + while kk < k { + let kb = TILE.min(k - kk); + let mut ii = 0; + while ii < m { + let ib = TILE.min(m - ii); + let mut jj = 0; + while jj < n { + let jb = TILE.min(n - jj); + for i in 0..ib { + for p in 0..kb { + let a_val = alpha * a[(ii + i) * lda + (kk + p)]; + let va = F64x8::splat(a_val); + let b_row = (kk + p) * ldb + jj; + let c_row = (ii + i) * ldc + jj; + let mut j = 0; + // Unfused `vc + va·vb` — bit-identical to the scalar tail. + while j + LANES <= jb { + let vb = F64x8::from_slice(&b[b_row + j..]); + let vc = F64x8::from_slice(&c[c_row + j..]); + (vc + va * vb).copy_to_slice(&mut c[c_row + j..]); + j += LANES; + } + while j < jb { + c[c_row + j] += a_val * b[b_row + j]; + j += 1; + } + } + } + jj += jb; + } + ii += ib; + } + kk += kb; + } +} + +#[cfg(test)] +mod gemm_f64_tiled_tests { + use super::gemm_f64_tiled; + + /// splitmix64 finalizer → f64 in [-1, 1) with full mantissa churn. + fn mix(mut z: u64) -> u64 { + z = z.wrapping_add(0x9E37_79B9_7F4A_7C15); + z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9); + z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB); + z ^ (z >> 31) + } + fn rnd(seed: u64) -> f64 { + (mix(seed) as f64 / u64::MAX as f64) * 2.0 - 1.0 + } + + /// Reference with the SAME documented op sequence: beta pass, then + /// `c += (alpha·a)·b` in ascending-p order. Bit-exact oracle. + #[allow(clippy::too_many_arguments)] + fn ref_gemm_ipj( + m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, + c: &mut [f64], ldc: usize, + ) { + for i in 0..m { + for x in &mut c[i * ldc..i * ldc + n] { + if beta == 0.0 { + *x = 0.0; + } else if beta != 1.0 { + *x *= beta; + } + } + } + for i in 0..m { + for p in 0..k { + let a_val = alpha * a[i * lda + p]; + for j in 0..n { + c[i * ldc + j] += a_val * b[p * ldb + j]; + } + } + } + } + + /// The classic naive form (local accumulator, i-j-k). Bit-identical to + /// the kernel for alpha=1, beta=0 — the certification shape. + #[allow(clippy::too_many_arguments)] + fn ref_naive_ijk( + m: usize, n: usize, k: usize, a: &[f64], lda: usize, b: &[f64], ldb: usize, ldc: usize, + ) -> Vec { + let mut c = vec![0.0f64; if m == 0 { 0 } else { (m - 1) * ldc + n }]; + for i in 0..m { + for j in 0..n { + let mut s = 0.0f64; + for p in 0..k { + s += a[i * lda + p] * b[p * ldb + j]; + } + c[i * ldc + j] = s; + } + } + c + } + + fn fill(len: usize, salt: u64) -> Vec { + (0..len).map(|i| rnd(i as u64 ^ salt)).collect() + } + + fn assert_bits_eq(got: &[f64], want: &[f64], tag: &str) { + assert_eq!(got.len(), want.len(), "{tag}: length"); + for (i, (g, w)) in got.iter().zip(want.iter()).enumerate() { + assert_eq!(g.to_bits(), w.to_bits(), "{tag}: bit mismatch at {i}: got {g:e}, want {w:e}"); + } + } + + const SHAPES: &[(usize, usize, usize)] = &[ + (1, 1, 1), + (2, 3, 4), + (3, 7, 5), + (7, 8, 9), + (8, 8, 8), + (16, 16, 16), + (9, 17, 33), + (16, 16, 64), + (5, 8, 128), + (65, 3, 65), + (64, 64, 64), + (70, 70, 70), + (1, 100, 1), + ]; + + #[test] + fn bit_exact_vs_both_references_shape_sweep() { + for &(m, n, k) in SHAPES { + let a = fill(m * k, 0xA); + let b = fill(k * n, 0xB); + let mut c = vec![0.0f64; m * n]; + gemm_f64_tiled(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c, n); + + let mut c_ref = vec![0.0f64; m * n]; + ref_gemm_ipj(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c_ref, n); + assert_bits_eq(&c, &c_ref, &format!("ipj ({m},{n},{k})")); + + let c_naive = ref_naive_ijk(m, n, k, &a, k, &b, n, n); + assert_bits_eq(&c, &c_naive, &format!("naive ({m},{n},{k})")); + } + } + + #[test] + fn bit_exact_general_alpha_beta_preloaded_c() { + for &(m, n, k) in &[(16, 16, 16), (9, 17, 33), (70, 70, 70)] { + for &(alpha, beta) in &[(0.75, 2.5), (-1.25, 1.0), (2.0, -0.5), (0.0, 3.0)] { + let a = fill(m * k, 0x1); + let b = fill(k * n, 0x2); + let c0 = fill(m * n, 0x3); + let mut c = c0.clone(); + gemm_f64_tiled(m, n, k, alpha, &a, k, &b, n, beta, &mut c, n); + let mut c_ref = c0.clone(); + ref_gemm_ipj(m, n, k, alpha, &a, k, &b, n, beta, &mut c_ref, n); + assert_bits_eq(&c, &c_ref, &format!("({m},{n},{k}) α={alpha} β={beta}")); + } + } + } + + #[test] + fn strided_leading_dims_padding_untouched() { + let (m, n, k) = (17usize, 9usize, 33usize); + let (lda, ldb, ldc) = (k + 5, n + 3, n + 7); + const SENTINEL: f64 = 12345.678; + let mut a = vec![SENTINEL; (m - 1) * lda + k + 4]; + let mut b = vec![SENTINEL; (k - 1) * ldb + n + 2]; + let mut c = vec![SENTINEL; (m - 1) * ldc + n + 6]; + for i in 0..m { + for p in 0..k { + a[i * lda + p] = rnd((i * k + p) as u64 ^ 0x11); + } + } + for p in 0..k { + for j in 0..n { + b[p * ldb + j] = rnd((p * n + j) as u64 ^ 0x22); + } + } + for i in 0..m { + for j in 0..n { + c[i * ldc + j] = rnd((i * n + j) as u64 ^ 0x33); + } + } + let mut c_ref = c.clone(); + gemm_f64_tiled(m, n, k, 1.5, &a, lda, &b, ldb, 0.5, &mut c, ldc); + ref_gemm_ipj(m, n, k, 1.5, &a, lda, &b, ldb, 0.5, &mut c_ref, ldc); + assert_bits_eq(&c, &c_ref, "strided"); + // Row padding (columns n..ldc) must be untouched. + for i in 0..m - 1 { + for j in n..ldc { + assert_eq!(c[i * ldc + j], c_ref[i * ldc + j], "padding clobbered at row {i} col {j}"); + assert_eq!(c[i * ldc + j].to_bits(), SENTINEL.to_bits(), "padding not sentinel at row {i} col {j}"); + } + } + } + + #[test] + fn beta_zero_overwrites_nan_c() { + let (m, n, k) = (8usize, 11usize, 5usize); + let a = fill(m * k, 0x4); + let b = fill(k * n, 0x5); + let mut c = vec![f64::NAN; m * n]; + gemm_f64_tiled(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c, n); + let c_naive = ref_naive_ijk(m, n, k, &a, k, &b, n, n); + assert_bits_eq(&c, &c_naive, "beta=0 over NaN"); + } + + #[test] + fn k_zero_applies_beta_only() { + let (m, n) = (6usize, 10usize); + let c0 = fill(m * n, 0x6); + // beta = 0 → zeros + let mut c = c0.clone(); + gemm_f64_tiled(m, n, 0, 1.0, &[], 0, &[], n, 0.0, &mut c, n); + assert!(c.iter().all(|x| x.to_bits() == 0.0f64.to_bits()), "beta=0, k=0 must zero c"); + // beta = 2.5 → scaled + let mut c = c0.clone(); + gemm_f64_tiled(m, n, 0, 1.0, &[], 0, &[], n, 2.5, &mut c, n); + for (got, want) in c.iter().zip(c0.iter()) { + assert_eq!(got.to_bits(), (want * 2.5).to_bits(), "beta=2.5, k=0 must scale c"); + } + // beta = 1 → untouched + let mut c = c0.clone(); + gemm_f64_tiled(m, n, 0, 1.0, &[], 0, &[], n, 1.0, &mut c, n); + assert_bits_eq(&c, &c0, "beta=1, k=0"); + } + + #[test] + fn m_or_n_zero_is_noop() { + let mut c = vec![7.5f64; 12]; + gemm_f64_tiled(0, 4, 3, 1.0, &[], 3, &[1.0; 12], 4, 0.0, &mut c, 4); + gemm_f64_tiled(3, 0, 4, 1.0, &[1.0; 12], 4, &[], 1, 0.0, &mut c, 1); + assert!(c.iter().all(|x| x.to_bits() == 7.5f64.to_bits()), "m==0 / n==0 must not touch c"); + } + + #[test] + fn denormals_and_negzero_bit_exact() { + let (m, n, k) = (4usize, 9usize, 7usize); + let mut a = fill(m * k, 0x7); + let mut b = fill(k * n, 0x8); + // Inject denormals, -0.0, and tiny/huge magnitudes. + a[0] = 5e-324; // smallest denormal + a[1] = -0.0; + a[2] = 1e-310; + b[0] = -5e-324; + b[1] = -0.0; + b[2] = 1e300; + let mut c = vec![0.0f64; m * n]; + gemm_f64_tiled(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c, n); + let c_naive = ref_naive_ijk(m, n, k, &a, k, &b, n, n); + assert_bits_eq(&c, &c_naive, "denormals"); + } +} From c880a5b881de4cca3c7850e7c95bb1779cacb556 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 21:25:16 +0000 Subject: [PATCH 4/7] simd: fold in adversarial-review findings on gemm_f64_tiled MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three-angle verify (IEEE semantics / W1a compliance / regression trace) returned PASS on all angles; this folds in the surviving findings: - P1: correct the simd.rs re-export comment — pub mod simd/simd_ops are std-gated in lib.rs, so the kernel (though alloc-free) is reachable only in std builds; the prior comment claimed no_std availability. - Scope the cross-backend bit-identity claim to non-NaN inputs (NaN payload propagation is backend-defined; WASM may canonicalize). - Document alpha == 0.0 non-short-circuit semantics (a/b still read, 0*Inf = NaN propagates, -0.0 in c can flip — unlike BLAS quick-return). - Panics doc: state that all checks are skipped at m/n == 0 and the a/b length checks at k == 0. - Overflow-check the length-extent asserts (checked_mul) so a wrapping (rows-1)*ld cannot spuriously pass the precondition in release. - Widen the parity corpus to 70+ kernel invocations (full 13-shape sweep x 4 alpha/beta combos) to meet W1a criterion 3's 50+ letter. - Cite the free-fn GEMM-family precedent in the doc-comment for the simd-savant pre-merge audit. Gates re-run: clippy -D warnings clean, 7/7 gemm_f64_tiled tests, doctest green. --- .claude/blackboard.md | 17 +++++++++++++++++ src/simd.rs | 5 +++-- src/simd_ops.rs | 43 ++++++++++++++++++++++++++++++++++--------- 3 files changed, 54 insertions(+), 11 deletions(-) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index ec4212f9..c0c4f64b 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -59,6 +59,23 @@ queued (Gotcha-14 assert guard, PackedBf16B throughput leg, positive- operand comment fix). Adversarial 3-angle verify workflow on this diff was in flight at commit time; findings (if any) land as follow-up. +**[VERIFY OUTCOME, same day]** 3-angle adversarial review (IEEE/bit- +exactness, W1a compliance, regression-trace): **PASS / PASS / PASS**. +One convergent P1 fixed in the follow-up commit: the simd.rs re-export +comment claimed no_std availability, but `pub mod simd`/`simd_ops` are +std-gated in lib.rs — the no_std build passes because the code is +compiled OUT (comment corrected; un-gating simd_ops for a genuinely +no_std kernel is possible future work, not this diff). P2s folded in: +cross-backend bit-identity scoped to non-NaN inputs (NaN payloads are +backend-defined, WASM may canonicalize); `alpha == 0.0` non-short- +circuit documented (0·Inf=NaN propagates, −0.0 can flip, unlike BLAS +quick-return); Panics doc states which checks are skipped (m/n==0, k==0); +length-extent asserts now overflow-checked (`checked_mul`); parity corpus +widened to 70+ invocations (full 13-shape sweep × 4 α/β combos) for the +W1a "50+" letter; free-fn-shape precedent sentence added to the doc +(GEMM family: bf16_tile_gemm_16x16, matmul_i8_to_i32). Deferred as +noise: x87-only i586 excess-precision footnote (tier-2, pre-SSE2). + ## 2026-07-06 — Review pass: health check green + 10 findings on subpel_tap_tile (#235) Review-only session (no code changes). **Health check:** `cargo fmt --check` diff --git a/src/simd.rs b/src/simd.rs index 552ce33d..59f868ca 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -571,8 +571,9 @@ pub use crate::simd_ops::{array_chunks, array_chunks_checked, array_windows, arr // on every backend (AVX-512/AVX2/NEON/WASM/scalar) and, at α=1 β=0, // bit-identical to the naive triple-loop reference. This is the in-crate // ground-truth GEMM for probes/certification — `backend::native::gemm_f64` -// delegates to the external `matrixmultiply` crate instead. Alloc-free, so -// available without `std` (hence outside the std-gated kernel list below). +// delegates to the external `matrixmultiply` crate instead. The kernel is +// alloc-free, but `pub mod simd`/`simd_ops` are std-gated in lib.rs, so +// like every kernel here it is reachable only in `std` builds today. pub use crate::simd_ops::gemm_f64_tiled; pub use crate::simd_soa::MultiLaneColumn; diff --git a/src/simd_ops.rs b/src/simd_ops.rs index d1cff84b..39291810 100644 --- a/src/simd_ops.rs +++ b/src/simd_ops.rs @@ -902,20 +902,33 @@ mod add_mul_tests { /// multiply and add **unfused** (the `*`/`+` operators on the SIMD /// types lower to plain `mul`/`add` on every backend — never /// FMA-contracted). The result is therefore **bit-identical across all -/// backends**, and for `alpha == 1.0, beta == 0.0` (the certification -/// shape) bit-identical to the naive triple-loop reference -/// `s = Σ_p A[i,p]·B[p,j]; C[i,j] = s`. +/// backends for non-NaN inputs** (NaN-*ness* propagates everywhere, +/// but NaN *payloads* are backend-defined — WASM in particular may +/// canonicalize), and for `alpha == 1.0, beta == 0.0` (the +/// certification shape) bit-identical to the naive triple-loop +/// reference `s = Σ_p A[i,p]·B[p,j]; C[i,j] = s`. /// /// Edge semantics: `m == 0` or `n == 0` is a no-op (`c` untouched); /// `k == 0` applies only the beta pass; `beta == 0.0` overwrites `c` /// without reading it (BLAS convention — pre-existing NaN/Inf in `c` /// do not propagate); `beta == 1.0` accumulates into `c` as-is; -/// non-finite values in `a`/`b` propagate per IEEE. +/// non-finite values in `a`/`b` propagate per IEEE. Unlike BLAS +/// `xGEMM`, `alpha == 0.0` does **not** short-circuit the product +/// term: `a`/`b` are still read and `(±0.0·a)·b` terms are added, so +/// NaN/Inf in `a`/`b` propagate into `c` (`0·Inf = NaN`) and a `-0.0` +/// in `c` can flip to `+0.0`. +/// +/// Free-function shape follows the simd-surface GEMM family precedent +/// ([`bf16_tile_gemm_16x16`], `matmul_i8_to_i32`) — an 11-argument +/// slice-level BLAS kernel has no typed-wrapper home (W1a). /// /// # Panics /// /// Panics if `lda < k`, `ldb < n`, `ldc < n`, or if `a`, `b`, or `c` -/// is shorter than `(rows − 1)·ld + cols` for its shape. +/// is shorter than `(rows − 1)·ld + cols` for its shape (computed +/// overflow-checked). All checks are skipped when `m == 0 || n == 0` +/// (the no-op path), and the `a`/`b` length checks are skipped when +/// `k == 0` (neither is read). /// /// # Examples /// @@ -940,11 +953,20 @@ pub fn gemm_f64_tiled( assert!(lda >= k, "lda ({lda}) must be >= k ({k})"); assert!(ldb >= n, "ldb ({ldb}) must be >= n ({n})"); assert!(ldc >= n, "ldc ({ldc}) must be >= n ({n})"); + // Overflow-checked extents: a wrapping (rows−1)·ld in release would let + // the precondition spuriously pass and surface later as a mid-compute + // index panic (after the beta pass mutated `c`). + let extent = |rows: usize, ld: usize, cols: usize| { + (rows - 1) + .checked_mul(ld) + .and_then(|x| x.checked_add(cols)) + .expect("matrix extent overflows usize") + }; if k > 0 { - assert!(a.len() >= (m - 1) * lda + k, "a too short for m×k with lda"); - assert!(b.len() >= (k - 1) * ldb + n, "b too short for k×n with ldb"); + assert!(a.len() >= extent(m, lda, k), "a too short for m×k with lda"); + assert!(b.len() >= extent(k, ldb, n), "b too short for k×n with ldb"); } - assert!(c.len() >= (m - 1) * ldc + n, "c too short for m×n with ldc"); + assert!(c.len() >= extent(m, ldc, n), "c too short for m×n with ldc"); const TILE: usize = 64; const LANES: usize = 8; // F64x8 @@ -1108,7 +1130,10 @@ mod gemm_f64_tiled_tests { #[test] fn bit_exact_general_alpha_beta_preloaded_c() { - for &(m, n, k) in &[(16, 16, 16), (9, 17, 33), (70, 70, 70)] { + // Full shape sweep × 4 (alpha, beta) combos = 52 kernel invocations, + // + 13 in the sweep above + edge tests → 70+ corpus entries (W1a + // criterion 3's "50+ test inputs"). + for &(m, n, k) in SHAPES { for &(alpha, beta) in &[(0.75, 2.5), (-1.25, 1.0), (2.0, -0.5), (0.0, 3.0)] { let a = fill(m * k, 0x1); let b = fill(k * n, 0x2); From 666403d6bd7cfe40ee87ed6fe3e13dc23e10ddd6 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 21:28:42 +0000 Subject: [PATCH 5/7] blackboard: correct AMX-flip note with measured amx_probe diagnosis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The earlier 'XTILEDATA enablement drift' guess was wrong. amx_probe shows CPUID TILE/INT8/BF16 all false and cpu_model() = OtherX86 (morning run: EMR 0xCF with AMX) — the session container was rescheduled onto non-AMX silicon. Not Gotcha 4 (has_amx()==true/available==false signature), not Gotcha 14 (corruption while available). Gotcha 9 tier-printing surfaced it. Rule recorded: re-probe before any AMX-tier claim in remote sessions; never carry amx_available() across runs. --- .claude/blackboard.md | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index c0c4f64b..e09e4ed7 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -48,11 +48,20 @@ std-gated in lib.rs — see reviewer note below if that changed). palette_codec.rs:806 needless_range_loop) — not touched here; the house gate (`cargo clippy -- -D warnings`, no --tests) is clean. -[OBSERVED] `amx_available()` flipped true→false on this VM between runs -two hours apart (subpel ran AMX TDPBF16PS earlier, F32x16 polyfill now; -identical BF16-class errors either way). Host/hypervisor XTILEDATA -enablement drift — Gotcha 9/14 adjacent; worth a note in AMX_GOTCHAS if -observed again. +[OBSERVED, then MEASURED] `amx_available()` flipped true→false between +runs two hours apart (subpel ran AMX TDPBF16PS earlier, F32x16 polyfill +later; identical BF16-class errors either way — tier ladder correct). +**Diagnosis via `examples/amx_probe` per AMX_GOTCHAS discipline** (initial +"enablement drift / Gotcha 14 adjacent" guess was WRONG): CPUID leaf-7 +TILE/INT8/BF16 bits all false, `cpu_model() = OtherX86`, `has_amx() = +false` — the **silicon identity itself changed** (morning: EMR 0xCF with +AMX). The session container was rescheduled onto non-AMX/CPUID-masked +silicon. NOT Gotcha 4 (that signature is has_amx()==true with +available==false) and NOT Gotcha 14 (corruption while available==true). +Gotcha 9's always-print-the-tier discipline is what surfaced the flip. +Consequence for remote sessions: the host under an ephemeral container +can change mid-session — re-run `amx_probe`/`amx_report()` before any +AMX-tier claim, never carry `amx_available()` results across runs. [LOOSE END] subpel_tap_tile findings 1-3 from the same-day review still queued (Gotcha-14 assert guard, PackedBf16B throughput leg, positive- From ec0373f19c6e39a30e567c55a494214fccf9c773 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 21:44:04 +0000 Subject: [PATCH 6/7] =?UTF-8?q?simd:=20complete=20the=20f64=20GEMM=20?= =?UTF-8?q?=E2=80=94=20FMA=20tier,=20register-resident=20C,=20own-engine?= =?UTF-8?q?=20swap?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Completes the f64 GEMM family on the canonical surface: 1. gemm_f64_tiled_fma — fast fused tier via const-generic gemm_f64_tiled_impl (monomorphized). Same tiling and ascending-p per-element order; each step is one fused rounding (F64x8::mul_add / f64::mul_add). Bit-identical to the reference tier on integer-valued operands (asserted across the full shape sweep); last-ulp-per-step on general floats (tolerance scaled k*eps to the summand magnitude). Cross-backend fusion caveats documented. 2. Register-resident C row-block (both tiers): the C block is loaded once into [F64x8; TILE/LANES] accumulators plus a scalar tail array, the whole kb-loop accumulates in registers, one store. f64 store/ reload never rounds, so the per-element op sequence — and every bit-exactness contract and test — is unchanged. Measured 2.2x: ref 4.6 -> 10.0 GFLOP/s, fma 4.7 -> 10.7 GFLOP/s. 3. backend::native::gemm_f64 engine swapped from matrixmultiply::dgemm to gemm_f64_tiled_fma: the f64 GEMM behind BlasFloat::backend_gemm, hpc::blas_level3::blas_gemm, and batched linalg is now entirely in-crate Rust. matrixmultiply remains for gemm_f32 and upstream Array::dot (untouched). Measured own-engine gap ~3.1x (was 6.6x before the restructure): 1024^3 fma 10.3 vs matrixmultiply 33.7 GFLOP/s, max |fma - matrixmultiply| <= 1.4e-13. Next rung recorded: B-panel register reuse + packing. Gates: clippy -D warnings clean, 2185/2185 lib tests green with the engine swap, doctests green, bit-equality suite green after the restructure. --- .claude/blackboard.md | 43 +++++++++++ src/backend/native.rs | 30 +++----- src/simd.rs | 5 +- src/simd_ops.rs | 161 +++++++++++++++++++++++++++++++++++++++--- 4 files changed, 205 insertions(+), 34 deletions(-) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index e09e4ed7..5a3b649b 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3,6 +3,49 @@ > **Read this first.** The "Polyglot Notebook" architecture below is a > separate/older program, not the current epoch. +## 2026-07-06 (latest) — F64 GEMM completed: FMA tier + register residency + native-engine swap + +Operator: "complete the F64 gemm… and pr". Three moves, all on the entry +below's foundation: + +1. **`gemm_f64_tiled_fma`** — fast fused tier via const-generic + `gemm_f64_tiled_impl` (monomorphized, no runtime + branch). Per-element ascending-p order preserved; fused step + `c = fma(α·a, b, c)`. Bit-identical to the reference tier on + integer-valued operands (products+sums < 2^53 — asserted with + `assert_eq!` on the full shape sweep); last-ulp-per-step differences + on general floats (tolerance test scaled k·ε to SUMMAND magnitude — + the initial result-scaled tolerance was wrong for cancellation-heavy + elements and failed honestly). Cross-backend caveat documented (WASM + without relaxed-simd has unfused vector lanes + fused scalar tail). +2. **Register-resident C row-block** (both tiers): C block loaded into + `[F64x8; TILE/LANES]` accumulators + scalar tail array ONCE per + (kk,ii,jj,i), whole kb-loop accumulates in registers, one store. + f64→f64 store/reload never rounds ⇒ per-element op sequence + unchanged ⇒ every bit-equality test green untouched. [MEASURED] + ref 4.6→10.0 GF, fma 4.7→10.7 GF (2.2×). +3. **Engine swap:** `backend::native::gemm_f64` now routes to + `gemm_f64_tiled_fma` — the f64 GEMM behind `BlasFloat::backend_gemm` + / `hpc::blas_level3::blas_gemm` / batched linalg is entirely own + Rust; matrixmultiply remains only in gemm_f32 and upstream + `Array::dot` (impl_linalg.rs, untouched at ~33 GF). + +[MEASURED, 3-engine, this VM (v3 compile → AVX2 arm, PREFERRED_F64_LANES=4; +host runtime has avx512f but committed .cargo config is v3)]: +256³/512³/1024³ — ref 11.2/10.0/9.6 GF | fma 11.9/10.7/10.3 GF | +matrixmultiply 34.1/32.3/33.7 GF | max|fma−mm| ≤ 1.4e-13. +**Own-engine gap: ~3.1×** (was 6.6× pre-restructure). Trade accepted per +operator priority (own reverse-engineered Rust in the path, auditable +numerics); blast radius = hpc BLAS surface only. 2185/2185 lib tests +green WITH the swap. + +[LOOSE END → next rung] Closing the 3× needs a real microkernel: B-panel +register reuse (i-tiling IR=2..4 × narrower j-block), then A/B packing — +the matrixmultiply Goto recipe, own-Rust edition. Also: `gemm_f32_tiled` +still dead in native.rs `mod scalar` (f32 sibling completion); +avx512f compile arm untested on CI (v3 config) — the F64x8=__m512d arm +runs only on local v4 builds. + ## 2026-07-06 (later) — `ndarray::simd::gemm_f64_tiled` surfaced (operator directive) The crate-native tiled f64 GEMM graduated from dead code diff --git a/src/backend/native.rs b/src/backend/native.rs index 836297da..df3c1d91 100644 --- a/src/backend/native.rs +++ b/src/backend/native.rs @@ -237,32 +237,18 @@ pub fn gemm_f32( } /// GEMM: C = alpha * A * B + beta * C (f64, row-major). +/// +/// Engine: the crate-native `simd_ops::gemm_f64_tiled_fma` (F64x8-vectorized +/// fused tier — entirely in-crate Rust, no external matmul engine). The +/// unfused sibling `ndarray::simd::gemm_f64_tiled` is the bit-exact +/// certification reference; this entry favors throughput. Swapped from +/// `matrixmultiply::dgemm` — results agree to BLAS tolerance (last-ulp +/// accumulation-order/fusion differences), not bit-for-bit. pub fn gemm_f64( m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, c: &mut [f64], ldc: usize, ) { - if m == 0 || n == 0 { - return; - } - // SAFETY: same as sgemm — valid slices, row-major strides. - unsafe { - matrixmultiply::dgemm( - m, - k, - n, - alpha, - a.as_ptr(), - lda as isize, - 1, - b.as_ptr(), - ldb as isize, - 1, - beta, - c.as_mut_ptr(), - ldc as isize, - 1, - ); - } + crate::simd_ops::gemm_f64_tiled_fma(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } // ─── GEMV dispatch ─────────────────────────────────────────────── diff --git a/src/simd.rs b/src/simd.rs index 59f868ca..66ce3fa4 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -574,7 +574,10 @@ pub use crate::simd_ops::{array_chunks, array_chunks_checked, array_windows, arr // delegates to the external `matrixmultiply` crate instead. The kernel is // alloc-free, but `pub mod simd`/`simd_ops` are std-gated in lib.rs, so // like every kernel here it is reachable only in `std` builds today. -pub use crate::simd_ops::gemm_f64_tiled; +// `gemm_f64_tiled_fma` is the fast fused tier (same tiling/order, one +// rounding per step) — the engine behind `backend::native::gemm_f64`; +// the unfused reference tier stays the certification ground truth. +pub use crate::simd_ops::{gemm_f64_tiled, gemm_f64_tiled_fma}; pub use crate::simd_soa::MultiLaneColumn; pub use crate::hpc::quantized::{ diff --git a/src/simd_ops.rs b/src/simd_ops.rs index 39291810..fcf07e65 100644 --- a/src/simd_ops.rs +++ b/src/simd_ops.rs @@ -946,6 +946,59 @@ mod add_mul_tests { pub fn gemm_f64_tiled( m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, c: &mut [f64], ldc: usize, +) { + gemm_f64_tiled_impl::(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); +} + +/// Fast tier of [`gemm_f64_tiled`]: same tiling, same ascending-`p` +/// per-element order, but each step is a **fused** multiply-add +/// (`c = fma(alpha·A[i,p], B[p,j], c)` — one rounding per step via +/// `F64x8::mul_add`: `vfmadd`/`vfmaq_f64` on x86/NEON, `f64::mul_add` +/// on the scalar tail and polyfill backends). +/// +/// Semantics relative to the reference tier: +/// - **Deterministic per build**, but NOT bit-stable across backends: +/// fusion quality differs (WASM without `relaxed-simd` has no fused +/// vector `mul_add`, so its vector lanes match the unfused reference +/// while the scalar tail is still fused). +/// - **Bit-identical to [`gemm_f64_tiled`]** whenever every +/// `a_val·b` product and running sum is exactly representable — +/// e.g. integer-valued data with products and accumulation inside +/// `2^53` (asserted in tests). On general float data each step +/// differs from the reference by at most the dropped intermediate +/// rounding (classic FMA-vs-mul+add last-ulp effects). +/// - Same argument contract, edge semantics, and panics as +/// [`gemm_f64_tiled`]. +/// +/// Use the reference tier for certification/ground truth; use this +/// tier when GEMM throughput matters. +/// +/// # Examples +/// +/// ``` +/// use ndarray::simd::{gemm_f64_tiled, gemm_f64_tiled_fma}; +/// +/// let a = [1.0, 2.0, 3.0, 4.0]; // integer-valued → tiers bit-identical +/// let b = [5.0, 6.0, 7.0, 8.0]; +/// let (mut c1, mut c2) = ([0.0f64; 4], [0.0f64; 4]); +/// gemm_f64_tiled(2, 2, 2, 1.0, &a, 2, &b, 2, 0.0, &mut c1, 2); +/// gemm_f64_tiled_fma(2, 2, 2, 1.0, &a, 2, &b, 2, 0.0, &mut c2, 2); +/// assert_eq!(c1, c2); +/// assert_eq!(c1, [19.0, 22.0, 43.0, 50.0]); +/// ``` +#[allow(clippy::too_many_arguments)] +pub fn gemm_f64_tiled_fma( + m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, c: &mut [f64], + ldc: usize, +) { + gemm_f64_tiled_impl::(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); +} + +#[allow(clippy::too_many_arguments)] +#[inline(always)] +fn gemm_f64_tiled_impl( + m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, c: &mut [f64], + ldc: usize, ) { if m == 0 || n == 0 { return; @@ -996,25 +1049,45 @@ pub fn gemm_f64_tiled( let mut jj = 0; while jj < n { let jb = TILE.min(n - jj); + let nvec = jb / LANES; // ≤ TILE/LANES full-lane groups + let ntail = jb % LANES; for i in 0..ib { + let c_row = (ii + i) * ldc + jj; + // Register-resident C row-block: load once, accumulate the + // whole kb-loop in registers, store once. f64→f64 + // store/reload never rounds, so this is bit-identical to + // per-step write-back — same per-element op sequence, the + // C-row traffic just drops from O(kb) to O(1). + let mut acc = [F64x8::splat(0.0); TILE / LANES]; + for (v, av) in acc.iter_mut().enumerate().take(nvec) { + *av = F64x8::from_slice(&c[c_row + v * LANES..]); + } + let mut tail = [0.0f64; LANES - 1]; + for (t, tv) in tail.iter_mut().enumerate().take(ntail) { + *tv = c[c_row + nvec * LANES + t]; + } for p in 0..kb { let a_val = alpha * a[(ii + i) * lda + (kk + p)]; let va = F64x8::splat(a_val); let b_row = (kk + p) * ldb + jj; - let c_row = (ii + i) * ldc + jj; - let mut j = 0; - // Unfused `vc + va·vb` — bit-identical to the scalar tail. - while j + LANES <= jb { - let vb = F64x8::from_slice(&b[b_row + j..]); - let vc = F64x8::from_slice(&c[c_row + j..]); - (vc + va * vb).copy_to_slice(&mut c[c_row + j..]); - j += LANES; + // Reference tier: unfused `acc + va·vb`, bit-identical + // to the scalar tail. FMA tier: fused, one rounding per + // step. Monomorphized — no runtime branch. + for (v, av) in acc.iter_mut().enumerate().take(nvec) { + let vb = F64x8::from_slice(&b[b_row + v * LANES..]); + *av = if FMA { va.mul_add(vb, *av) } else { *av + va * vb }; } - while j < jb { - c[c_row + j] += a_val * b[b_row + j]; - j += 1; + for (t, tv) in tail.iter_mut().enumerate().take(ntail) { + let bv = b[b_row + nvec * LANES + t]; + *tv = if FMA { a_val.mul_add(bv, *tv) } else { *tv + a_val * bv }; } } + for (v, av) in acc.iter().enumerate().take(nvec) { + av.copy_to_slice(&mut c[c_row + v * LANES..]); + } + for (t, tv) in tail.iter().enumerate().take(ntail) { + c[c_row + nvec * LANES + t] = *tv; + } } jj += jb; } @@ -1239,4 +1312,70 @@ mod gemm_f64_tiled_tests { let c_naive = ref_naive_ijk(m, n, k, &a, k, &b, n, n); assert_bits_eq(&c, &c_naive, "denormals"); } + + // ── FMA tier ── + + /// Integer-valued corpus: values in [-9, 9], so every a·b product + /// (≤ 81) and running sum (≤ 70·81·2) is exactly representable — + /// fused and unfused steps round identically → tiers bit-identical. + fn fill_int(len: usize, salt: u64) -> Vec { + (0..len) + .map(|i| ((mix(i as u64 ^ salt) % 19) as i64 - 9) as f64) + .collect() + } + + #[test] + fn fma_bit_identical_to_reference_on_integer_operands() { + for &(m, n, k) in SHAPES { + for &(alpha, beta) in &[(1.0, 0.0), (2.0, 1.0)] { + let a = fill_int(m * k, 0xF1); + let b = fill_int(k * n, 0xF2); + let c0 = fill_int(m * n, 0xF3); + let mut c_ref = c0.clone(); + let mut c_fma = c0.clone(); + super::gemm_f64_tiled(m, n, k, alpha, &a, k, &b, n, beta, &mut c_ref, n); + super::gemm_f64_tiled_fma(m, n, k, alpha, &a, k, &b, n, beta, &mut c_fma, n); + assert_bits_eq(&c_fma, &c_ref, &format!("fma int ({m},{n},{k}) α={alpha} β={beta}")); + } + } + } + + #[test] + fn fma_matches_reference_within_fused_rounding_on_floats() { + for &(m, n, k) in SHAPES { + let a = fill(m * k, 0xF4); + let b = fill(k * n, 0xF5); + let mut c_ref = vec![0.0f64; m * n]; + let mut c_fma = vec![0.0f64; m * n]; + super::gemm_f64_tiled(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c_ref, n); + super::gemm_f64_tiled_fma(m, n, k, 1.0, &a, k, &b, n, 0.0, &mut c_fma, n); + // FMA-vs-unfused divergence scales with the SUMMAND magnitudes + // (each step drops one intermediate rounding relative to the + // running sum), not with the possibly-cancelled final value: + // ≤ k steps × ~2ε per step × O(1) summands for inputs in [-1,1). + let tol = (k as f64) * 4.0 * f64::EPSILON; + for (i, (f, r)) in c_fma.iter().zip(c_ref.iter()).enumerate() { + let tol = tol.max(1e-13 * r.abs()); + assert!( + (f - r).abs() <= tol, + "fma float ({m},{n},{k}) idx {i}: fma={f:e} ref={r:e} (Δ={:e} > tol={tol:e})", + (f - r).abs() + ); + } + } + } + + #[test] + fn fma_edge_cases_match_reference() { + // k=0 beta pass, m/n=0 no-op — same paths as the reference tier. + let c0 = fill(24, 0xF6); + let mut c_ref = c0.clone(); + let mut c_fma = c0.clone(); + super::gemm_f64_tiled(4, 6, 0, 1.0, &[], 0, &[], 6, 2.5, &mut c_ref, 6); + super::gemm_f64_tiled_fma(4, 6, 0, 1.0, &[], 0, &[], 6, 2.5, &mut c_fma, 6); + assert_bits_eq(&c_fma, &c_ref, "fma k=0"); + let mut c = c0.clone(); + super::gemm_f64_tiled_fma(0, 6, 4, 1.0, &[], 4, &[1.0; 24], 6, 0.0, &mut c, 6); + assert_bits_eq(&c, &c0, "fma m=0 no-op"); + } } From 8d85637e5f6b0e35e912b89aa70ee4a208f41d2e Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 6 Jul 2026 21:55:59 +0000 Subject: [PATCH 7/7] =?UTF-8?q?simd:=20verify-pass=20fixes=20=E2=80=94=20b?= =?UTF-8?q?ackend=20engine=20on=20the=20unfused=20reference=20tier?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three-angle adversarial verify (numerics/swap-trace/docs) on the completion diff. Substantive finding acted on: the fused tier's scalar polyfill (AVX2 arm F64x8::mul_add = per-lane f64::mul_add) lowers to a libm fma() call on baseline x86-64 builds without the fma target feature — downstream consumers do not inherit this repo's .cargo pin, so backend::native::gemm_f64 now routes through the UNFUSED reference tier instead: no libm dependence on any backend, ~7% slower than fused on pinned builds (10.0 vs 10.7 GFLOP/s), and the backend engine is now bit-identical to the certification reference. gemm_f64_tiled_fma stays public for FMA-pinned consumers. Doc fixes from the same pass: removed two stale 'gemm_f64 delegates to matrixmultiply' claims that contradicted the swap (simd.rs re-export comment, gemm_f64_tiled rustdoc); added # Panics to gemm_f64 (checked preconditions vs the old wrapper's silent UB on short slices, matching CBLAS xerbla); scoped fma determinism to per-(build,runtime) — wasm relaxed-simd fusion is implementation-defined; corrected the AVX2 vfmadd naming (per-lane polyfill, fused semantics); fixed the integer-corpus bound comment (k_max=128). Gates: clippy -D warnings clean, 2185/2185 lib tests, 4 gemm doctests. --- .claude/blackboard.md | 33 ++++++++++++++++++++++++++++----- src/backend/native.rs | 27 ++++++++++++++++++++------- src/simd.rs | 14 ++++++++------ src/simd_ops.rs | 33 +++++++++++++++++++++------------ 4 files changed, 77 insertions(+), 30 deletions(-) diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 5a3b649b..d311673a 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -24,11 +24,22 @@ below's foundation: f64→f64 store/reload never rounds ⇒ per-element op sequence unchanged ⇒ every bit-equality test green untouched. [MEASURED] ref 4.6→10.0 GF, fma 4.7→10.7 GF (2.2×). -3. **Engine swap:** `backend::native::gemm_f64` now routes to - `gemm_f64_tiled_fma` — the f64 GEMM behind `BlasFloat::backend_gemm` - / `hpc::blas_level3::blas_gemm` / batched linalg is entirely own - Rust; matrixmultiply remains only in gemm_f32 and upstream - `Array::dot` (impl_linalg.rs, untouched at ~33 GF). +3. **Engine swap:** `backend::native::gemm_f64` now routes to the + crate-native tiled kernel — the f64 GEMM behind + `BlasFloat::backend_gemm` / `hpc::blas_level3::blas_gemm` / batched + linalg is entirely own Rust; matrixmultiply remains only in gemm_f32 + and upstream `Array::dot` (impl_linalg.rs, untouched at ~33 GF). + **[REVISED post-verify] Engine = the UNFUSED reference tier**, not + fma: the verify pass surfaced a cliff — the AVX2-polyfill/scalar + `mul_add` lowers to a libm `fma()` call on baseline x86-64 builds + (consumers do NOT inherit this repo's `.cargo` target-cpu pin; CI + lands exactly there). The unfused tier has no libm dependence on any + backend, costs only ~7% vs fused on pinned builds (10.0 vs 10.7 GF), + and makes the backend engine bit-identical to the certification + reference. `gemm_f64_tiled_fma` stays public for FMA-pinned + consumers. New panic contract documented on `gemm_f64` (# Panics — + checked preconditions vs the old wrapper's silent-UB on short + slices; matches CBLAS xerbla behavior). [MEASURED, 3-engine, this VM (v3 compile → AVX2 arm, PREFERRED_F64_LANES=4; host runtime has avx512f but committed .cargo config is v3)]: @@ -46,6 +57,18 @@ still dead in native.rs `mod scalar` (f32 sibling completion); avx512f compile arm untested on CI (v3 config) — the F64x8=__m512d arm runs only on local v4 builds. +**[VERIFY OUTCOME]** 3-angle adversarial pass on the completion diff: +numerics PASS / swap-trace PASS / docs FAIL→fixed. Substantive P1 acted +on: baseline-x86 libm-fma cliff → engine revised to the unfused +reference tier (see #3 REVISED above). Doc P1s fixed: two stale +"gemm_f64 delegates to matrixmultiply" claims (simd.rs comment + +gemm_f64_tiled rustdoc) contradicted the swap in the same diff. P2s +fixed: fma determinism scoped to per-(build,runtime) (wasm relaxed-simd +fusion is implementation-defined); AVX2 vfmadd naming corrected (per- +lane f64::mul_add polyfill, fused semantics); integer-corpus bound +comment corrected (k_max=128, ≈2.1e4). All gates re-run green after +fixes: clippy -D warnings, 2185/2185 lib tests, 4 gemm doctests. + ## 2026-07-06 (later) — `ndarray::simd::gemm_f64_tiled` surfaced (operator directive) The crate-native tiled f64 GEMM graduated from dead code diff --git a/src/backend/native.rs b/src/backend/native.rs index df3c1d91..42559c2d 100644 --- a/src/backend/native.rs +++ b/src/backend/native.rs @@ -238,17 +238,30 @@ pub fn gemm_f32( /// GEMM: C = alpha * A * B + beta * C (f64, row-major). /// -/// Engine: the crate-native `simd_ops::gemm_f64_tiled_fma` (F64x8-vectorized -/// fused tier — entirely in-crate Rust, no external matmul engine). The -/// unfused sibling `ndarray::simd::gemm_f64_tiled` is the bit-exact -/// certification reference; this entry favors throughput. Swapped from -/// `matrixmultiply::dgemm` — results agree to BLAS tolerance (last-ulp -/// accumulation-order/fusion differences), not bit-for-bit. +/// Engine: the crate-native `simd_ops::gemm_f64_tiled` (F64x8-vectorized, +/// unfused reference tier — entirely in-crate Rust, no external matmul +/// engine for f64; `gemm_f32` below still delegates to `matrixmultiply`). +/// The unfused tier is chosen deliberately: it is bit-identical to the +/// certification reference on every backend AND has no dependence on the +/// `fma` target feature (the fused tier's scalar polyfill can lower to a +/// libm `fma()` call on baseline x86-64 builds that don't inherit this +/// repo's `.cargo` target-cpu pin). Consumers that pin an FMA-capable +/// target and want the fused tier call `ndarray::simd::gemm_f64_tiled_fma` +/// directly. Swapped from `matrixmultiply::dgemm` — results agree to BLAS +/// tolerance (last-ulp accumulation-order differences), not bit-for-bit. +/// +/// # Panics +/// +/// Unlike the previous `matrixmultiply` wrapper (which accepted `lda < k` +/// silently and was UB on short slices), this entry inherits the tiled +/// kernel's checked preconditions: panics if `lda < k`, `ldb < n`, +/// `ldc < n`, or a slice is shorter than `(rows − 1)·ld + cols` — matching +/// the CBLAS backends, which reject `lda < max(1, k)` via `xerbla`. pub fn gemm_f64( m: usize, n: usize, k: usize, alpha: f64, a: &[f64], lda: usize, b: &[f64], ldb: usize, beta: f64, c: &mut [f64], ldc: usize, ) { - crate::simd_ops::gemm_f64_tiled_fma(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); + crate::simd_ops::gemm_f64_tiled(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } // ─── GEMV dispatch ─────────────────────────────────────────────── diff --git a/src/simd.rs b/src/simd.rs index 66ce3fa4..4c7b95f1 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -570,13 +570,15 @@ pub use crate::simd_ops::{array_chunks, array_chunks_checked, array_windows, arr // contract: unfused mul+add in ascending-k order per element → bit-identical // on every backend (AVX-512/AVX2/NEON/WASM/scalar) and, at α=1 β=0, // bit-identical to the naive triple-loop reference. This is the in-crate -// ground-truth GEMM for probes/certification — `backend::native::gemm_f64` -// delegates to the external `matrixmultiply` crate instead. The kernel is -// alloc-free, but `pub mod simd`/`simd_ops` are std-gated in lib.rs, so -// like every kernel here it is reachable only in `std` builds today. +// ground-truth GEMM for probes/certification AND the engine behind +// `backend::native::gemm_f64` (own Rust in the f64 BLAS path; the f32 +// sibling still delegates to the external `matrixmultiply` crate). // `gemm_f64_tiled_fma` is the fast fused tier (same tiling/order, one -// rounding per step) — the engine behind `backend::native::gemm_f64`; -// the unfused reference tier stays the certification ground truth. +// rounding per step) for consumers on FMA-pinned targets — not the +// backend engine, because its scalar polyfill can lower to libm `fma()` +// on baseline builds. Both kernels are alloc-free, but `pub mod +// simd`/`simd_ops` are std-gated in lib.rs, so they are reachable only +// in `std` builds today. pub use crate::simd_ops::{gemm_f64_tiled, gemm_f64_tiled_fma}; pub use crate::simd_soa::MultiLaneColumn; diff --git a/src/simd_ops.rs b/src/simd_ops.rs index fcf07e65..2c0e213b 100644 --- a/src/simd_ops.rs +++ b/src/simd_ops.rs @@ -890,10 +890,11 @@ mod add_mul_tests { /// The crate-native f64 GEMM: a TILE=64-blocked loop nest with the /// innermost j-run vectorized on the dispatched [`F64x8`] (AVX-512 / /// AVX2 / NEON / WASM-SIMD128 / scalar — one source, every backend). -/// Unlike `backend::native::gemm_f64`, which delegates to the external -/// `matrixmultiply` crate, this kernel is entirely in-crate and carries -/// a bit-exactness contract, so probes can certify other kernels -/// against it as ground truth. +/// Entirely in-crate with a bit-exactness contract, so probes can +/// certify other kernels against it as ground truth — and since it IS +/// the engine behind `backend::native::gemm_f64`, the native f64 BLAS +/// path carries the same auditable numerics (the f32 sibling still +/// delegates to the external `matrixmultiply` crate). /// /// # Bit-exactness contract /// @@ -953,14 +954,21 @@ pub fn gemm_f64_tiled( /// Fast tier of [`gemm_f64_tiled`]: same tiling, same ascending-`p` /// per-element order, but each step is a **fused** multiply-add /// (`c = fma(alpha·A[i,p], B[p,j], c)` — one rounding per step via -/// `F64x8::mul_add`: `vfmadd`/`vfmaq_f64` on x86/NEON, `f64::mul_add` -/// on the scalar tail and polyfill backends). +/// `F64x8::mul_add`; fused semantics on every backend, though the +/// lowering varies: native `vfmadd` on the AVX-512 arm, `vfmaq_f64` on +/// NEON, per-lane `f64::mul_add` on the AVX2 polyfill and the scalar +/// tail — which can become a libm `fma()` call on targets built +/// without the `fma` feature: correct, but slow. Prefer the reference +/// tier, or an FMA-capable target pin, in that situation). /// /// Semantics relative to the reference tier: -/// - **Deterministic per build**, but NOT bit-stable across backends: -/// fusion quality differs (WASM without `relaxed-simd` has no fused -/// vector `mul_add`, so its vector lanes match the unfused reference -/// while the scalar tail is still fused). +/// - **Deterministic per (build, runtime)**, but NOT bit-stable across +/// backends: fusion differs (WASM without `relaxed-simd` has no +/// fused vector `mul_add`, so its vector lanes match the unfused +/// reference while the scalar tail is still fused; WASM WITH +/// `relaxed-simd` uses `f64x2_relaxed_madd`, whose fusion is +/// implementation-defined — the same wasm binary may round +/// differently on different runtimes). /// - **Bit-identical to [`gemm_f64_tiled`]** whenever every /// `a_val·b` product and running sum is exactly representable — /// e.g. integer-valued data with products and accumulation inside @@ -1316,8 +1324,9 @@ mod gemm_f64_tiled_tests { // ── FMA tier ── /// Integer-valued corpus: values in [-9, 9], so every a·b product - /// (≤ 81) and running sum (≤ 70·81·2) is exactly representable — - /// fused and unfused steps round identically → tiers bit-identical. + /// (≤ 81) and running sum (≤ k_max·81·|α| + 9 = 128·81·2 + 9 ≈ 2.1e4, + /// far below 2^53) is exactly representable — fused and unfused + /// steps round identically → tiers bit-identical. fn fill_int(len: usize, salt: u64) -> Vec { (0..len) .map(|i| ((mix(i as u64 ^ salt) % 19) as i64 - 9) as f64)