Merge pull request #299 from AdaWorldAPI/claude/ndarray-simd-tract-o3jfrn

Correct the GEMM/array_windows attribution, and start consolidating onto array_chunks

Author: AdaWorldAPI

README.md

@@ -111,9 +111,16 @@ Detection happens once on first access via `LazyLock<SimdCaps>` — 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::<N>()`, 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::<N>()`). 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
 

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::<f32, 8>(&a[..n]).zip(array_chunks::<f32, 8>(&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);

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() {