Performance — go-ndarray vs NumPy

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.9.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 ppc64le section (POWER9, against the scalar code and OpenBLAS’s POWER9 kernel) and the loong64 section (Loongson 3C5000L, against the scalar code only). On the Zen 3:

• Faster: whole-array reductions (Sum/Mean/Max at 4 Mi: 2.3–4.5×), row reductions (SumAxis(1) 2.2×, MaxAxis(1) 1.7×), Exp from 256 Ki elements (2.4–5.6×), and MatMul against single-threaded OpenBLAS (3.7–5.9× from 512²); in a Workspace, x + y and sqrt(x*y + x) match or 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 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 and mat·vec against NumPy’s 16-thread BLAS (0.69× and 0.52× since v0.6.0, 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×).

On Apple silicon the GEMM reached parity with tuned BLAS at 1024² (OpenBLAS in an arm64 VM, single-threaded vecLib on an M4 Max) and beats the pure-Go gonum 4–10× (BENCHMARKS.md). On a POWER9 the VSX kernels added in v0.2.3 double the serial GEMM (7.9 → 15.8 GFLOP/s at 512²), which is 0.69× OpenBLAS 0.3.34 on one thread; MatMul at 1024² on 32 threads is at 0.56× (see ppc64le: VSX kernels). On a Loongson 3C5000L the LASX kernels added in v0.2.5 lift the GEMM micro-kernel from 7.0 to 27.5 GFLOP/s and the serial GEMM from 5.3 to 20.0 GFLOP/s at 512², against the scalar code; no NumPy or BLAS reference was available on that host (see loong64: LASX kernels). On a RISC-V SpacemiT X60 (VLEN 256) the RVV kernels added in v0.7.0 make Sum up to 5.8× faster and, with v0.7.1’s 6×8 micro-kernel, cut MatMul time by 31–45%, one core against the scalar code. On an IBM z17 (a z/VM guest on real hardware) the s390x kernels cut MatMul 1024² from 340 to 70 ms (v0.9.0) and Max by 92–96% (v0.9.1), likewise. Neither host had a NumPy or BLAS reference (see SIMD coverage). SIMD kernels exist on amd64, arm64, ppc64le, loong64, riscv64 and s390x (LASX, RVV and the s390x vector facility used only when the CPU reports them); the 32-bit targets run the pure-Go reference code; the s390x kernels are measured on a real IBM z17 (see SIMD coverage at the end).

How a pure-Go library can beat NumPy

NumPy’s elementwise and reduction ufuncs (+, *, sqrt, sum, max, …) are single-threaded SIMD C. So a pure-Go library wins them by combining two levers NumPy does not use for these ops:

1. Multicore — fan the contiguous loop across GOMAXPROCS goroutines above a size threshold (small arrays stay serial so scheduling never dominates).
2. SIMD via go-asmgen — for the sum reduction, whose scalar s += v loop the Go compiler will not auto-vectorize (it would change FP grouping), a hand-written NEON/SSE2 kernel with four independent accumulators recovers the lane parallelism NumPy’s C has.

A third lever — eliminating allocation — turned out to be the single biggest factor: the old elementwise path materialised two full broadcast copies and a destination per op (3×N). The same-shape contiguous fast path now passes the operands’ live backing slices to the kernel and allocates only the result.

NumPy’s matmul is BLAS, a different fight: go-ndarray answers it with its own panel-packed, cache-blocked GEMM + go-asmgen SIMD-FMA micro-kernel (the OpenBLAS/BLIS structure). With the by-lane-FMLA micro-kernel it reaches parity with tuned multi-threaded OpenBLAS 0.3.29 at n=1024 (≈0.99×, ~203 GFLOP/s) and ~0.97× at n=512; only small n=256 still trails (per-call overhead). See the matmul section below.

Test bench

• Machine: debian arm64 (Apple-silicon Tart VM), 4 vCPU.
• NumPy: 2.2.4. Elementwise/reduction rows use Debian reference BLAS 3.12.1 (irrelevant to those ufuncs). The matmul rows are measured against OpenBLAS 0.3.29 (pthread, multi-threaded) — the tuned BLAS — selected via update-alternatives --set libblas.so.3 …/openblas-pthread/libblas.so.3 and confirmed loaded by reading /proc/self/maps (libopenblasp-r0.3.29.so). This is the meaningful bar: the earlier reference-BLAS matmul comparison is kept only as a footnote.
• Go: 1.26.4, CGO_ENABLED=0, GOWORK=off.
• go-ndarray bench: go test -bench=. -benchtime=1s (bench_test.go), ns/op. Multi-core = default GOMAXPROCS=4; single-core = GOMAXPROCS=1.
• NumPy bench: docs/bench_numpy.py (matched ops/sizes), min-of-7 ns/op.

Reproduce: rsync the repo to the VM, then go test -bench=. -benchtime=1s . and python3 docs/bench_numpy.py.

Elementwise and reductions (ns/op; lower is better)

go × is NumPy ÷ go-ndarray: > 1 means go-ndarray wins. Two go forms are shown for the binary/sqrt ops:

• alloc — a + b / np.sqrt(a): returns a fresh array (NumPy x + y).
• into — a.AddInto(out, b) / a.SqrtInto(out): writes a caller-supplied buffer, the no-allocation form (NumPy np.add(x, y, out=z) / np.sqrt(x, out=z)). This is the apples-to-apples small-N comparison: it removes the per-op result allocation, which is the only thing NumPy beats us on at small n (see below).

(All numbers measured back-to-back in one VM session, 2026-06; go is -benchtime=1s, NumPy min-of-runs in docs/bench_numpy.py. Add/Sub/Mul/Div now run a SIMD kernel — packed SSE2 ADD/SUB/MUL/DIVPD on amd64, NEON VFMLA on arm64 — so even the alloc form wins ~2× at large n. Sum/Max are reductions with no result array, so they have no separate into row.)

Reading the table

• Add / Mul on large arrays: go-ndarray wins ~2×. Allocation elimination + multicore. Note even single-core large Mul edges out NumPy now (the kernel is memory-bandwidth bound and the old copies were the bottleneck).
• Sum: go-ndarray wins everywhere — up to ~5× on small, 2.4× on large. The SIMD 4-accumulator kernel already beats NumPy single-core (493 µs vs 506 µs at 4 M); multicore then roughly halves it again. Small sums are dominated by NumPy’s per-call Python/dispatch overhead, which Go does not pay.
• Sqrt: go-ndarray now WINS 1.75× at 4 M (was a ~1.3× loss). The fix was to stop routing Sqrt through Map’s func(float64) float64: that indirection blocked both the Go compiler’s FSQRTD intrinsic (arm64) and a packed SSE2 SQRTPD (amd64), making the old path pay a real function call per element (~2.7× slower single-core). A dedicated sqrt kernel behind the kernels seam — packed SQRTPD on amd64, the intrinsic-lowered scalar FSQRTD loop on arm64/others — plus the existing multicore fan-out is what turns the loss into a win. The kernel is bit-identical to a scalar math.Sqrt loop (validated per-arch and bit-for-bit against NumPy’s np.sqrt). arm64 now also has a packed NEON FSQRT V.2D kernel (two doubles per instruction; emitted by WORD until v0.5.2, VFSQRT since Go 1.27 names it), replacing the prior one-lane-at-a-time scalar FSQRTD loop. Small n=1 024: the alloc form loses to NumPy on allocation cost (as Add/Mul), but a.SqrtInto(out) wins 1.8×.
• Max: go-ndarray now WINS — 4.0× small, 1.34× at 4 M (was 0.62 / 0.31 losses). Two changes. (1) The scalar oracle was switched from if v > m (which silently ignored NaNs, diverging from NumPy) to the builtin max/min, which is NaN-propagating and lowers to the hardware FMAXD/ FMIND intrinsic — math.Max would also propagate NaN but is an un-intrinsified call ~12× slower in this hot loop. (2) The SIMD max/min reduction uses four independent accumulators to break the dependency chain (~3.6× over one), exactly like the sum kernel: packed MAXPD/MINPD + a NaN scan on amd64, register-unrolled builtin max/min on arm64/others. See the NaN convention note below.
• Small arrays (n≈1 024) for Add/Mul/Sqrt — the alloc form: NumPy wins; the into form: go-ndarray wins ~2×. The entire gap is the result allocation. make([]float64, 1024) measures ~900–1 080 ns on this VM’s allocator (it zeroes the 8 KiB the kernel then fully overwrites), while NumPy serves the temp from a cached free-list for ~150 ns. So a + b at n=1024 is allocator-bound (≈1 µs, of which the SIMD compute is only ~150 ns) and loses 0.4×; a.AddInto(out, b) — the no-allocation form, NumPy’s np.add(x, y, out=z) — removes that cost and wins 1.8–2.4× at the same size. This is a Go-allocator ceiling, not a kernel one: Go’s make always zeroes and offers no unsafe “give me raw bytes” escape, and returning pooled memory would break the caller-owns-the-result contract, so the allocating small-N form cannot match NumPy’s free-list — the *Into parity path is the supported way to hit/beat NumPy there (and at large n even the alloc form wins, because the kernel time then dwarfs the one allocation). (Sum and Max are the exceptions — their reduction kernels need no result allocation and win even serially.) Superseded for loops by Workspace, see below.

