test: parallel stress of the helper pool; docs: Zen 3 re-measurement

README.md

@@ -88,7 +88,7 @@ is wrong on subnormals
   one core.
 - **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); `Dot`/mat·vec against NumPy's threaded BLAS (0.16×, measured before v0.6.0's helper pool, which took the 8-core dot from 625 to 58 µs on POWER9; not yet re-measured on the Zen 3); `Exp`
+  1.44× faster); `Dot` and mat·vec against NumPy's 16-thread BLAS (0.69× and 0.52× since v0.6.0's helper pool, measured in one run; 0.16× before); `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×: Go's allocator and garbage collector, which

SECURITY.md

@@ -57,6 +57,18 @@ standard library's known vulnerabilities; `govulncheck ./...` reports none.
 
 ## Audit log
 
+**2026-10-07 (after v0.6.1, tests only).** `TestParallelStress` runs the
+parallel kernels from 6 goroutines at once while another keeps changing
+`GOMAXPROCS`, with the thresholds lowered so that tiny arrays go parallel. The
+pool is contended, and some calls take the fork fallback. Every result is
+compared with the serial kernel on integer data, so the comparison is exact.
+A lost block hangs the call rather than corrupting it, so the test also fails,
+by name, when a call does not return in time. Both kinds of sabotage were
+caught: a pool that skips the first block of each run (hang) and a fork
+fallback that drops the last element (wrong value). It passes under `-race`
+and on real POWER9 (59,154 calls in 5 runs), X-Gene arm64 (65,798) and
+Loongson 3C5000L (36,174).
+
 **2026-10-07 (v0.6.1).** An audit of the helper pool added in v0.6.0 found
 two defects. Both are fixed, each with a test that fails without the fix.
 

docs/perf.md

@@ -590,8 +590,29 @@ operations speeding up (Sqrt of 1 Ki, 1.5–2.2×). Run alone, those are equal.
 That gain belonged to whatever state the earlier benchmarks had left, not to
 this change, so only fresh-process figures are reported here. The cost: after
 an operation, the helpers spend up to 200 µs of CPU polling. They are never
-stopped (see the README, *Goroutines*). Not yet measured on amd64 or on Apple
-M, whose hosts were unreachable or loaded at the time.
+stopped (see the README, *Goroutines*).
+
+**Zen 3 (cfarm421, back up on 2026-10-07)**, each case in a fresh process,
+6 interleaved rounds, v0.6.1 against v0.5.2 (median):
+
+| | 2 workers | 8 workers | 16 workers | 16: before → after |
+|---|---|---|---|---|
+| Dot 2^20 | 1.24× | 2.06× | 2.08× | 120 → 57 µs |
+| MatVec 1024² | 1.41× | 2.26× | 1.92× | 94 → 49 µs |
+| Sum 256 Ki | 1.61× | 2.12× | 2.31× | 37 → 16 µs |
+| Sum 4 Mi | 1.10× | 1.95× | 2.07× | 184 → 89 µs |
+| Max 256 Ki | 1.47× | 2.38× | 2.36× | 41 → 18 µs |
+| SumAxis(1) | 1.25× | 2.20× | 1.63× | 93 → 57 µs |
+| AddInto 256 Ki | 1.39× | 1.22× | 1.29× | 74 → 57 µs |
+| Exp 256 Ki | 1.05× | **0.85×** | 1.13× | |
+| Exp 4 Mi | 1.28× | 1.09× | **0.94×** | |
+| Chain 4 Mi | 1.74× | 1.04× | **0.88×** | |
+
+In the same run, NumPy 2.5.3 with 16 OpenBLAS threads took 39 µs for the 2^20
+dot and 25 µs for the 1024² mat·vec, so go-ndarray is at **0.69×** and **0.52×**
+of it (0.16× before the pool). The losses at the top of each table are long,
+memory-bound chains, where all workers stream at once; not investigated further.
+Apple M not measured (shared machine).
 
 ## Log and Log10 (2026-10-04)
 

internal/kernels/stress_test.go

@@ -0,0 +1,128 @@
+package kernels
+
+import (
+	"math"
+	"math/rand"
+	"runtime"
+	"sync"
+	"sync/atomic"
+	"testing"
+	"time"
+)
+
+// TestParallelStress is a differential stress test of the parallel paths and
+// the helper pool behind them. Several goroutines run parallel kernels at once
+// (so the pool is contended and some calls take the fork fallback) while
+// another keeps changing GOMAXPROCS, and every result is compared with the
+// serial kernel. The data are small integers, so every sum, product and dot is
+// exact and the comparison is equality; Exp runs the same per-element code in
+// both, so it is bit-identical too.
+func TestParallelStress(t *testing.T) {
+	dur := 1500 * time.Millisecond
+	if testing.Short() {
+		dur = 200 * time.Millisecond
+	}
+	op, og, omv := ParThreshold, GemmThreshold, matVecThreshold
+	ParThreshold, GemmThreshold, matVecThreshold = 64, 64, 64 // parallel from tiny sizes
+	defer func() { ParThreshold, GemmThreshold, matVecThreshold = op, og, omv }()
+	defer runtime.GOMAXPROCS(runtime.GOMAXPROCS(0))
+
+	stop := make(chan struct{})
+	var procs sync.WaitGroup
+	procs.Add(1)
+	go func() { // churn GOMAXPROCS under the callers
+		defer procs.Done()
+		r := rand.New(rand.NewSource(1))
+		for {
+			select {
+			case <-stop:
+				return
+			default:
+				runtime.GOMAXPROCS(1 + r.Intn(16))
+				time.Sleep(time.Duration(r.Intn(300)) * time.Microsecond)
+			}
+		}
+	}()
+
+	var calls atomic.Int64
+	var wg sync.WaitGroup
+	deadline := time.Now().Add(dur)
+	for c := 0; c < 6; c++ {
+		wg.Add(1)
+		go func(seed int64) {
+			defer wg.Done()
+			r := rand.New(rand.NewSource(seed))
+			ints := func(n int) []float64 {
+				a := make([]float64, n)
+				for i := range a {
+					a[i] = float64(r.Intn(7) - 3)
+				}
+				return a
+			}
+			for time.Now().Before(deadline) {
+				n := 1 + r.Intn(5000)
+				a, b := ints(n), ints(n)
+				got, want := make([]float64, n), make([]float64, n)
+				switch r.Intn(9) {
+				case 0:
+					AddP(got, a, b)
+					Add(want, a, b)
+				case 1:
+					MulP(got, a, b)
+					Mul(want, a, b)
+				case 2:
+					got[0], want[0] = SumP(a), Sum(a)
+				case 3:
+					got[0], want[0] = MaxP(a), Max(a)
+				case 4:
+					got[0], want[0] = Dot1DP(a, b), Dot1D(a, b)
+				case 5:
+					m := 1 + r.Intn(64)
+					k := n / m
+					if k == 0 {
+						continue
+					}
+					got, want = got[:m], want[:m]
+					MatVecP(got, a[:m*k], b[:k], m, k)
+					for i := 0; i < m; i++ {
+						want[i] = Dot1D(a[i*k:i*k+k], b[:k])
+					}
+				case 6:
+					outer, inner := 1+r.Intn(8), 1+r.Intn(700)
+					axis := n / (outer * inner)
+					if axis == 0 {
+						continue
+					}
+					src := a[:outer*axis*inner]
+					got, want = make([]float64, outer*inner), make([]float64, outer*inner)
+					RunAxisP(SumAxis, got, src, outer, axis, inner)
+					SumAxis(want, src, outer, axis, inner, 0, inner)
+				case 7:
+					ExpP(got, a)
+					Exp(want, a)
+				default:
+					SqrtP(got, b)
+					sqrtSIMD(want, b)
+				}
+				for i := range want {
+					if got[i] != want[i] && !(math.IsNaN(got[i]) && math.IsNaN(want[i])) {
+						t.Errorf("seed %d n=%d: [%d] parallel %v, serial %v", seed, n, i, got[i], want[i])
+						return
+					}
+				}
+				calls.Add(1)
+			}
+		}(int64(c + 2))
+	}
+	// A lost block hangs the call rather than corrupting it: fail by name.
+	finished := make(chan struct{})
+	go func() { wg.Wait(); close(finished) }()
+	select {
+	case <-finished:
+	case <-time.After(dur + 10*time.Second):
+		t.Fatal("a parallel call never returned: a block was lost")
+	}
+	close(stop)
+	procs.Wait()
+	t.Logf("%d parallel calls checked", calls.Load())
+}