From 6cf6b170068b72553d8580875b7c664788308d83 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 4 Sep 2026 14:36:58 +0000 Subject: [PATCH 1/4] Correct the GEMM/array_windows attribution in the README PR #298 claimed the 512x512/1024x1024 GEMM numbers depend on simd_ops::array_windows/array_chunks. Tracing the actual call graph (backend::native::gemm_f32 -> matrixmultiply crate's Goto algorithm; simd_ops::gemm_f64_tiled -> fixed TILE=64 blocking with F64x8 register accumulation) shows neither path calls array_windows or array_chunks at all -- their only production call site is hpc::blake3's 64-byte block chunking. Corrected the README to describe the GEMM kernels' real mechanism and to attribute array_windows/array_chunks to where they're actually used, without losing the real (separately verified) array_windows-vs-JIT latency finding. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_012wrzeZAdwGYTCKoxamwQht --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 599246a5..c0d31ff4 100644 --- a/README.md +++ b/README.md @@ -111,9 +111,9 @@ Detection happens once on first access via `LazyLock` — a single CPU | 1024 x 1024 | ~13 GFLOPS | 139 GFLOPS | ~120 GFLOPS | ~3,500 GFLOPS | | 2048 x 2048 | ~13 GFLOPS | ~150 GFLOPS | ~140 GFLOPS | ~5,000 GFLOPS | -Upstream hits a cache cliff at 1024 x 1024: no tiling, no threading, no microkernel. The fork uses the Goto algorithm with cache blocking (L1/L2/L3) and achieves 10.5x throughput — on par with NumPy's decades-old OpenBLAS. +Upstream hits a cache cliff at 1024 x 1024: no tiling, no threading, no microkernel. The fork's f32 GEMM (`backend::native::gemm_f32`) delegates to the `matrixmultiply` crate's Goto algorithm with cache blocking (L1/L2/L3); f64 GEMM (`ndarray::simd::gemm_f64_tiled`) is an in-crate fixed-tile (`TILE=64`) kernel with `F64x8`-vectorized register accumulation. Either way the fork achieves 10.5x throughput over upstream — on par with NumPy's decades-old OpenBLAS. -The 512 x 512 and 1024 x 1024 numbers specifically depend on `simd_ops::array_windows`/`array_windows_checked` (a stable-Rust, const-generic reimplementation of nightly `slice::array_windows::()`, giving overlapping `&[T; N]` references with no bounds check per step) paired with `array_chunks` (the non-overlapping counterpart) and the polyfilled `mul_add`/`add_mul_f32`/`add_mul_f64` FMA primitives. That combination is the same blocking/window discipline the original C blasgraph kernels relied on for cache-line-exact reuse at these matrix sizes, and it was benchmarked per-call against the original C blasgraph kernel and an actual Cranelift-JIT-compiled inner loop as alternatives: roughly 7 ns (original C blasgraph) vs 12 ns (Cranelift JIT) vs 17 ns (static `array_windows` microkernel) — landing close to the JIT without paying for one (compile latency, codegen complexity, or the dependency). Remove `array_windows` and the fork's GEMM falls back toward upstream's unblocked cache-cliff behavior at exactly this size range. +`simd_ops::array_windows`/`array_windows_checked` (a stable-Rust, const-generic reimplementation of nightly `slice::array_windows::()`, giving overlapping `&[T; N]` references with no bounds check per step) and `array_chunks` (the non-overlapping counterpart) are separate fork primitives, not part of the GEMM path above — their production call site today is block-chunked hashing (`hpc::blake3`'s 64-byte `BLOCK_LEN` iteration via `array_chunks`). They're available to any consumer needing overlapping- or chunked-window traversal with the same cache-line-exact discipline the original C blasgraph kernels used, and were benchmarked standalone against a Cranelift-JIT-compiled inner loop as an alternative, landing close to the JIT's per-call latency without paying for one (compile latency, codegen complexity, or the dependency). ### Data Types Beyond f32/f64 From ce3563103cb929dd81098352d54c401aa62577b2 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 4 Sep 2026 14:46:05 +0000 Subject: [PATCH 2/4] Use array_chunks in the f32->f64 cosine SIMD loop cosine_f32_to_f64_simd indexed a[off + j] / b[off + j] per element, so each of the 16 loads per 8-element chunk carried a bounds check against a dynamic length. array_chunks:: yields &[f32; 8], where the inner index is provably in bounds, and from_array consumes the widened buffer without a second length check. The op order is unchanged, so results are bit-identical; verified standalone across every length 0..300 including unequal-length inputs, since the crate's lib tests currently do not compile (pre-existing missing U16x8 imports in simd_avx512.rs test modules, unrelated). This is the first of the mechanism-level consolidations: hand-rolled chunk index bookkeeping replaced by the const-generic window primitive that simd_ops.rs documents as the foundation for exactly these kernels. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_012wrzeZAdwGYTCKoxamwQht --- src/heel_f64x8.rs | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/src/heel_f64x8.rs b/src/heel_f64x8.rs index 87ff42bb..76572567 100644 --- a/src/heel_f64x8.rs +++ b/src/heel_f64x8.rs @@ -9,7 +9,7 @@ //! Scalar: [f64; 8] fallback //! Consumer writes `crate::simd::F64x8`. The polyfill handles the rest. -use crate::simd::F64x8; +use crate::simd::{array_chunks, F64x8}; /// Compute weighted dot product of 8 HEEL plane distances. /// @@ -158,14 +158,15 @@ pub fn cosine_f32_to_f64_simd(a: &[f32], b: &[f32]) -> f64 { let mut buf_a = [0.0f64; 8]; let mut buf_b = [0.0f64; 8]; - for i in 0..chunks { - let off = i * 8; + for (ca, cb) in + array_chunks::(&a[..n]).zip(array_chunks::(&b[..n])) + { for j in 0..8 { - buf_a[j] = a[off + j] as f64; - buf_b[j] = b[off + j] as f64; + buf_a[j] = ca[j] as f64; + buf_b[j] = cb[j] as f64; } - let va = F64x8::from_slice(&buf_a); - let vb = F64x8::from_slice(&buf_b); + let va = F64x8::from_array(buf_a); + let vb = F64x8::from_array(buf_b); dot_acc = va.mul_add(vb, dot_acc); na_acc = va.mul_add(va, na_acc); nb_acc = vb.mul_add(vb, nb_acc); From 6580310208e8997459af40fe80f495c565e746cf Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 4 Sep 2026 14:47:59 +0000 Subject: [PATCH 3/4] Fix rustfmt, and stop attributing the GEMM table to unverified paths The array_chunks loop was hand-wrapped and failed format/stable. The GEMM Performance table has now been mis-attributed twice, so this stops guessing and marks its provenance unverified instead. Array::dot() calls matrixmultiply::sgemm/dgemm, which for f32 is the same engine backend::native::gemm_f32 uses, so that path has no fork-vs-upstream engine difference to credit; matrixmultiply also has Goto blocking and microkernels, making the "upstream has no tiling/microkernel" line false. gemm_f64_tiled is the one fork-local kernel and Array::dot() does not reach it. Also narrows the window-primitive claim: blake3 and the cosine kernel call array_chunks; array_windows and the checked variants have no in-crate production caller. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_012wrzeZAdwGYTCKoxamwQht --- README.md | 13 ++++++++++--- src/heel_f64x8.rs | 4 +--- 2 files changed, 11 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index c0d31ff4..eb9a6d17 100644 --- a/README.md +++ b/README.md @@ -111,9 +111,16 @@ Detection happens once on first access via `LazyLock` — a single CPU | 1024 x 1024 | ~13 GFLOPS | 139 GFLOPS | ~120 GFLOPS | ~3,500 GFLOPS | | 2048 x 2048 | ~13 GFLOPS | ~150 GFLOPS | ~140 GFLOPS | ~5,000 GFLOPS | -Upstream hits a cache cliff at 1024 x 1024: no tiling, no threading, no microkernel. The fork's f32 GEMM (`backend::native::gemm_f32`) delegates to the `matrixmultiply` crate's Goto algorithm with cache blocking (L1/L2/L3); f64 GEMM (`ndarray::simd::gemm_f64_tiled`) is an in-crate fixed-tile (`TILE=64`) kernel with `F64x8`-vectorized register accumulation. Either way the fork achieves 10.5x throughput over upstream — on par with NumPy's decades-old OpenBLAS. - -`simd_ops::array_windows`/`array_windows_checked` (a stable-Rust, const-generic reimplementation of nightly `slice::array_windows::()`, giving overlapping `&[T; N]` references with no bounds check per step) and `array_chunks` (the non-overlapping counterpart) are separate fork primitives, not part of the GEMM path above — their production call site today is block-chunked hashing (`hpc::blake3`'s 64-byte `BLOCK_LEN` iteration via `array_chunks`). They're available to any consumer needing overlapping- or chunked-window traversal with the same cache-line-exact discipline the original C blasgraph kernels used, and were benchmarked standalone against a Cranelift-JIT-compiled inner loop as an alternative, landing close to the JIT's per-call latency without paying for one (compile latency, codegen complexity, or the dependency). +> **Provenance of this table is unverified.** The numbers predate the current +> tree and the benchmark that produced them is not in the repository, so the +> API and element type they measured cannot be identified. Do not cite them as +> a fork-vs-upstream result until they are reproduced. What the code does say: +> +> - `Array::dot()` calls `matrixmultiply::sgemm`/`dgemm` (`src/linalg/impl_linalg.rs:503,522`) — for f32 this is the **same engine** `backend::native::gemm_f32` uses (`src/backend/native.rs:220`), so on that path there is no fork-vs-upstream engine difference to attribute a speedup to. +> - `matrixmultiply` implements Goto-style cache blocking with microkernels, so "upstream has no tiling/microkernel" is false. +> - The one genuinely fork-local GEMM kernel is `simd::gemm_f64_tiled` (`simd_ops.rs:947`) — fixed `TILE=64`, `F64x8` register accumulation, reached via `backend::native::gemm_f64` / `BlasLevel3::blas_gemm`, **not** via `Array::dot()`. + +`simd_ops::array_chunks` walks a slice as non-overlapping `&[T; N]` windows; `array_windows` is the overlapping counterpart (a stable-Rust equivalent of nightly `slice::array_windows::()`). Both pin the window size at the call site so it feeds `F32x16::from_array` / `F64x8::from_array` directly, and both drop the per-element bounds check a dynamically-indexed loop pays. Current in-crate call sites: `hpc::blake3` (64-byte block chunking) and `heel_f64x8::cosine_f32_to_f64_simd`, both via `array_chunks`; `array_windows`, `array_windows_checked`, and `array_chunks_checked` are exported but have no in-crate production caller yet. They are the traversal primitive the hand-rolled BLAS-graph/bgz17 kernels are built on, where the const-generic window landed close to a Cranelift-JIT'd inner loop without paying for a JIT — see `src/simd_ops.rs` module docs. ### Data Types Beyond f32/f64 diff --git a/src/heel_f64x8.rs b/src/heel_f64x8.rs index 76572567..9fd79540 100644 --- a/src/heel_f64x8.rs +++ b/src/heel_f64x8.rs @@ -158,9 +158,7 @@ pub fn cosine_f32_to_f64_simd(a: &[f32], b: &[f32]) -> f64 { let mut buf_a = [0.0f64; 8]; let mut buf_b = [0.0f64; 8]; - for (ca, cb) in - array_chunks::(&a[..n]).zip(array_chunks::(&b[..n])) - { + for (ca, cb) in array_chunks::(&a[..n]).zip(array_chunks::(&b[..n])) { for j in 0..8 { buf_a[j] = ca[j] as f64; buf_b[j] = cb[j] as f64; From 506f7d74e7b5866a0d50e419063b17e31422f2ef Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 4 Sep 2026 14:56:40 +0000 Subject: [PATCH 4/4] Restore the AVX-512 test modules to a compiling state int_simd_tests referenced I8x16, U64x8 and U16x8 without importing them, so the lib test target failed to compile whenever AVX-512 was enabled at build time. CI does not pass -C target-cpu, so those modules were cfg'd out there and CI stayed green while the tests were broken in the configuration .cargo/config.toml actually ships (x86-64-v4). All three types are already exported from crate::simd; adding them to the existing use list is the whole fix. The suite now compiles and runs under target-cpu=native: 2322 passed, 0 failed, 30 ignored. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_012wrzeZAdwGYTCKoxamwQht --- src/simd_avx512.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simd_avx512.rs b/src/simd_avx512.rs index 89dc104b..71f92cc9 100644 --- a/src/simd_avx512.rs +++ b/src/simd_avx512.rs @@ -4654,7 +4654,7 @@ mod tier3_tests { #[cfg(all(test, target_feature = "avx512f"))] mod int_simd_tests { - use crate::simd::{I16x16, I16x32, I8x32, I8x64}; + use crate::simd::{I16x16, I16x32, I8x16, I8x32, I8x64, U16x8, U64x8}; #[test] fn i8x64_add_pair_to_constant() {