NaN convention for Max / Min (and the SIMD max kernel)

Max/Min are NaN-propagating: if any element is NaN the result is NaN. This matches numpy.max/numpy.min, Go’s builtin max/min, and IEEE-754 maximum (not the C fmax / IEEE maxNum “ignore NaN” rule); on signed zeros it returns +0 for max(-0,+0) and -0 for min(-0,+0), also matching NumPy. The previous if v > m scalar form silently skipped NaNs (and a leading NaN poisoned the scan), so it disagreed with NumPy; this is the corrected, documented semantics, and the entire stack — scalar oracle, the four-accumulator reducer, the amd64 MAXPD/MINPD kernel, and the multicore fold — is held bit-for-bit identical to it (NaN at any position → NaN; exact extreme otherwise), validated across every length residue and NaN placement in CI. The amd64 kernel gets this despite MAXPD’s own non-propagating NaN rule by carrying a parallel CMPPD-unordered OR-mask and forcing the result to NaN if any lane ever held one — so the value MAXPD leaves in a NaN lane is irrelevant.

Matrix multiply vs TUNED BLAS (OpenBLAS) — the hard frontier

This is the honest fight: go-ndarray’s pure-Go (+ generated asm) GEMM against OpenBLAS 0.3.29 pthread, decades-tuned multi-threaded assembly. Both run on the same 4-vCPU Apple-silicon VM, back-to-back. GFLOP/s = 2·n³ / time.

All rows measured back-to-back, interleaved, in the same VM session (the only fair way on a shared host whose load drifts), 4 samples each, go ÷ OB reported as the per-sample ratio (robust to the absolute time drifting between sessions).

Verdict: at n=1024 the pure-Go (+ generated asm) GEMM reaches PARITY with tuned OpenBLAS (~0.99×, individual samples straddling 1.00×), sustaining ~203 GFLOP/s — up from the prior kernel’s ~23 GFLOP/s and from this kernel’s own 0.76× before the by-lane-FMLA fix below. n=512 is ~0.97×; only small n=256, where fixed per-call overhead (panel packing + the 4-way goroutine fan-out) dominates the ~250 µs of compute, still trails at 0.67× — a small-matrix overhead ceiling, not a kernel-throughput one (the same kernel hits parity once the matrix is large enough to amortise the setup).

The fix that reached parity: the by-lane FMLA

The single change that took n=1024 from 0.76× to ≈0.99× was the arm64 micro-kernel’s inner FMA. The kernel broadcasts one packed A value into every column of its row; the prior version did this with MOVD (A double → GP register)

• VDUP (GP → all lanes of a vector) before a plain vector VFMLA — an extra instruction and a GP↔SIMD register-file round-trip per A value, every k-step. The fix uses the indexed-element double FMLA FMLA Vd.2D, Vn.2D, Vm.D[i] that every tuned arm64 dgemm uses: load the four packed A doubles of a k-step as two D2 vectors and have each row’s FMA read its A scalar straight from a lane — no GP detour, no VDUP. Go’s arm64 assembler cannot name this indexed form (“illegal combination … ELEM”), so it is emitted by its raw 32-bit WORD encoding (see fmlaElem in asmgen/arm64/gen.go, cross-checked against objdump); this is still pure Go assembly, CGO=0.

Measured on this Apple-silicon core: the L1-resident micro-kernel went ~44 → ~58 GFLOP/s/core (~1.28×), and the full GEMM tracked it to parity. (An earlier note on this page claimed the indexed form was slower on Apple silicon — that was wrong; it had never actually been measured against a correct encoding.)

Two further levers were prototyped and measured to NOT help on this core, so they were not shipped (documented here so they are not re-attempted blindly):

• Software pipelining (prefetch + double-buffer the next k-step’s A/B into shadow registers, rotate after the FMAs): ~2% on the L1 micro-kernel but 1–9% slower on the full GEMM — the M-series out-of-order window already hides the L1/L2 load latency, and the extra VMOV rotations + register pressure cost more than they save.
• Wider 6×8 tile (24 accumulators vs 16): same GFLOP/s as the 4×8 — the loop is FMA-latency-bound, not register-pressure-bound, so a wider tile adds no throughput.
• MC/KC/NC retuning: a 12-point sweep confirmed the shipped 256/256/512 is already the optimum on this core.

What shipped: a panel-packed, cache-blocked GEMM

The earlier register-blocked attempt read A and B straight from the source matrices; on power-of-two row strides its destination rows collided in L1 and it lost to the autovectorized ikj kernel at large n, so it was capped to n ≤ 224 and packing was prototyped but not shipped (the previous note here). This round implements the full OpenBLAS/BLIS structure and ships it for all sizes:

• Packing — A is copied into MR-tall row panels and B into NR-wide column panels in contiguous, unit-stride scratch (pooled, allocation-free per call). The micro-kernel then streams conflict-free memory regardless of the source stride — this is the precise fix for the power-of-two L1 set-conflicts that defeated the unpacked attempt. Edge tiles are zero-padded so the kernel always sees a full MR×NR block; a ragged right/bottom tile runs the same micro-kernel into a scratch tile, and only its valid corner is added to C.
• Cache blocking — an NC (columns) → KC (contraction) → MC (rows) loop nest keeps the packed B panel L2-resident and each packed A panel in L1. Defaults MC=256, KC=256, NC=512 (tuned on this VM).
• SIMD-FMA micro-kernel (go-asmgen) — NEON 4×8 on arm64 (16 D2 accumulators, the source of the win), AVX2/FMA 6×8 on amd64 (12 YMM accumulators, picked at run time by a CPUID probe; an SSE2 6×4 pair without FMA is the fallback), scalar 4×4 on the four arches without vector-double asm — which still gain from contiguous packed data + blocking + multicore.
• Parallelism — the M rows are split into MR-aligned bands, one per core, each writing a disjoint region of C.

The arithmetic is the ikj order dst[i][j] += a[i][p]·b[p][j], summed per KC block; for k ≤ KC that is the scalar oracle’s order exactly, and above it the block partials regroup the sum (a valid reordering, as in every blocked BLAS). At 128×128 the result matched NumPy’s A@B against OpenBLAS with max abs diff 0.0.

Packed GEMM vs the PRIOR kernel (why it ships)

Measured back-to-back in one binary under identical load (so the ratio is robust to the noisy VM); the prior kernel is the register-blocked/ikj path this replaces:

The packed kernel wins at every size, so it is shipped uniformly — no size-routing — and the old n ≤ 224 blocked/ikj split is removed.

Where the remaining gap is (small n only)

At n=1024 there is effectively no kernel-throughput gap left: the by-lane FMLA inner loop runs ~58 GFLOP/s/core, the full GEMM hits ~203 GFLOP/s, and go ÷ OB straddles 1.00×. The only sizes still short are small matrices, and there the gap is fixed per-call overhead, not compute:

