tannevaled · Oct 4, 2026 · Oct 4, 2026
diff --git a/docs/performance.md b/docs/performance.md
@@ -4,7 +4,7 @@ Honest, reproducible head-to-head of `go-ndarray/ndarray` against **NumPy
 2.2.4** on identical hardware. "On n'a pas le droit de se tromper": every number
 here is measured, every win is real, and where NumPy still leads it says so.
 
-> **Where it stands at v0.2.6.** This page grew round by round: the sections
+> **Where it stands at v0.4.1.** This page grew round by round: the sections
 > up to *Where go-ndarray still loses* were measured on an arm64 Tart VM
 > (4 vCPU, NumPy 2.2.4, OpenBLAS 0.3.29); the dated sections after it on an
 > AMD Zen 3, 16 cores, against **NumPy 2.5.3 + OpenBLAS 0.3.34**, except the
@@ -19,8 +19,9 @@ here is measured, every win is real, and where NumPy still leads it says so.
 >   beat NumPy at every size measured.
 > - **At parity:** `SumAxis(0)`, `Inner` against 16-thread OpenBLAS, `Log` on
 >   one core.
-> - **Slower:** `MatMul` against 16-thread OpenBLAS (0.55× at 1024², 0.3× at
->   256² and below); `Dot`/mat·vec against NumPy's threaded BLAS (0.16×); `Exp`
+> - **Slower:** `MatMul` against 16-thread OpenBLAS at 512² and up (0.42× at
+>   512², 0.51× at 1024², measured in one run with OpenBLAS; at 256² ours is
+>   1.44× faster, see *Parallel GEMM: start once, wait on work*); `Dot`/mat·vec against NumPy's threaded BLAS (0.16×); `Exp`
 >   below 256 Ki elements (0.82×); the *allocating* forms of elementwise ops up
 >   to 256 Ki elements, `Concatenate`/`Stack` and slice copies outside a
 >   `Workspace` (0.2–0.5×).
@@ -447,9 +448,13 @@ below 256²). Where the parallel loss is, measured on the Zen 3:
   each ran at 41 GFLOP/s (656 in total), while the parallel 1024² GEMM ran its
   kernel at 15–17 GFLOP/s per core. Per-core efficiency already drops from 40
   to ~30 GFLOP/s at 2 cores.
