perf: VSX kernels on ppc64le (sum, dot, sqrt, elementwise, 8x8 GEMM)

.github/workflows/ci.yml

@@ -14,8 +14,9 @@ permissions:
 # A single `go test -coverpkg` run cannot reach 100% statement coverage over a
 # SIMD-dispatched package: each arch compiles exactly one sumSIMD wrapper — the
 # amd64 SSE2 wrapper (reduce_amd64.go -> sum_amd64.s), the arm64 NEON wrapper
-# (reduce_arm64.go -> sum_arm64.s), or the scalar-alias wrapper
-# (reduce_generic.go, on loong64/ppc64le/riscv64/s390x). No single build sees
+# (reduce_arm64.go -> sum_arm64.s), the ppc64le VSX wrapper (reduce_ppc64le.go
+# -> sum_ppc64le.s), or the scalar-alias wrapper (reduce_generic.go, on
+# loong64/riscv64/s390x and 32-bit). No single build sees
 # more than one, so no single run covers them all.
 #
 # The honest split:
@@ -36,10 +37,9 @@ permissions:
 #
 # So every line is either statement-covered by the gate or executed by a
 # per-arch job; the gate is never lowered and there is no coverage-gaming knob.
-# amd64/arm64 ship a SIMD sum kernel; the other four keep the validated scalar
-# path (Go's loong64/ppc64le assemblers expose no vector double arithmetic, and
-# riscv64's V extension is optional) and the multicore path still beats
-# single-threaded numpy on large arrays.
+# amd64/arm64/ppc64le ship SIMD kernels (ppc64le's run under qemu here and were
+# run on real POWER8/POWER9); loong64/riscv64/s390x keep the validated scalar
+# path for now.
 jobs:
   test:
     name: test (${{ matrix.os }})
@@ -123,6 +123,13 @@ jobs:
         with:
           go-version: stable
       - run: sudo apt-get update && sudo apt-get install -y qemu-user-static
+      - name: Regenerate committed assembly and check it is current
+        if: matrix.goarch == 'ppc64le'
+        run: |
+          ( cd internal/kernels/asmgen/ppc64le && go run . )
+          mv internal/kernels/asmgen/ppc64le/sum_ppc64le.s internal/kernels/sum_ppc64le.s
+          git diff --exit-code -- internal/kernels/sum_ppc64le.s || {
+            echo "::error::generated ppc64le assembly is stale; run go generate and commit"; exit 1; }
       - name: Vet + test under qemu (${{ matrix.goarch }})
         run: |
           GOOS=linux GOARCH=${{ matrix.goarch }} go vet ./...

README.md

@@ -59,14 +59,14 @@ 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 products run **multicore** (across
-`GOMAXPROCS`). On **amd64 and arm64** the hot loops are
+`GOMAXPROCS`). On **amd64, arm64 and ppc64le** 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
+add/sub/mul/div, sqrt and the dot product (SSE2 or AVX2/FMA on amd64, NEON on
+arm64, VSX on ppc64le), max/min too on amd64, 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; VSX 8×8), the OpenBLAS/BLIS structure. The
+other 64-bit targets (riscv64, loong64, 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

docs/perf.md

@@ -565,6 +565,44 @@ this log.
 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 `WORD`s, 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):
+
+| op | before | after | gain |
+|----|--:|--:|:--:|
+| micro-kernel, 1 core | 8.3 GFLOP/s | 16.2 GFLOP/s | ×1.96 |
+| serial GEMM 512², 1 core | 7.9 GFLOP/s | 15.8 GFLOP/s | ×2.0 |
+| MatMul 1024², 32 threads | 39.2 ms | 23.6 ms | ×1.66 |
+| MatMul 128² / 512² | | | ×1.8 / ×1.5 |
+| MatVec 1024² | 823 µs | 267 µs | ×3.1 |
+| Sum 1 Ki / 16 Ki / 256 Ki | | | ×6.6 / ×7.3 / ×3.2 |
+| AddInto 1 Ki–256 Ki | | | ×2.7–4.1 |
+| Sqrt, 1 core | | | ×1.05–1.18 |
+
+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.
+
 ## SIMD coverage
 
 - **amd64 (SSE2)** ships hand-vectorized `sum` (4-accumulator `ADDPD`), `sqrt`
@@ -588,13 +626,17 @@ Speed, 16 Ki elements, one core: M4 37 µs (`math.Log` 72), Zen 3 73 µs
   `FMAXD`/`FMIND`); the C-accumulate fold still uses `VFMLA` against `1.0` (no
   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
+- **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.
+- The other three 64-bit Go targets — **riscv64, loong64, s390x** — 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`);
+  **ppc64le** has no vector-double arithmetic (no `XVADDDP`/`XVMADDADP`), which
+  go-asmgen v0.13.0 now encodes as `WORD`s, 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`) but no vector FMA;
   **s390x** has them, FMA included (`VFADB`, `VFMADB`); **riscv64** has them

internal/kernels/asmgen/go.mod

@@ -6,4 +6,4 @@ module github.com/go-ndarray/ndarray/internal/kernels/asmgen
 
 go 1.26.4
 
-require github.com/go-asmgen/asmgen v0.11.0
+require github.com/go-asmgen/asmgen v0.13.0

internal/kernels/asmgen/go.sum

@@ -1,2 +1,2 @@
-github.com/go-asmgen/asmgen v0.11.0 h1:BQKTBAPEZ1LDS/gdvSemHxUzUGKx9pfHNQtbcln8q8g=
-github.com/go-asmgen/asmgen v0.11.0/go.mod h1:Hj29fXioLRe9oSnQ18SuDBiOcCbu8Qwxba6/5TMr4+Y=
+github.com/go-asmgen/asmgen v0.13.0 h1:O3DIZtgTgwHJ2uof6NhMCO1AnTZqub24r3DS83IjOb8=
+github.com/go-asmgen/asmgen v0.13.0/go.mod h1:Hj29fXioLRe9oSnQ18SuDBiOcCbu8Qwxba6/5TMr4+Y=