• At n=256 the matmul is ~250 µs of real work, but every call still pays the full panel-pack pass and a 4-way goroutine fan-out. OpenBLAS amortises that setup far better (smaller pack cost, a thread pool that is already warm), so it leads 0.67× here even though the per-core FMA rate is the same. As n grows the setup is amortised away — n=512 is already 0.97×.

This is a small-matrix overhead ceiling. Closing it further would mean a serial-vs-parallel crossover tuned per size and a cheaper pack for tiny panels; it is not a Go-assembler or algorithm limit, and it does not affect the large-matrix parity result. (amd64 now has its AVX2/FMA kernel too; see the real-hardware section below.)

Footnote: vs reference BLAS

Against Debian’s single-threaded reference BLAS 3.12.1 (the prior comparison on this page), the packed GEMM wins by a wide margin at every size — but that is not a meaningful bar for a matmul claim, so the table above uses OpenBLAS.

Correctness

Every result is validated bit-for-bit against NumPy (shapes and values), not just timed — see the cross-check in the work log: Sqrt (bit-identical to np.sqrt via a uint64 view), Max, Min (exact, incl. NaN propagation), Add, Mul and the full MatMul (zero max-abs-diff against A @ A) match NumPy at n = 100 000 / 80×80.

The SIMD kernels are also validated against the scalar oracle per-arch in CI:

• Sqrt — SQRTPD (amd64) is held bit-identical to a scalar math.Sqrt loop across every length residue mod 8 and the IEEE edge cases (negatives → NaN, ±Inf, signed zeros, explicit NaN); sqrt is a single correctly-rounded operation, so bit-identity (not mere closeness) is the contract.
• Max / Min — the four-accumulator reducer and the amd64 MAXPD/MINPD kernel are held bit-identical to the builtin-max/min oracle across every length and with a NaN injected at the start, middle, end, and as the sole element (all must return NaN); ±Inf and signed zeros bit-match too.
• Add / Sub / Mul / Div — the elementwise SIMD kernels are held bit-identical to the scalar oracle across every length residue mod 8 and the IEEE edge inputs (signed zeros, ±Inf, NaN, extremes): each lane is a single independent correctly-rounded IEEE op (no reduction grouping), and on arm64 the FMA-against-an-exact-constant form (b+a·1, a−b·1, 0+a·b) rounds identically to the plain op, so bit-identity is the contract.
• Sum — validated to a tight relative tolerance (not bit-identity): its lane-parallel grouping is a valid reordering, the same kind NumPy’s pairwise summation uses, so closeness — not bit-identity — is the contract for a floating-point reduction.

Where go-ndarray still loses, and the plan

Allocation: Workspace (2026-10-04)

On real amd64 hardware the allocating forms were far worse than on the M4: on a Zen 3 (16 cores), x + y lost to NumPy 0.13–0.25× at every size up to 256 Ki elements, while AddInto won. Slicing copies, Concatenate and broadcasting adds (whose results also allocate) lost 5–9×. Measured causes, on M4:

• The garbage collector. A benchmark’s live heap is tiny, so a 128 KiB result per call triggered a GC cycle every ~22 calls (912 cycles in 20 000). With 256 MiB–2 GiB of live heap, which is closer to a real program, the overhead halves: 8.3 → 3.6 µs at 16 Ki, against 2.4 µs for AddInto.
• The rest: make zeroes memory the kernel overwrites anyway, and fresh memory is cold or page-faulted. GOGC=off does not help: it trades the cycles for page faults (9.4 µs).

A library must not set GOGC, and Go has no unzeroed allocation. What NumPy’s reference counting gets for free is scoped reuse. ND4J, the Java ndarray library, has the same problem and answers it with workspaces: an arena that results of one pass come from, recycled at the end of the pass. Workspace is that. ws.Use(x) binds an input, every result computed from a bound array is carved from the arena (64-byte aligned, capped so an append cannot spill), results stay bound so a chain stays in it, and Reset recycles everything, folding a pass that needed several blocks into one block. Sites whose kernel writes every element (elementwise, ufuncs, copies, broadcasts, concatenation, Where, Clip, scans, Take, Outer, mat·vec) skip the zeroing. Those that accumulate (GEMM, vec·mat, axis reductions) clear their slot.

Zen 3, 16 cores pinned, NumPy 2.5.3 single-threaded, best of two runs:

On M4 the same chain gains ×5 at 1 024 and ×1.6 at 4 Mi. Correctness: every operation is checked bit for bit against its heap result, including on passes whose recycled memory was first filled with NaN (a site that skipped zeroing it needed would leak the NaN; planting that bug in the GEMM made the test fail).

amd64 on real hardware (2026-10-04)

Every figure above was taken on Apple silicon. This round measured amd64 for the first time, on an AMD EPYC 7773X (Zen 3, GCC Compile Farm cfarm421), with go-ndarray and NumPy both pinned to the same 16 cores of one NUMA node (taskset -c 16-31; Go’s GOMAXPROCS follows the affinity mask). NumPy 2.5.3 with its bundled OpenBLAS 0.3.34 (DYNAMIC_ARCH, Haswell kernels on Zen 3), once with OPENBLAS_NUM_THREADS=1 and once with 16. go rows are the minimum of three runs interleaved with the previous version, so drift on a shared host hits both alike.

Four changes, each measured on its own before the next:

1. AVX2/FMA 6×8 micro-kernel (the BLIS Haswell shape). In L1 it runs at 43 GFLOP/s on one Cascade Lake core (16 flops per cycle at about 2.7 GHz, if the virtual machine’s unknown clock is near that); the SSE2 4×4 it replaces had no FMA and half the vector width.
2. Edge tiles through the micro-kernel. With MR=6, every 256-row MC block ends in a 4-row edge panel, and the scalar edge loop took 40% of one profile covering serial 64² and 512² products. A partial tile now runs the kernel into a scratch tile.
3. Prefetch of the C tile at kernel entry and of the A panel eight steps ahead: serial 768² went from 27.5 to 34.7 GFLOP/s (+26%). 90% of the time was already in the kernel, which ran at 28 GFLOP/s there against 43 with its panels in L1, so it was waiting on memory. A sweep of MC (48–252) and KC (128–512) moved nothing (±3%), so the blocking is unchanged.
4. B packed by all workers. MatMulP packs the shared B block once per (jc, pc) round, and one goroutine did it while the others waited: +30–70% in parallel. Below 64 Ki elements the launch costs more than the copy (M4: −8–15% at 96³–128²), so small blocks are still packed by one.

Inner gained the most because it used to materialise bᵀ before multiplying. The GEMM now reads any 2-D view through its strides while packing (the copy it makes anyway), so a transposed, sliced, reversed or broadcast operand costs no extra pass.

Still behind OpenBLAS with 16 threads (0.54–0.66× on large products, 0.3× below 256²). Where the parallel loss is, measured on the Zen 3:

• It is not the hardware: 16 independent serial 512² GEMMs, one per core, 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. 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. 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.
• Taller row bands (more reuse of each B micro-panel) do not help: interleaved over 5 runs they were worse at 1024² and equal at 2048².
• A 2-D tiling of C with private packing (no barrier, nothing shared) was built and measured, and is not shipped: ×1.24 at 1024² on the Zen 3, but ×0.62 at 512² and Inner, and ×0.70–0.93 at every size on M4, whose cluster L2 holds the shared panel well. A tile repacks (rows + cols)·k elements for rows·cols·k multiply-adds, and a strided pack costs far more per element than an FMA, so tiles under ~256×256 lose.

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) and loong64 (cfarm401). The s390x host (LinuxONE) did not answer; s390x stays covered by the qemu lane in CI.

Axis reductions (2026-10-04)

On the Zen 3 the axis reductions of a 1024×1024 matrix lost to NumPy 0.23–0.43×. Two causes, both measured:

• An axis-0 reduction ran on one core. It has a single outer slab, and the driver only split outer slabs; the comment justifying that said a column split would need a strided gather. It does not: a band of columns is a contiguous run of every row. RunAxisP now splits the columns into bands of at least 512 when there are fewer slabs than workers, and each row run is added with the elementwise SIMD kernel (exact per element, so the sum is still in axis order).
• A row reduction (axis 1, inner = 1) was a sequential scalar loop, latency-bound at one add per 4 cycles. Each row now goes through the SIMD sum/max/min. Max and min are exact; the sum is regrouped lane-parallel, as NumPy’s own pairwise sum is.

Zen 3, 16 cores pinned, minimum of three interleaved runs:

On M4: MaxAxis(1) ×6–8, SumAxis(0) ×1.9.

Broadcasting a row without materialising it (2026-10-04)

M + row (1024×1024 plus 1×1024) lost to NumPy 0.26× on the Zen 3 even with a Workspace: the broadcast path copied the row out to a full 8 MiB matrix before adding. When one operand has the full shape and the other, leading 1s dropped, is a suffix of it (a row, a 2-D mean subtracted from a 3-D stack), the second simply repeats in blocks of its own length, so binOp now streams the live slices block by block through the SIMD kernel (kernels.RepeatP), split across cores. Below 64 elements per block the old path is kept.

Exp (2026-10-04)

Exp went through Map(math.Exp), and math.Exp itself was the cost: 13.2 ns per element on one Zen 3 core, 3.9 ns on M4 (the func-pointer call added only 15–20%). It is now a port of the double-precision exp of Arm’s optimized-routines (Szabolcs Nagy; glibc’s exp since 2.28; MIT licence): x = k·ln2/128 + r, 2^(k/128) from a 128-entry table as scale·(1+tail), exp(r)−1 from a degree-5 polynomial. The common case is written out in the loop, because the function is too large to inline and a call per element cost a third of the time.

• Accuracy: worst 0.504 ULP over 80 000 inputs across the whole finite range, the region near 0, the overflow edge and the subnormal range, measured against a 300-bit reference (math/big, ln 2 from its atanh series), against Arm’s documented 0.509; the committed test samples 20 000 of them and finds 0.501. math.Exp measured 0.879 ULP on the 80 000 inputs on arm64 (0.841 on the test’s 20 000). The table was not trusted as copied: TestExpTable re-derives all 256 entries from 2^(k/128) at 300 bits.
• A Go bug it removes, reported as golang/go#81995: on amd64, math.Exp returns +Inf for 709.436139303104 ≤ x < 709.782712893384, although the result is finite there (NumPy, glibc and mpmath: exp(709.5) = 1.3549863193146328e+308). Its assembly rounds k = x·log2(e) to 1024 from x = 1023.5·ln 2, and the biased-exponent check treats that as overflow although the reduced factor is below 1. Reproduced with Go 1.26.8 and 1.27.1 on AMD Zen 3, Intel Cascade Lake and Intel Haswell, with and without the FMA path; the assembly is unchanged on master. go-ndarray’s Exp inherited it until v0.2.0; TestExpTopOfRange pins the fix.
• Just above ln(2^-1075), exp rounds up to the smallest subnormal, as glibc and NumPy do. Go’s arm64 math.Exp returns 0 there (amd64 does not); that is one unit of the last subnormal place, not a contract violation.

Speed, kernel only, 16 Ki elements: M4 30 µs (math.Exp 65 µs); one Zen 3 core 62 µs (math.Exp 215 µs, NumPy single-threaded 78 µs). Whole Exp on Zen 3, 16 cores: 1 Ki ×2.7, 16 Ki ×3, 256 Ki ×1.6, 4 Mi ×1.65; from 256 Ki on it is 2–5× NumPy. A SIMD version (AVX2 gathers for the table) is the next step.

Dot and mat·vec (2026-10-04)

dotRange, under Dot (1-D), MatVec and Dot1DP, was a four-chain scalar loop (gc does not vectorise it): on one Zen 3 core 2^20 elements took 464 µs, 2× NumPy’s single thread. Two generated kernels replace it: dotFMA on amd64 (four YMM accumulators, b read straight from memory by VFMADD231PD; two loads per FMA make it load-bound, which four chains cover), gated on the FMA probe, and dotNEON on arm64 (eight D2 accumulators for the four FP pipes of an Apple core). The other targets keep dotRange. The lane-parallel sum is a regrouping, held to the n·ε·Σ|aᵢbᵢ| bound and exact on integer data.

NumPy’s multi-threaded dot still wins on this part (24 µs). The go-ndarray version stops scaling at about 90 GB/s from 4 cores on. A likely cause, not verified, is cache affinity: OpenMP keeps each thread on one core, so it rereads its own chunk from its own L3 slice on every repetition, while goroutines move between core complexes.

Measured 2026-10-07 (POWER9 bare metal, cfarm29, one thread per core): the cache-affinity guess does not hold. The pattern is the fork/join of parallelFor, the same as the GEMM’s (see Parallel GEMM: start once, wait on work). A 2^20 dot took 212 µs on 4 cores and 553 µs on 8, close to the 688 µs of one core. With 8 workers, only 2.4 CPUs were busy (perf stat). Three replacements for parallelFor were built and measured against it (8–10 interleaved rounds); none is shipped:

The second row keeps exactly the static chunks of the original and still loses at 4 workers, so affinity is not the cause. perf stat places the loss in scheduling instead: with the caller working, 2.1 of 4 CPUs were busy against 3.4 for the original. Its goroutines wait in the run queue of the caller’s processor until another one is woken to steal them. With the caller blocking, the 4-worker case recovers, but the 8-worker gains go with it. The code is on the experiment/parallel-for-scheduling branch.

Solved in v0.6.0 with a persistent helper pool. The cost being paid was waking threads, so the fix is to stop putting them to sleep between back-to-back operations, as OpenBLAS and OpenMP do. Up to GOMAXPROCS-1 helper goroutines are started once. After an operation they poll for the next one for 200 µs, then park. The work keeps owned runs with stealing (row two of the table above). Reductions now use four partials per worker: with one each, a helper that saw the job late had its only block taken by the caller, who then summed two in a row (Sum of 4 Mi on 4 cores: 2.8 of 4 CPUs busy). On POWER9 bare metal, each case in a fresh process, 6 interleaved rounds, against v0.5.2 (median):

A first combined run, with every benchmark in one process, also showed serial 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).

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):

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)

Log went through Map(math.Log): 11.6 ns per element on a Zen 3 core against NumPy’s 4.2. It is now a port of the double-precision log of Arm’s optimized-routines (the companion of Exp, also glibc’s since 2.28; MIT licence): x = 2^k·z, z split into 128 subintervals, log(x) = k·ln2 + log(c) + log1p(z/c − 1) with 1/c and log(c) from a table and a degree-6 polynomial; inputs within 2⁻⁴ of 1 take a degree-12 polynomial. z/c − 1 is one math.FMA (exact on every target), so the second table the non-FMA variant needs is not used. Log10 is log(x)·(1/ln10), the formula of math.Log10, on this log.

• Accuracy: worst 0.508 ULP over the test’s 10 000 inputs (whole range, near 1, subnormals), against a 300-bit reference (log m = 2·atanh((m−1)/(m+1))); math.Log measured 0.724 on the same inputs on arm64. TestLogTable checks each of the 128 entries against the properties log_data.c documents (1/invc inside its subinterval, logc = log(c) within the rounding of invc, 0x1.8p9 + logc exact).
• A second Go bug it removes, known since 2022 as golang/go#56600 (fix pending in CL 448216): on amd64, math.Log is wrong for every subnormal input. math.Log(5e-324) returns −709.09 instead of −744.44, and math.Log(1e-310) −709.09 instead of −713.80 (NumPy, glibc: −744.44, −713.80). Its assembly inlines Frexp with bit masks that assume a normal input. Reproduced with Go 1.26.8 and 1.27.1 on three x86 CPUs. Log and Log10 inherited it; TestLogSubnormal pins the fix. (Log2 already normalised through math.Frexp and was right.)

Speed, 16 Ki elements, one core: M4 37 µs (math.Log 72), Zen 3 73 µs (math.Log 175, NumPy 69).

