diff --git a/docs/performance.md b/docs/performance.md index c6786ac..ef97ae4 100644 --- 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 index 13313c1..ab4903c 100644 --- 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