-- It is not scheduling overhead: 83% of the parallel profile is the
-  micro-kernel itself, so the kernel runs slower when cores share work. Why is
-  **not established**. The shared packed B panel is the suspect (on this part
+- ~~It is not scheduling overhead: 83% of the parallel profile is the
+  micro-kernel itself, so the kernel runs slower when cores share work.~~
+  **Wrong, corrected the same day:** a profile of cycles samples only the cores
+  that run, never the ones that wait. `perf stat` shows 7.1 of 16 CPUs busy
+  during a 1024² product: the loss was mostly waiting, and the fix is in
+  [Parallel GEMM: start once, wait on work](#parallel-gemm-start-once-wait-on-work-2026-10-04).
+  Why the kernel itself also runs slower in parallel is still **not established**. The shared packed B panel is the suspect (on this part
   the last-level cache is split between core complexes, so the panel is partly
   remote), but the drop already shows at 2 cores of one complex, which that
   explanation does not cover.
@@ -462,9 +467,8 @@ below 256²). Where the parallel loss is, measured on the Zen 3:
   rows·cols·k multiply-adds, and a strided pack costs far more per element than
   an FMA, so tiles under ~256×256 lose.
 
-The next step is to find where the kernel stalls at 2 cores (hardware
-counters: cache misses per FMA, serial against parallel) before choosing a
-lever; then an AVX-512 kernel for the hosts that have it.
+The hardware counters were read the same day; see the next-but-one
+section. An AVX-512 kernel for the hosts that have it remains open.
 
 Correctness was checked on real hardware for five of the six 64-bit targets:
 amd64 (Zen 3, FMA path), arm64 (cfarm185), ppc64le (cfarm120), riscv64 (cfarm94)
@@ -673,6 +677,86 @@ parallel 4 Mi sum is bimodal on this host (1.6 ms or 7 ms in either version),
 so only the minima are compared. No NumPy or BLAS is installable there (the
 host has no route to package mirrors), so there is no external reference.
 
+## Parallel GEMM: start once, wait on work (2026-10-04)
+
+**Diagnosis, from hardware counters** (`perf stat`, Zen 3 guest cfarm421,
+steady state, 1024², 3 s per configuration):
+
+| workers | CPUs busy | IPC | instructions |
+|---|---|---|---|
+| 1 | 1.00 | 3.35 | 37.5 G |
+| 2 | 1.91 | 2.00 | 37.6 G |
+| 8 | 5.58 | 2.16 | 172 G (192 products) |
+| 16 | 7.08 | 2.28 | 193 G (226 products) |
+
+The instruction count does not grow, so nobody busy-waits. L1 misses do not
+grow either. What changes is that at 16 workers the CPUs are busy only 44% of
+the time. The driver forked and joined goroutines twice per (jc, pc) block of
+the product, once to pack B and once for the row bands, 16 times in a 1024²
+product. A fork/join wakes parked threads one after another. A micro-benchmark
+on the same guest (16 goroutines, 700 µs of work each) measures its cost at
+**26 µs with 2 workers, 165 µs with 8, 330 µs with 16**. 16 × 330 µs is ~40%
+of a 12.5 ms product.
+
+Two hypotheses were tested and **refuted** first:
+
+- SMT siblings. The guest has 128 vCPUs, as many as the host has threads.
+  Two serial GEMMs pinned on vCPU 16 and any other vCPU slow each other by
+  0–10%, never by half.
+- A block of A larger than L2. With MC = 256, the A block is 512 KiB, all of
+  Zen 3's L2, against BLIS's MC = 72. Sweeping MC over {72, 120, 144, 192, 256}
+  changed nothing beyond the noise at 2 workers, and gave +17% at 16 workers
+  for 1024² but −8% for 512².
+
+**The change.** The workers start once per product. They synchronise on
+**work done**, never on workers arrived: a block's bands start once all of its
+B panels are packed and all bands of the previous block are done, and the
+packing is shared in claims of four NR-wide panels. A first version used a
+barrier that waited for every worker, and it put the wake-up chain straight
+back on the critical path for small products: 96³ on 8 POWER9 cores ran at
+0.57×. B is double-buffered, so block i+1 is packed while block i is still
+being multiplied. Each C tile is still accumulated in pc order, so the result
+is bit-identical to the serial product. The pack buffers moved to two pools,
+one for A and one for B: a pair per worker had held 16 MiB of unused B buffers
+at 16 workers, which the collector empties and the next call zeroes again
+(`memclr` page faults, 20% of the cycles in a KVM spinlock in one profile).
+
+**Measured.** Every figure is the median of interleaved runs against the
+previous driver.
+
+| | 2 workers | 8 workers | 16 workers |
+|---|---|---|---|
+| Zen 3 guest, 1024² | 1.06× | 1.33× | 1.39× |
+| Zen 3 guest, 512² | 1.20× | 1.20× | 1.33× |
+| Zen 3 guest, 256² | 0.92× | 1.19× | 1.24× |
+| Zen 3 guest, 128² | **0.89×** | 1.50× | 1.16× |
+| Zen 3 guest, (1024×256)·(256×1024) | **0.82×** | 1.07× | 0.89× (p25 1.00×) |
+| POWER9 bare metal, 512² | 1.34× | 1.37× | |
+| POWER9 bare metal, 256² | 1.61× | 1.78× | |
+| POWER9 bare metal, (1024×256)·(256×1024) | 1.13× | 1.18× | |
+
+The Zen 3 runs were 12 rounds and the POWER9 runs 6, except 24 for 128² and
+256² at 8 workers. Those two are bimodal on that machine, in both drivers, and
+6 rounds had shown 256² at 0.72×. POWER9 used one thread per core (CPUs 0, 4,
+8, …).
+
+**Still behind.** On the Zen 3 guest, two workers lose 8–18% on mid-size
+products, where a fork/join costs only 26 µs. Bare-metal POWER9 gains at two
+workers, so this is not a property of the scheme, and no second code path was
+added for it. Against OpenBLAS with 16 threads, measured three ways in one run:
+
+| n | OpenBLAS-16 | previous driver | new driver |
+|---|---|---|---|
+| 256 | 708 µs | 1.06× | **1.44×** |
+| 512 | 789 µs | 0.29× | 0.42× |
+| 1024 | 4.60 ms | 0.43× | 0.51× |
+
+(× = OpenBLAS time / ours, so higher is faster). On this guest, absolute times
+drift through the day: the previous driver measured 0.54–0.66× of OpenBLAS in
+the morning. Only ratios taken in the same run are comparable. Not measured on
+Apple M-series, whose machine was shared with other jobs at the time. The bands
+stay dynamic, so its slower efficiency cores are handled as before.
+
 ## SIMD coverage
 
 - **amd64 (SSE2)** ships hand-vectorized `sum` (4-accumulator `ADDPD`), `sqrt`

diff --git a/docs/security.md b/docs/security.md
@@ -1,6 +1,6 @@
 # Security
 
-> Mirrored from [`SECURITY.md`](https://github.com/go-ndarray/ndarray/blob/main/SECURITY.md) at v0.3.0.
+> Mirrored from [`SECURITY.md`](https://github.com/go-ndarray/ndarray/blob/main/SECURITY.md) at v0.4.1.
 
 ## Reporting
 
@@ -41,11 +41,40 @@ bounds checks do not cover, so they are held to the following:
   `TestEmptyArraysWithHugeAxesReturnPromptly` covers concatenation, the axis
   reductions and the scans.
 
+- **Memory is the caller's to bound.** Validation guarantees a shape is
+  representable, not that it fits in memory, and a Go program that runs out of
+  memory dies (NumPy raises `MemoryError`). Broadcasting, `Outer`, `MatMul` and
+  summing the empty axis of `(0, n)` produce more than their operands hold;
+  with untrusted shapes, check the result's size first.
+- **Fuzzed against an oracle.** `FuzzOps` ([`fuzz_test.go`](https://github.com/go-ndarray/ndarray/blob/main/fuzz_test.go)) drives slices,
+  reshapes, transposes, squeezes, broadcasting arithmetic, axis reductions,
+  `Take`, `MaskSelect`, concatenation, `MatMul` and scans from arbitrary bytes,
+  and compares every result element by element with an oracle written in the
+  test from NumPy's definitions (Python's `slice.indices`, the broadcasting
+  rule, ...). The library must reject exactly the calls the oracle rejects and
+  never panic.
+
 The minimum Go version (`go` in `go.mod`) is kept on a release that fixes the
 standard library's known vulnerabilities; `govulncheck ./...` reports none.
 
 ## Audit log
 
+**2026-10-04 (v0.4.0).**
+
+- `FuzzOps` ran 47 million executions without a divergence. A first version
+  that only checked an array against its own `Copy` was **blind**: with
+  `Squeeze` deliberately given wrong strides it found nothing in 31 million
+  executions, because a wrong view reads the same wrong values both ways. The
+  oracle version finds the same sabotage in 11 seconds; its failing input is
+  kept in `testdata/fuzz/FuzzOps`.
+- It found two places where the library and NumPy differ. `Slice` with fewer
+  indices than axes was an error; it now takes the rest whole, as NumPy does.
+  `MaskSelect` accepts a mask that broadcasts to the array, which NumPy
+  rejects; this is kept as a documented extension.
+- `SumAxis`/`ProdAxis` along a zero-length axis now return 0s/1s as NumPy does
+  (they were errors, unlike `Sum`/`Prod` of an empty array). This makes the
+  memory caveat above concrete, and it is documented there.
+
 **2026-10-04 (v0.3.0).**
 
 - `govulncheck` found no reachable vulnerability. At `go 1.26.4`, 10 standard