ppc64le: VSX kernels (2026-10-04)

Go’s ppc64 assembler has the VSX loads, stores and permutes but no vector double arithmetic, so ppc64le ran the scalar code. go-asmgen v0.13.0 encodes xvadddp, xvsubdp, xvmuldp, xvdivdp, xvmaddadp, xvmaxdp, xvmindp and xvsqrtdp as WORDs, pinned against GNU as 2.44 and run bit-identical to Go’s arithmetic and math.FMA on POWER8 and POWER9. With them, go-ndarray has VSX kernels for sum, dot, sqrt, add/sub/mul/div and the GEMM.

The GEMM micro-kernel (8×8, 32 accumulators) follows OpenBLAS’s POWER9 dgemm: no splats. Each A pair and its swap multiply each B column pair, so one accumulator holds a diagonal (C[r][c], C[r+1][c+1]) and another the anti-diagonal; the store rebuilds the rows with xxpermdi. On one POWER9 core, measured in one binary: 13.3 GFLOP/s with lxvdsx splats, 14.8 with xxpermdi splats, 15.4 with those unrolled by two, 16.9 with the diagonal scheme unrolled by two (shipped).

POWER9 (cfarm29, 8 cores × SMT4), against the scalar code (best of 3–4 interleaved runs):

At 4 Mi elements sum and add are memory-bound (×1.1). Against OpenBLAS 0.3.34 (its POWER9 kernel, loaded through scipy-openblas64 on the same host; the ppc64le NumPy wheel there has no BLAS): one thread 22.8 GFLOP/s at 512², so the serial GEMM is at 0.69×; 32 threads 13.1 ms at 1024², so MatMul is at 0.56×. OpenBLAS unrolls further and uses POWER9-only lxv loads; Go’s ppc64le baseline is POWER8.

loong64: LASX kernels (2026-10-04)

Go’s loong64 assembler has LASX vector float64 add/sub/mul/div, square root, 256-bit loads and the broadcast load (spelled as an arrangement, XVMOVQ off(R), X.V4, which assembles to xvldrepl.d), but no fused multiply-add; go-asmgen v0.14.0 encodes xvfmadd.d as a WORD (pinned against GNU as 2.43, bit-identical to math.FMA on a Loongson 3C5000L; transitional until cmd/asm names it). go-ndarray uses them for sum, dot, sqrt, add/sub/mul/div and an 8×8 GEMM micro-kernel (16 four-lane accumulators, A broadcast by xvldrepl.d). (v0.2.5 said Go had no broadcast load and emitted it as a WORD too; it has one, and v0.2.6 emits the mnemonic, the same machine code.) LASX is not on every LoongArch CPU, so the kernels run only when the kernel reports it in AT_HWCAP (/proc/self/auxv), with the scalar code as the fallback; the parsing is a pure function tested on every target.

Loongson 3C5000L (cfarm401, 32 cores), against the scalar code, best of 3 interleaved runs (the host’s load average of ~146 is processes stuck in uninterruptible sleep; its CPUs measured 100% idle):

At 4 Mi elements in parallel, sum and sqrt are memory-bound (parity); the 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):

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.

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:

