diff --git a/README.md b/README.md index 307bf4f..b66ecc1 100644 --- a/README.md +++ b/README.md @@ -57,33 +57,56 @@ from 1.18 ms to 0.39 ms (7.5× NumPy). ⚠ An array allocated from the workspace is invalid after `Reset`; `Detach` is how a result leaves it. The `*Into` methods (`AddInto`, `SqrtInto`, …) remain the way to reuse one named buffer. -The numeric inner loops are kept behind a narrow kernel API. Behind it, -**large** elementwise ops, reductions and matmul run **multicore** (across -`GOMAXPROCS`); the sum reduction uses a **go-asmgen SIMD kernel** (NEON on -arm64, SSE2 on amd64); and `MatMul` is a **panel-packed, cache-blocked GEMM** -with a go-asmgen SIMD-FMA micro-kernel (NEON 4×8 on arm64; AVX2/FMA 6×8 on -amd64, chosen at run time, with an SSE2 fallback) — -the OpenBLAS/BLIS structure. So go-ndarray **beats single-threaded NumPy** on -its core ops on large arrays (Add/Mul ~2×, Sum ~2.4×), and `MatMul` reaches -**tuned-BLAS parity at 1024²**: ≈0.99× vs multi-threaded OpenBLAS 0.3.29 on an -arm64 VM (~203 GFLOP/s, up from an earlier ~156 GFLOP/s / 0.76× before the -by-lane-FMLA micro-kernel fix — see **[docs/perf.md](docs/perf.md)**), and -≈1.00× of single-threaded vecLib on an Apple M4 Max (~373 GFLOP/s), where it -also decisively **beats the pure-Go `gonum` 4–10×** at every size — see -**[BENCHMARKS.md](BENCHMARKS.md)**. Where a faster reference exists — -`Sqrt`/`Max` small sizes, and tuned BLAS at small matmul sizes — it says so. -It is a **standalone, reusable** module and the cgo-free ndarray backend behind +The numeric inner loops are kept behind a narrow kernel API. Behind it, large +elementwise ops, reductions and products run **multicore** (across +`GOMAXPROCS`). On **amd64 and arm64** the hot loops are +[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels: sum, +max/min, add/sub/mul/div, sqrt and the dot product (SSE2 or AVX2/FMA on amd64, +NEON on arm64), and a **panel-packed, cache-blocked GEMM** with an SIMD-FMA +micro-kernel (NEON 4×8; AVX2/FMA 6×8 chosen at run time, SSE2 fallback), the +OpenBLAS/BLIS structure. The other four 64-bit targets (riscv64, loong64, +ppc64le, s390x) and the 32-bit ones run the same pure-Go code those kernels +are tested against. `Exp` and `Log` are ports of Arm's optimized-routines +(≤ 0.51 ULP measured, against up to 0.88 and 0.72 for Go's +`math.Exp`/`math.Log` on arm64), +and they are correct where amd64's `math.Exp` returns +Inf +([golang/go#81995](https://github.com/golang/go/issues/81995)) and `math.Log` +is wrong on subnormals +([golang/go#56600](https://github.com/golang/go/issues/56600)). + +**Measured against NumPy 2.5.3 (OpenBLAS 0.3.34) on an AMD Zen 3, 16 cores** +(full tables in **[docs/perf.md](docs/perf.md)**): + +- **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×), the `*Into` forms of `Add`/`Mul`/`Sqrt` (faster or at + parity at every size), 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 (0.55× at 1024², 0.3× at + 256² and below); `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×: Go's allocator and garbage collector, which + `Workspace` largely removes). + +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 every product beats the +pure-Go `gonum` 4–10× (**[BENCHMARKS.md](BENCHMARKS.md)**). It is a +**standalone, reusable** module and the cgo-free ndarray backend behind [go-embedded-ruby](https://github.com/go-embedded-ruby/ruby)'s `NDArray` class. -> ⚠️ **Status: float64 NumPy parity for the core surface, and faster than NumPy -> on it.** Creation, slicing/views, broadcasting elementwise + ufuncs, -> reductions (incl. arg/cumulative/clip/where), manipulation, and linear algebra -> are complete, **100%-covered**, six-arch CI-green, and differentially checked -> against numpy. The hot paths are **multicore + SIMD** and beat single-threaded -> NumPy on large arrays (**[docs/perf.md](docs/perf.md)**). See -> **[docs/plan-ndarray.md](docs/plan-ndarray.md)** for the roadmap (more SIMD -> kernels, more dtypes). The Ruby binding has shipped — go-embedded-ruby's -> `NDArray` class binds this module. +> ⚠️ **Status: float64 NumPy parity for the core surface.** Creation, +> slicing/views, broadcasting elementwise + ufuncs, reductions (incl. +> arg/cumulative/clip/where), manipulation, linear algebra and `Workspace` are +> complete, **100%-covered**, and differentially checked against NumPy. CI runs +> the suite on amd64, arm64, 386, and riscv64/loong64/ppc64le/s390x/arm under +> qemu, on Linux, macOS and Windows, and compiles it for every `GOOS/GOARCH` +> pair Go supports; v0.1.0 and v0.2.0 were also run on real amd64, arm64, +> ppc64le, riscv64 and loong64 hardware. See **[docs/plan-ndarray.md](docs/plan-ndarray.md)** +> for the roadmap (more dtypes, more SIMD targets). ## Why this module? @@ -91,10 +114,10 @@ It is a **standalone, reusable** module and the cgo-free ndarray backend behind **amd64-only**. More broadly, **Ruby has no cgo-free ndarray** (`Numo::NArray`, `NMatrix` are C extensions). A pure-Go core whose kernels are generated for every arch is therefore a durable foundation. The numeric loops live in -`internal/kernels`; **Phase 1** replaces them with -[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels across all six -64-bit Go targets (amd64, arm64, riscv64, loong64, ppc64le, s390x), selected at -runtime, behind the same API and the same tests. +`internal/kernels`, each with a pure-Go reference implementation; +[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels replace them +on amd64 and arm64 today, behind the same API and held to the same tests, and +the other 64-bit targets are next. ## Example diff --git a/docs/perf.md b/docs/perf.md index 07039e2..e4f957f 100644 --- a/docs/perf.md +++ b/docs/perf.md @@ -340,7 +340,7 @@ 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, median of two runs: +Zen 3, 16 cores pinned, NumPy 2.5.3 single-threaded, best of two runs: | op | n | heap | **Workspace** | `*Into` | NumPy | Workspace vs NumPy | |----|--:|--:|--:|--:|--:|:--:| @@ -381,14 +381,17 @@ 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, the FMA ports' peak at that clock; the - SSE2 4×4 it replaces had no FMA and half the vector width. + 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 a 512² - product (profiled). A partial tile now runs the kernel into a scratch tile. + 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, stalled on memory. A sweep of MC (48–252) and KC + 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% @@ -408,9 +411,11 @@ below 256²). Where the parallel loss is, measured on the Zen 3: 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, which stalls on the shared packed B panel. On this - part the last-level cache is split between core complexes, so the panel is - partly remote. + 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 + 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 @@ -420,9 +425,9 @@ 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 lever is therefore a shared B panel per core complex (or per cache -domain), not a different split of C, then an AVX-512 kernel for the hosts that -have it. +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. 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) @@ -482,18 +487,23 @@ 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. `math.Exp` measured 0.879 ULP on - the same inputs on arm64. The table was not trusted as copied: + 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:** on amd64, Go 1.26's `math.Exp` returns **+Inf for - x ≥ 1023.5·ln 2 ≈ 709.436**, although the result is finite up to - ln(MaxFloat64) ≈ 709.7827 (NumPy: `exp(709.5)` = 1.3549863193146328e+308). - Its assembly rounds k = x·log2(e) to 1024 there, and the biased-exponent - check treats that as overflow although the reduced factor is below 1. Seen on - real Zen 3 hardware, not only Rosetta. go-ndarray's `Exp` inherited it until - now; `TestExpTopOfRange` pins the fix. -- Also just above ln(2^-1075), exp now rounds up to the smallest subnormal as - it should, where `math.Exp` returns 0. +- **A Go bug it removes**, reported as + [golang/go#81995](https://github.com/golang/go/issues/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, @@ -517,9 +527,12 @@ held to the n·ε·Σ|aᵢbᵢ| bound and exact on integer data. | Dot 2^20, 16 cores | 159 µs | 99 µs | — | | MatVec 1024², 16 cores | 165 µs | 89 µs | — | -NumPy's multi-threaded dot still wins on this part (24 µs): OpenMP pins its -threads, so each one rereads its own chunk from its own L3 slice on every -repetition, while goroutines move between core complexes. +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. + ## Log and Log10 (2026-10-04) `Log` went through `Map(math.Log)`: 11.6 ns per element on a Zen 3 core @@ -532,17 +545,20 @@ inputs within 2⁻⁴ of 1 take a degree-12 polynomial. z/c − 1 is one needs is not used. `Log10` is log(x)·(1/ln10), the formula of `math.Log10`, on this log. -- **Accuracy:** worst 0.508 ULP over 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.727 on arm64. `TestLogTable` checks each of the 128 +- **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:** on amd64, Go 1.26's `math.Log` is wrong for - subnormal inputs. `math.Log(5e-324)` returns −709.09 instead of −744.44, and - `math.Log(1e-310)` −709.09 instead of −713.80 (NumPy: −744.44, −713.80); the - error measured 7·10¹³ ULP. Its assembly does not normalise them. Confirmed on - a real Zen 3. `Log` and `Log10` inherited it; `TestLogSubnormal` pins the fix. +- **A second Go bug it removes**, known since 2022 as + [golang/go#56600](https://github.com/golang/go/issues/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 @@ -572,12 +588,18 @@ Speed, 16 Ki elements, one core: M4 37 µs (`math.Log` 72), Zen 3 73 µs plain vector FP add exists). All beat/parity NumPy via SIMD + multicore, no `func`-pointer indirection. - The other four 64-bit Go targets — **riscv64, loong64, ppc64le, s390x** — keep - the validated scalar oracles (Go's loong64/ppc64le assemblers expose no - vector-double arithmetic; riscv64's V extension is optional), 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** win, so they also beat single-threaded NumPy on large - arrays. (s390x additionally exercises the big-endian path in CI.) + 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`); + **loong64** has vector-double add/sub/mul/div (`VADDD` assembles to + `vfadd.d`, `VMULD` to `vfmul.d`, `XVADDD` to `xvfadd.d`) but no vector FMA; + **s390x** has them, FMA included (`VFADB`, `VFMADB`); **riscv64** has them + (`VFADDVV`, `VFMACCVV`), but the V extension is optional and needs a run-time + check. So loong64, s390x and riscv64 kernels are work not yet done, not a + toolchain wall. (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 diff --git a/internal/kernels/exp_test.go b/internal/kernels/exp_test.go index e48901d..3c06be7 100644 --- a/internal/kernels/exp_test.go +++ b/internal/kernels/exp_test.go @@ -144,8 +144,8 @@ func TestExpSpecial(t *testing.T) { } } // Against the reference: the subnormal range (just above ln(2^-1075) exp - // rounds up to the smallest subnormal; math.Exp gives 0), and the top of - // the finite range, where Go 1.26's amd64 math.Exp returns +Inf from + // rounds up to the smallest subnormal; arm64's math.Exp gives 0), and the + // top of the finite range, where amd64's math.Exp (golang/go#81995) returns +Inf from // x >= 1023.5*ln2 ~ 709.436 (its k = round(x*log2e) reaches 1024 and the // biased exponent check calls that overflow, though fr*2^1024 is finite). for _, x := range []float64{-745.13, -745.1332, -740, -720.5, -708.5, -708.39, diff --git a/internal/kernels/gemm_generic.go b/internal/kernels/gemm_generic.go index cd3a82f..ba202f5 100644 --- a/internal/kernels/gemm_generic.go +++ b/internal/kernels/gemm_generic.go @@ -4,9 +4,9 @@ package kernels // Generic GEMM micro-kernel for the four targets without a vector-double SIMD // kernel (loong64, ppc64le, riscv64, s390x): a portable scalar 4x4 register tile -// over the packed panels. There is no hand-vectorized .s here — Go's loong64/ -// ppc64le assemblers expose no vector-double arithmetic and riscv64's V extension -// is optional — but the kernel still gets the two structural wins of the packed +// over the packed panels. There is no hand-vectorized .s here yet (ppc64le's Go +// assembler has no vector-double arithmetic; loong64, s390x and riscv64 have it, +// see docs/perf.md) — but the kernel still gets the two structural wins of the packed // GEMM: the panels are contiguous unit-stride (so the loads are conflict-free, // the L1-collision fix), and the cache-blocking loop nest plus the multicore // fan-out (MatMulP) keep it well ahead of an unblocked scalar loop and beat diff --git a/internal/kernels/reduce_generic.go b/internal/kernels/reduce_generic.go index 8ed8ee7..6a3d38a 100644 --- a/internal/kernels/reduce_generic.go +++ b/internal/kernels/reduce_generic.go @@ -10,9 +10,10 @@ package kernels // amd64 (SSE2) and arm64 (NEON sum + intrinsic scalar sqrt/max/min) ship kernels // (see reduce_amd64.go / reduce_arm64.go), so this file is the // loong64/ppc64le/riscv64/s390x fallback. For those four no hand-vectorized -// kernel is shipped: Go's loong64/ppc64le assemblers expose no vector-double -// arithmetic (the same wall go-fft documented), riscv64's V extension is -// optional, and the per-arch qemu jobs still exercise this dispatch. The scalar +// kernel is shipped yet. Go's ppc64le assembler has no vector-double +// arithmetic; loong64 has vector add/mul (VADDD, VMULD) but no vector FMA, +// s390x has both, and riscv64's V extension is optional (a run-time check); +// see docs/perf.md. The per-arch qemu jobs still exercise this dispatch. The scalar // math.Sqrt/math.Max/math.Min still lower to the hardware sqrt/max where the // target has it, and the multicore path beats single-threaded numpy on large // arrays regardless. diff --git a/ufunc.go b/ufunc.go index daf6de6..f04e935 100644 --- a/ufunc.go +++ b/ufunc.go @@ -50,8 +50,8 @@ func (a *Array) SqrtInto(out *Array) error { // // It runs a table-driven port of Arm's optimized-routines exp (0.51 ULP worst // case, measured), split across cores: several times faster than math.Exp per -// element, and finite up to ln(MaxFloat64) on amd64, where Go 1.26's math.Exp -// returns +Inf from x ~ 709.436. +// element, and finite up to ln(MaxFloat64) on amd64, where math.Exp returns +// +Inf from x ~ 709.436 (golang/go#81995; Go 1.26.8 and 1.27.1). func (a *Array) Exp() *Array { src := a.contiguousData() dst := a.alloc(len(src), false) @@ -64,8 +64,8 @@ func (a *Array) Exp() *Array { // // It runs a table-driven port of Arm's optimized-routines log (0.51 ULP worst // case, measured), split across cores: about twice as fast as math.Log, and -// correct for subnormal inputs on amd64, where Go 1.26's math.Log is not -// (log(5e-324) gives -709.09 instead of -744.44). +// correct for subnormal inputs on amd64, where math.Log is not (log(5e-324) +// gives -709.09 instead of -744.44; golang/go#56600). func (a *Array) Log() *Array { return a.unary(kernels.LogP) } // Log2 returns the elementwise base-2 logarithm, within 2 ULP of numpy's. diff --git a/ufunc_test.go b/ufunc_test.go index ca5263b..22ff359 100644 --- a/ufunc_test.go +++ b/ufunc_test.go @@ -104,8 +104,8 @@ func TestComparisonBroadcastAndError(t *testing.T) { } } -// TestExpTopOfRange: exp is finite up to ln(MaxFloat64) on every arch. Go -// 1.26's amd64 math.Exp returns +Inf from x ~ 709.436, which Exp used to +// TestExpTopOfRange: exp is finite up to ln(MaxFloat64) on every arch. +// amd64's math.Exp returns +Inf from x ~ 709.436 (golang/go#81995), which Exp used to // inherit; numpy gives 1.3549863193146328e+308 for 709.5. The input is a // strided view, so the contiguous copy path is exercised too. func TestExpTopOfRange(t *testing.T) { @@ -118,7 +118,7 @@ func TestExpTopOfRange(t *testing.T) { } // TestLogSubnormal: Log and Log10 are right for subnormal inputs on every -// arch. Go 1.26's amd64 math.Log returns -709.09 for both 5e-324 and 1e-310, +// arch. amd64's math.Log returns -709.09 for both 5e-324 and 1e-310 (golang/go#56600), // which Log and Log10 used to inherit; numpy gives -744.44 and -713.80. func TestLogSubnormal(t *testing.T) { x := mustArr(t, ok(FromData([]float64{5e-324, 1e-310, 1, math.E}, 4)))