diff --git a/README.md b/README.md index 599246a5..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 uses the Goto algorithm with cache blocking (L1/L2/L3) and achieves 10.5x throughput — 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. +> **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 87ff42bb..9fd79540 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,13 @@ 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); 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() {