diff --git a/.claude/blackboard.md b/.claude/blackboard.md index 61d0e1fc..d311673a 100644 --- a/.claude/blackboard.md +++ b/.claude/blackboard.md @@ -3,6 +3,208 @@ > **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 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)]: +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. + +**[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 +(`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, 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- +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` +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). + +**[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` 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..42559c2d 100644 --- a/src/backend/native.rs +++ b/src/backend/native.rs @@ -237,32 +237,31 @@ pub fn gemm_f32( } /// GEMM: C = alpha * A * B + beta * C (f64, row-major). +/// +/// 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, ) { - 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(m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } // ─── GEMV dispatch ─────────────────────────────────────────────── @@ -468,52 +467,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..4c7b95f1 100644 --- a/src/simd.rs +++ b/src/simd.rs @@ -566,6 +566,20 @@ 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 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) 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; pub use crate::hpc::quantized::{ diff --git a/src/simd_ops.rs b/src/simd_ops.rs index 455f915e..2c0e213b 100644 --- a/src/simd_ops.rs +++ b/src/simd_ops.rs @@ -880,3 +880,511 @@ 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). +/// 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 +/// +/// 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 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. 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 (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 +/// +/// ``` +/// 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, +) { + 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`; 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, 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 +/// `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; + } + 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() >= 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() >= extent(m, 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); + 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; + // 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 }; + } + 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; + } + 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() { + // 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); + 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"); + } + + // ── FMA tier ── + + /// Integer-valued corpus: values in [-9, 9], so every a·b product + /// (≤ 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) + .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"); + } +}