(× = 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.

Sums of the other dtypes (2026-10-18)

Measured first, because nothing had: BenchmarkDTypeSum and BenchmarkDTypeAdd (benchmarks/dtype_bench_test.go) run the same operation on float64, float32, int64 and int32. On POWER9 (cfarm29, one thread pinned to the idlest CPU by hostload -pick 1, the host 3% busy before and after), a float32 Sum of 2²⁰ elements took 6.85 ms against 0.23 ms for float64, 6.8 ns an element, while Add was only 1.5× slower than float64. A missing SIMD kernel explains a factor of 4 to 8, not 30.

The cause was the reduction loop. Every non-float64 reduction walked the [outer][axisLen][inner] layout with the reduced axis outside and the kept elements inside, so a whole-array sum (inner = 1) took a slice and ran a one-iteration loop per element. A contiguous run is now folded directly, as NumPy’s add and multiply reductions do: the first element, then the rest added to it, pairwise for float32 and the complex dtypes (NumPy’s FLOAT/CFLOAT_pairwise_sum, loops_utils.h.src at v2.5.3: eight accumulators, blocks of 128, halved above), in order for the integers (exact in any order) and for products (NumPy multiplies in order).

Six runs each, old and new binaries alternating, benchstat p = 0.002. The float32 sum is also more accurate: on 2¹⁶ values of mixed sign and magnitude its error against the exact sum (math/big) fell from 75 to 0.87, and 2²⁴ followed by 1000 ones no longer saturates at 2²⁴ (sum_pairwise_test.go).

The narrow integer sums still converted the whole array to int64 (NumPy’s accumulator) before adding: an extra pass and 8 bytes an element. Since v0.14.0 a sum or product of bools or integers narrower than 64 bits along a contiguous axis widens each element as it is read (foldWidened, same arithmetic, wrap-around included): int32 Sum of 2²⁰ elements 2.872 → 1.548 ms (−46%), 2¹⁴ −51%, with float64, float32 and int64 unchanged as controls (POWER9, same protocol, p = 0.002).

float32 sums on SSE2 and NEON (2026-10-18)

NumPy keeps eight accumulators in its pairwise sum so that the inner loop can be vectorised “without changing summation ordering”. kernels.AddBlocks8F32 is that loop: blocks of eight float32 added lane by lane, one IEEE addition per lane per block, in order. It is SSE2 on amd64 (ADDPS on two XMM registers) and NEON on arm64 (VFADD on two .S4 registers), and the scalar loop elsewhere, so float32 sums are unchanged bit for bit. The kernel is tested bit for bit against the scalar loop (NaN, the infinities, −0, subnormals, overflow), and fenced with PROT_NONE pages on both operands. A kernel adding from the wrong register fails the first test, and one reading a block too many faults in the fence.

Measured twice. First with the bench workflow on the hosted runners (base and head alternated ten times in fresh processes on one thread, benchstat, p ≤ 0.002), because the development Mac ran at a load of 97 and most cfarm hosts refuse their keys since the host CA was revoked. Then on real, idle hardware reached through a known_hosts file whose keys were confirmed over a second network path: cfarm424 (Neoverse N1, 0.0% busy) and cfarm420 (Zen 3, EPYC 7773X, 1.7% busy), same protocol, pinned to the idlest CPU. On those, float32 Sum was −39% / −32% (N1) and −43% / −41% (Zen 3) at 2¹⁴ / 2²⁰. int32 was −15% / −11% on the N1, and on the Zen 3 within the noise (±10% at that host’s variance). The controls were unchanged.

The integer rows come from a defect the first run exposed. With the kernel added, int32 Sum was 48–56% slower on amd64 although its code had not changed: the disassembly of both binaries showed the same loop, byte for byte, moved 224 bytes and now straddling a 64-byte line. The loop also reloaded the operation from the stack every element. The widening folds now use one loop per operation and four accumulators (exact: integer addition and multiplication wrap modulo 2⁶⁴, so the order changes nothing). They also consume the slice instead of indexing i+1…i+3, which had left four bounds checks per iteration in the amd64 code; now there are none.

riscv64 and loong64 (2026-10-18). RVV: VSETIVLI $8, E32, M2 gives vl = 8 on every VLEN ≥ 128, and VFADDVV adds a block to the two-register group. LSX: the eight accumulators sit in two 128-bit registers, and VFADD.S (VADDF) adds a block. Each runs behind the feature bit its float64 kernels use. Both were run on the real hardware: the kernel tests passed there, and a sabotaged kernel failed them, so the kernel path is the one taken. Same A/B protocol, idle hosts:

The first loong64 kernel used one 256-bit LASX register for all eight accumulators, and was 26–47% slower than the scalar loop. Every block’s addition waited on the previous one, a single dependency chain, where the scalar loop has eight. Two registers give two chains, as SSE2 does on amd64.

ppc64le (2026-10-18). Go’s assembler has no xvaddsp, so go-asmgen emits it as a WORD. The encoding is pinned against GNU as 2.44 and llvm-mc on a POWER9 (go-asmgen/asmgen#27). The cmd/asm spelling is CL 846285, and go-asmgen’s registry check assembles every entry with that CL applied. The eight accumulators sit in two VSX registers, loaded and stored with LXVW4X/STXVW4X, so every word keeps its lane. VSX is the ppc64le baseline, so there is no feature check. The kernel tests passed on the POWER9, and a sabotaged kernel (one xvsubsp) failed them.

s390x (2026-10-18). Go’s assembler lacks vfasb (vector-enhancements facility 1, z14), which CL 846286 adds. go-asmgen v0.19.0 emits it as WORD+BYTE+BYTE, the 6 bytes pinned against GNU as (go-asmgen/asmgen#29). The eight accumulators sit in two vector registers; VL/VST keep memory order on this big-endian machine. The kernel runs behind the useVXE gate of max/min, so a z13 takes the scalar loop. The kernel tests passed on the z machine, and a sabotaged kernel (one vfssb) failed them. Same A/B protocol:

The float64 control at 2²⁰ moved +2.0% with its kernel byte-identical: the change moved sumVX by 32 bytes (0x18f920 → 0x18f940). It is attributed to layout, not proven.

Every 64-bit target now adds the float32 block with SIMD.

float32 elementwise arithmetic (2026-10-18)

Measured against NumPy 2.5.3 on one POWER9 core, float32 Add was 7.1× slower at 2¹⁴ (33.8 µs against 4.74 µs) and 5.0× slower at 2²⁰. A profile put 57% of the time in the generic loop binNum[float32]: a scalar add per element, with two bounds checks, because nothing told the compiler that x and y were as long as d.

Two changes:

• Bounds checks. The elementwise loops reslice their operands to len(d) once, so the compiler drops the per-element checks (verified with -d=ssa/check_bce). On the POWER9, float32 and int32 Add gain 7–12%.

• Kernels. kernels.Bin32 adds, subtracts, multiplies or divides 16 floats per iteration on every 64-bit target:

  target	instructions
  amd64	SSE2 ADDPS…
  arm64	NEON .S4
  ppc64le	VSX xvaddsp…
  s390x	z14 vfasb…
  loong64	LSX vfadd.s…
  riscv64	RVV E32, strip-mined

  Each lane is one IEEE operation, so the result is the scalar loop’s bit for bit. The kernel tests (special values, every length to 70, in place) and the PROT_NONE fence pass on Zen 3, Apple M, POWER9, IBM Z, SpacemiT X60 and Loongson 3C5000L. A sabotaged kernel (one vfsub in add32NEON) failed them.

AddInto and the other *Into methods also wrote a typed result into a new array and then copied it into out. When out has the result dtype they now write into it directly.

POWER9, one core, against NumPy 2.5.3:

float64 AddInto against np.add(out=): 5.70 against 7.75 µs (2¹⁴), and 724 against 768 µs (2²⁰).

float32 survey against NumPy, and Sqrt (2026-10-20)

The rest of the float32 surface, on the same POWER9 core against NumPy 2.5.3:

Read each row against the one below it:

• float32 Max has no kernel anywhere, and on ppc64le neither has float64: xvmaxdp does not propagate NaN as NumPy does.
• float32 Sqrt called math.Sqrt through a function value for every element.
• float32 MatMul is 9× slower than our own float64 MatMul: float32 has no packed GEMM.

Sqrt is fixed here. kernels.Sqrt32 runs SQRTPS, NEON .S4, xvsqrtsp, vfsqsb (z14), LSX vfsqrt.s or RVV E32. Each is the correctly rounded float32 square root, so float32(math.Sqrt(float64(x))) bit for bit: float64 carries more than twice float32’s precision.

The kernel tests (special values, every length to 70, in place) and the PROT_NONE fence pass on Zen 3, Apple M, POWER9, IBM Z, SpacemiT X60 and Loongson 3C5000L. Measured with fleettools’ abbench (base and head alternated, 8 rounds, pinned core, --max-busy 15):

The controls (float64 Sqrt, Add, AddInto) are flat. The gap left to NumPy is the allocation, the open defect below.

float32 Max/Min (2026-10-26)

Two changes:

• Every architecture. float32 Max/Min converted the whole array to a new float64 array and reduced that. They now read the float32 elements (kernels.Max32/Min32), with four independent chains of Go’s max/min. The value is the same, since the conversion is exact.
• amd64 and arm64: vector kernels.
  • amd64: MAXPS/MINPS with a CMPPS-unordered NaN mask. A zero result goes back to the scalar loop for its sign, the v0.19.1 rule.
  • arm64: NEON VFMAX/VFMIN on .S4, which propagate NaN and order −0 below +0 themselves.

Each is bit for bit one sequential pass of Go’s max/min. That is tested with NaN, real signed zeros and infinities at every position up to length 70, and a sabotaged kernel fails it on both architectures. Measured with abbench (alternated, 6 rounds, pinned core):

On the POWER9 float32 Max is still about 7× slower than NumPy, because Go’s NaN-correct max is scalar there.

Still to do:

• ppc64le: done in v0.24.0, below, without a registry entry.
• loong64: done in v0.25.0, below, without a registry entry.

loong64 Max/Min, and negative NaNs (2026-10-26)

Go’s loong64 assembler has no vector float max, min or compare, so the LSX kernels work entirely in integer arithmetic:

• Order: t = x XOR ((x >>a N−1) >>l 1) flips the magnitude bits of the negatives. The signed order of t is the float order, with −0 below +0 (Go’s rule), and the map is its own inverse.
• Extremum: max(acc, t) = acc XOR ((acc XOR t) AND (acc < t)), with vslt.
• NaN: a lane is NaN iff (x AND abs) >u +Inf's bits (vsltu), OR-accumulated.

What the sabotage run found. A sabotaged kernel, with the NaN mask of max32LSX’s first lane dropped, passed every test in CI’s loong64 lane. In that integer order a positive NaN sorts above +Inf, so a max finds it without any NaN test. Every NaN in the tests was math.NaN(), which is positive. A negative NaN sorts below everything, and only the NaN test catches it. That is no corner case: x86’s default NaN (from 0×∞ or ∞−∞) is negative.

The extremum tests now cover both signs, for float32 and float64 on every architecture. With them the sabotage failed CI’s loong64 lane (special Max32 n=16: -0, sequential NaN), and every other architecture’s kernels passed.

Validated under emulation only. The kernels pass on qemu la464 (LSX and LASX, required since #67). The Loongson they would be measured on (cfarm401) accepted TCP but never sent its SSH banner on 2026-10-26, so their speed is not measured yet.

ppc64le Max/Min (2026-10-26)

xvmaxdp/xvmindp (and the sp forms) do not propagate NaN as Go and NumPy do, and Go’s assembler has no VSX float compare to build a NaN mask from. The mask is built in integer arithmetic instead:

• A lane is NaN iff its magnitude bits (x AND 0x7FF…F) exceed +Inf’s bits as an unsigned integer. vcmpgtud/vcmpgtuw test that on the VMX registers, against constants loaded from a read-only DATA table, and the results are OR-accumulated.
• The kernel works only vertically: four accumulators folded into one vector stored to memory, plus the NaN flag.
• Go folds the lanes and the tail with its max/min. A zero result goes to the scalar oracle for its sign, as on amd64.

Max is now faster than NumPy on the POWER9 for both widths, and the Sum controls are flat. A sabotaged kernel, with one of the four NaN-mask ORs removed, fails the float32 and float64 oracle tests there.

• s390x: done in v0.23.0, below.

s390x float32 Max/Min, and the z14 float32 kernels measured (2026-10-26)

The s390x kernels use VFMAXSB/VFMINSB in mode 1 (Java’s Math.max/min): a NaN operand gives NaN, and −0 orders below +0. Four four-lane accumulators are folded by an 8-byte VSLDB rotation, then lane 1 onto lane 0, behind the useVXE gate (z14).

The IBM Z (LinuxONE) had looked unreachable since 2026-10-25, but the access had not expired. After the macOS upgrade the SSH agent held no key, and the id_ed25519 the server accepts needs a passphrase. With the key loaded again, everything was measured on the real machine (abbench, 6 rounds, alternated, pinned core):

The float64 controls are flat in both runs.

On the machine, a sabotaged kernel (one VFMINSB lane in max32VXE) fails the sequential-oracle test. The kernel suite passes with NDARRAY_REQUIRE_VXE=1, and so does the root suite. Since #68 the s390x qemu lane in CI also runs VX and VXE, on qemu’s max CPU.

• riscv64: done in v0.22.0, below.

riscv64 Max/Min and a sum defect (2026-10-26)

The RVV kernels for float64 and float32 work in three steps:

• Accumulate: each strip is folded into a register group with vfmax.vv (or vfmin.vv). The tail-undisturbed policy keeps a short last strip from touching the other lanes.
• Count NaNs: RISC-V’s fmax returns the non-NaN operand, so each strip also counts its NaN lanes (vmfne.vv x, x, then vcpop.m), and a non-zero count gives NaN.
• Fold: vfredmax.vs reduces the group at the end. fmax already orders −0 below +0.

On a SpacemiT X60 (abbench, 6 rounds, alternated):

Running the root suite on that real X60 found a defect: MeanAxis(1) of the int64 row [MaxInt64, MinInt64, −3] was 0, NumPy’s −1.

• Cause: the strip-mined sum kernel took even 3 elements as one vector, and its unordered reduction added MaxInt64 + (−3) first. In float64 that rounds back to MaxInt64 and cancels.
• Fix: below 8 elements the sum is now sequential, the order NumPy and the other targets’ kernels (whose scalar tails take n < 8 whole) use there.
• Why CI never saw it: its riscv64 lane runs under qemu, whose default CPU has no V, so no RVV kernel had ever run in CI. The lane now enables V (QEMU_CPU=rv64,v=true,vlen=256). It also sets NDARRAY_REQUIRE_RVV, so a CPU without V fails the lane instead of passing on the scalar paths.

float32 MatMul on the float64 GEMM (2026-10-26)

float32 MatMul ran a generic triple loop, 9× slower than our float64 GEMM. It now runs on the float64 GEMM:

• the operands are widened, which is exact and costs O(n²);
• the products are accumulated in float64;
• each element is rounded once to float32.

That is at least as accurate as a BLAS sgemm, which accumulates in float32, and every architecture gets the packed, cache-blocked SIMD kernels at once. With integer-valued operands the result is exact; with random ones it agrees with the float32 triple loop to float32 precision (TestMatMulFloat32).

float32 is now within the widening of float64: 86.7 against 54.0 µs at 64², 2.56 against 2.25 ms at 256² on the POWER9. A native float32 GEMM, twice the lanes per register and half the bytes, could go further; it is the next step for products.

A native float32 GEMM on amd64 and arm64 (2026-10-26)

The packed GEMM driver (packing, cache blocking, the parallel bands) is now generic over float32 and float64. Each element type brings its micro-kernel, its tile and its pack-buffer pools. float32 gets two kernels, generated by go-asmgen:

• arm64: gemmMicro8x8S, a NEON 8×8 tile. It has 16 accumulators and uses the by-element FMLA Vd.4S, Vn.4S, Vm.S[i], encoded as a WORD and checked against go tool objdump.
• amd64: gemmMicro6x16FMA, the BLIS Haswell sgemm shape: a 6×16 tile, VBROADCASTSS and VFMADD231PS into 12 YMM accumulators. It runs only where the CPU has FMA.

float32 MatMul (and Dot, which goes through it) now accumulates in float32, as NumPy’s sgemm does. Without a kernel (no FMA, or the other four targets) it still widens to the float64 GEMM. The tests check it against a float32 ikj oracle, bit for bit on integer data, in every block configuration, on a transposed and offset B view. On random data they bound it by k·eps32·Σ|ab|. A sabotaged kernel (a wrong A lane on arm64, a wrong A row on amd64) failed every job that runs one.

Profiling the first version at 64² found packA at 33% of the time, ahead of the micro-kernel. It gathered A column by column and paid its loop and bounds checks every MR elements. A row-major A is now read one contiguous row at a time, which also speeds up float64.

Against NumPy 2.5.4 (OpenBLAS, OPENBLAS_NUM_THREADS=1, the same core):

v0.26.1 packs full B panels with one fixed-size array move per row (8 or 16 elements, lowered to vector moves), where packB had copied element by element. Profiling had put it at 21% of a 64² float32 product.

Reading A in place (2026-10-27)

v0.27.0 no longer packs a row-major A (unit column stride, positive row stride). Unpacked-A kernels read it in place: the “rv” gemmsup kernels of BLIS (kernels/haswell/3/sup/bli_gemmsup_rv_haswell_asm_s6x16m.c). B stays packed, so its last panel is still zero-padded. There is one kernel per tile height, so a band’s bottom rows never read past A’s last row.

• amd64 (FMA): each A value is broadcast straight from its row through a scaled index (VBROADCASTSS (SI)(R8*2)), as BLIS does.
• arm64: the same idea, but broadcast loads (VLD1R) made it 26–78% slower on a Neoverse N1, apparently because they compete with the FMAs for the vector pipes. What shipped instead loads 4 floats (2 doubles) of each row with one plain vector load, and the next k steps read them by lane with the by-element FMLA.

Against NumPy 2.5.4 (OpenBLAS, one core), the ratio ours / NumPy:

Open defect: products are still up to 1.2× slower than OpenBLAS on one core. The exception is float64 at 1024² on Zen 3, where we are 3% faster. The kernels for ppc64le, s390x, riscv64 and loong64 are still to write.

Open: the allocating operations. Add is still 1.5–3.2× slower than NumPy’s a+b, for float64 as well. The kernel is no longer the difference, as AddInto shows; the allocation is. It splits like this, measured 2026-10-26 (float64 Add, 2¹⁴ elements, one POWER9 core, 3 runs each):

Neither part is the library’s to remove. Go zeroes every new allocation, and the GC’s pace belongs to the program. With a live heap this small, a GC starts every few dozen results, so the benchmark is the worst case. A program that holds more data runs the GC far less often per operation.

What does remove it:

• Reuse buffers: the *Into methods with a reused out, or a Workspace for chains of operations.
• Tune the GC: GOGC or GOMEMLIMIT for a program dominated by short-lived results.

SIMD coverage

• amd64 (SSE2) ships hand-vectorized sum (4-accumulator ADDPD), sqrt (packed SQRTPD), max/min (MAXPD/MINPD + CMPPD NaN scan), the elementwise add/sub/mul/div (packed ADD/SUB/MUL/DIVPD, 8 doubles/iter + scalar tail), all at the GOAMD64=v1 baseline, and the GEMM micro-kernel: gemmMicro6x8FMA (6×8, 12 YMM accumulators, VBROADCASTSD + VFMADD231PD, C and A prefetched) when go-asmgen’s CPUID probe reports FMA with the OS saving YMM state, else a pair of SSE2 6×4 MULPD+ADDPD tiles over the same packing. Generated by go-asmgen and validated per-arch in CI; the FMA kernel on real Cascade Lake and Zen 3 hosts. float32 products have their own FMA kernel, gemmMicro6x16FMA (6×16, VFMADD231PS).

• arm64 (NEON) ships hand-vectorized sum, dot, packed sqrt (VFSQRT), the elementwise add/sub/mul/div (VFADD/VFSUB/ VFMUL/VFDIV), max/min (VFMAX/VFMIN, four D2 accumulators from 8 elements up) and the GEMM micro-kernel (gemmMicro4x8: a 4×8 tile, 16 D2 accumulators, using the by-lane FMLA FMLA Vd.2D, Vn.2D, Vm.D[i]), with a float32 one, gemmMicro8x8S (8×8, FMLA Vd.4S, Vn.4S, Vm.S[i]). Only the indexed-element FMLA is still a raw WORD (16 of them): Go’s assembler does not name it yet (pending as CL 764224).

  ⚠ Corrected in v0.5.2. Up to v0.5.1 this paragraph said Go’s arm64 assembler had no vector FSQRT, no plain vector FP add, and no vector form for div, so add/sub/mul went through an exact FMA against 1.0, sqrt was a WORD, and div and max/min stayed scalar. Go 1.27 (the floor since v0.5.0) names all of them; the ISA always had them. Measured on one core, interleaved against v0.5.1 (8 rounds, median; Sum, SqrtInto and AddInto, unchanged, measured at 1.00× as the control):

  one core	DivInto	Max	MulInto
  Apple M (local)	1.34–2.19×	1.87–2.24×	1.21× at 1 Ki, 1.00× above
  X-Gene (cfarm185)	1.00×	1.18–1.41×	1.08–1.35×

  X-Gene’s vector divide is no faster than its scalar one. The sum, dot and GEMM C-accumulate still fold with an exact VFMLA against 1.0, written before VFADD existed in Go; the result is the same.

• ppc64le (VSX), since v0.2.3: sum, dot, sqrt, add/sub/mul/div and an 8×8 GEMM micro-kernel; max/min stay scalar (the ISA’s xvmaxdp NaN rule is not NumPy’s). See the ppc64le section above.

• loong64 (LASX), since v0.2.5, when the kernel reports LASX in AT_HWCAP (not every LoongArch CPU has it): sum, dot, sqrt, add/sub/mul/div and an 8×8 GEMM micro-kernel; max/min stay scalar. See the loong64 section above.

• riscv64 (RVV), since v0.7.0, when the kernel reports V in AT_HWCAP (the vector extension is optional): sum, dot, sqrt and add/sub/mul/div, in vector-length-agnostic loops (VSETVLI strip-mining, LMUL=8, no scalar tail), every instruction a Go mnemonic. max/min stay scalar: RVV’s vfmax/vfmin implement IEEE 754-2019 maximumNumber, which drops a NaN that NumPy propagates. Since v0.7.1 the GEMM micro-kernel is a 6×8 RVV tile (see below). On a SpacemiT X60 (cfarm95, VLEN = 256), one core pinned to an idle CPU, 6 interleaved rounds against the scalar code (median):

  	1 Ki	16 Ki	256 Ki
  Sum	5.81×	2.51×	1.84×
  AddInto / MulInto	3.2–3.3×	1.11–1.13×	1.07–1.12×
  DivInto	2.27×	2.30×	2.28×
  SqrtInto	2.43×	2.44×	2.47×

  MatVec 1024² is 2.02× faster and Dot 2^20 1.29×. On a riscv64 CPU without V (cfarm94, SiFive U74), the same build runs the scalar path and passes the same tests.

  GEMM micro-kernel (v0.7.1). NR is fixed by the packing, while VLEN varies between RISC-V cores. The V extension guarantees VLEN ≥ 128, so a group of four registers (LMUL=4) holds at least eight float64, and VSETIVLI $8 gives vl = 8 on every core: a 6×8 tile, six accumulator groups, the B row loaded once per step and each A value fused in with vfmacc.vf. Without V, the same 6×8 tile runs in Go. Measured with hostload -pick 1 to choose an idle CPU on a host otherwise 64% busy, and benchstat over 8 interleaved runs each against v0.7.0’s scalar 4×4 (all p = 0.000):

  one core, X60	64²	128²	256²	512²	1024²
  MatMul time	−31%	−32%	−45%	−43%	−44% (1.79 s → 1.00 s)

  Inner is −44%, and the small rectangular products are −22% to −38%. At about 2.1 GFLOP/s on one core this is still far from the hardware: VLEN = 256 with LMUL=4 uses half of each register group, which keeps one tile portable to every VLEN at that cost.

• s390x (vector facility), since v0.8.0, when AT_HWCAP reports it: sum, dot, sqrt and add/sub/mul/div, eight float64 per loop in four two-lane registers plus a scalar tail, every instruction a Go mnemonic. Since v0.9.0 the GEMM micro-kernel is a 4×8 tile: sixteen two-lane accumulators, the B row in one VLM, each A value replicated with VLREPG and fused in with VFMADB. Since v0.9.1, max/min are vector too on a z14 or later: VFMAXDB/VFMINDB in mode 1, Java’s Math.max/min, which gives NaN for a NaN operand and orders −0 below +0, NumPy’s rule. They come with the vector-enhancements facility 1 (AT_HWCAP VXRS_EXT), beyond the z13 that Go also supports, where the scalar reducers remain. On the z17, one core, benchstat over 8 interleaved runs, all at p = 0.000: Max takes 92% to 96% less time (1 Ki: 2.58 µs → 112 ns; 4 Mi: 10.7 ms → 0.81 ms), and MaxAxis(1) 94% less. Sum, untouched, measured within ±0.6%: the control. A sabotaged copy in mode 0 (IEEE maxNum, which drops a NaN) failed TestMaxMinSIMD on the z17. The operand order of VFSDB/VFDDB/VFMADB was measured, not read: each run once on known inputs (VFSDB V1, V2, V3 is V3 = V2 − V1). AT_HWCAP is read in the machine’s byte order, which on big-endian s390x is not the little-endian the parser first assumed.

  v0.8.0 shipped these kernels validated under emulation only. Since then they run on real hardware: an IBM z17 (Type 9175), through the LinuxONE Community Cloud, as a z/VM 7.3.0 guest with 8 vCPUs in a shared LPAR. There the whole suite passes, as do the fence tests (a kernel reading one element too far faults), TestParallelStress (27,752 calls) and two sabotaged GEMM and elementwise kernels, which fail as they should. Measured with hostload (0% busy, 0% steal) and benchstat, 8 interleaved runs each, every difference at p = 0.000 unless marked:

  IBM z17, time vs the scalar code	1 core	8 cores
  Sum	−85% to −87%	−82% to −87%
  AddInto / MulInto	−44% to −61%	−36% to −60%
  DivInto	−6% to −18%	−6% to −22%
  SqrtInto	−2% to −11%	−5% to −9% at 1–16 Ki, ~ above
  MatVec 1024²	−44%	−40%
  Dot 2^20	−11%	~ (p = 0.08)
  MatMul 1024² (v0.9.0 GEMM)	340 → 70 ms (−79%)	43 → 10 ms (−76%)
  GEMM geomean, 11 shapes	−74%	−66%

  Exp, which none of this touches, measured within ±0.3% on one core: the control.

• The 32-bit targets keep the validated scalar oracles, using the same four-accumulator max/min, direct sqrt loop, and a scalar 4×4 GEMM micro-kernel over the packed panels, and still get the packing + cache blocking + multicore structure (they have not been measured against NumPy). What the Go assembler offers them, checked on Go 1.26.4 and 1.27.1 by assembling and disassembling (2026-10-04): ppc64le has no vector-double arithmetic (no XVADDDP/XVMADDADP), which go-asmgen v0.13.0 now encodes as WORDs, so ppc64le has kernels; loong64 has vector-double add/sub/mul/div (VADDD assembles to vfadd.d, VMULD to vfmul.d, XVADDD to xvfadd.d) and the broadcast load (XVMOVQ off(R), X.V4 = xvldrepl.d), but no vector FMA, which go-asmgen v0.14.0 encodes (transitionally), so loong64 has kernels; s390x has them, FMA included (VFADB, VFMADB); riscv64 has them (VFADDVV, VFMACCVV) and has had kernels since v0.7.0, s390x since v0.8.0. (s390x additionally exercises the big-endian path in CI.)

All six are exercised in CI (native amd64/arm64 + qemu for the rest); each per-arch job regenerates the committed .s, fails if it is stale, vets (asmdecl), builds (cmd/asm encodes), and runs the bit/NaN-correctness suite. The multicore path and the packed/cache-blocked GEMM driver are architecture-independent; only the GEMM micro-kernel is per-arch (NEON 4×8 on arm64, AVX2/FMA 6×8 or SSE2 on amd64, scalar 4×4 on the other four), each held to the scalar ikj oracle in CI: exactly on integer-valued data, within tolerance on general data, since above KC the block partials regroup the sum and the FMA kernels round once per step where the scalar oracle may round twice.