simd(s390x): vector-facility kernels; native-endian AT_HWCAP

.github/workflows/ci.yml

@@ -17,7 +17,8 @@ permissions:
 # (reduce_arm64.go -> sum_arm64.s), the ppc64le VSX wrapper (reduce_ppc64le.go
 # -> sum_ppc64le.s), the loong64 LASX wrapper (reduce_loong64.go ->
 # sum_loong64.s), the riscv64 RVV wrapper (reduce_riscv64.go -> sum_riscv64.s),
-# or the scalar-alias wrapper (reduce_generic.go, on s390x and 32-bit). No
+# the s390x vector wrapper (reduce_s390x.go -> sum_s390x.s), or the
+# scalar-alias wrapper (reduce_generic.go, on the 32-bit targets). No
 # single build sees
 # more than one, so no single run covers them all.
 #
@@ -42,8 +43,8 @@ permissions:
 # amd64/arm64/ppc64le/loong64 ship SIMD kernels (ppc64le's and loong64's were
 # run on real POWER8/POWER9 and a Loongson 3C5000L; under qemu here loong64
 # takes its scalar fallback unless the emulated CPU reports LASX);
-# riscv64 ships RVV kernels (used when the CPU reports V); s390x keeps the
-# validated scalar path for now.
+# riscv64 ships RVV kernels (used when the CPU reports V) and s390x
+# vector-facility kernels (used when AT_HWCAP reports it).
 jobs:
   test:
     name: test (${{ matrix.os }})
@@ -128,7 +129,7 @@ jobs:
           go-version: 1.27.1
       - 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' || matrix.goarch == 'loong64' || matrix.goarch == 'riscv64'
+        if: matrix.goarch == 'ppc64le' || matrix.goarch == 'loong64' || matrix.goarch == 'riscv64' || matrix.goarch == 's390x'
         run: |
           ( cd internal/kernels/asmgen/${{ matrix.goarch }} && go run . )
           mv internal/kernels/asmgen/${{ matrix.goarch }}/sum_${{ matrix.goarch }}.s internal/kernels/sum_${{ matrix.goarch }}.s

README.md

@@ -59,15 +59,16 @@ 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, arm64, ppc64le, loong64 and riscv64** the hot loops
-are [go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels: sum,
+`GOMAXPROCS`). On **amd64, arm64, ppc64le, loong64, riscv64 and s390x** the hot
+loops are [go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels: sum,
 add/sub/mul/div, sqrt and the dot product (SSE2 or AVX2/FMA on amd64, NEON on
-arm64, VSX on ppc64le, LASX on loong64 and RVV on riscv64 when the CPU has
-them), max/min too on
+arm64, VSX on ppc64le, LASX on loong64, RVV on riscv64 and the vector facility
+on s390x, the last three when the CPU has them), max/min too on
 amd64 and 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; VSX
-8×8; LASX 8×8; RVV 6×8), the OpenBLAS/BLIS structure. s390x and the 32-bit targets run
-the same pure-Go code those kernels are tested against. `Exp` and `Log` are ports of Arm's optimized-routines
+8×8; LASX 8×8; RVV 6×8), the OpenBLAS/BLIS structure. The 32-bit targets run
+the same pure-Go code those kernels are tested against. The s390x kernels are
+validated under emulation only (no s390x host was reachable to measure them). `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
@@ -119,8 +120,8 @@ pure-Go `gonum` 4–10× (**[BENCHMARKS.md](BENCHMARKS.md)**). It is a
 arch is therefore a durable foundation. The numeric loops live in
 `internal/kernels`, each with a pure-Go reference implementation;
 [go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels replace them
-on amd64, arm64, ppc64le, loong64 and riscv64 today, behind the same API and
-held to the same tests; s390x is next.
+on every 64-bit target Go supports (amd64, arm64, ppc64le, loong64, riscv64,
+s390x), behind the same API and held to the same tests.
 
 ## Goroutines
 

SECURITY.md

@@ -10,7 +10,7 @@ issue for it.
 
 go-ndarray is pure Go (`CGO_ENABLED=0`) except for its SIMD kernels, which are
 Plan 9 assembly generated by [go-asmgen](https://github.com/go-asmgen/asmgen)
-for amd64, arm64, ppc64le, loong64 and riscv64. Those kernels are the only code Go's
+for amd64, arm64, ppc64le, loong64, riscv64 and s390x. Those kernels are the only code Go's
 bounds checks do not cover, so they are held to the following:
 
 - **Every wrapper checks the reach of its kernel.** Before an assembly call, the
@@ -57,6 +57,16 @@ standard library's known vulnerabilities; `govulncheck ./...` reports none.
 
 ## Audit log
 
+**2026-10-08 (v0.8.0).** New s390x vector-facility kernels, held to the same
+guards (wrappers check each kernel's reach; the PROT_NONE fence tests now
+include s390x). ⚠ They were validated **under emulation only** (Docker's
+qemu-s390x locally, qemu-s390x in CI): the whole suite and the fence tests
+pass, and two sabotaged kernels (a subtraction with swapped operands, a dot
+that does not advance `b`) both failed. No s390x host was reachable. Found on
+the way: `hwcapFrom` read `/proc/self/auxv` as little-endian; on big-endian
+s390x that would have read the wrong word, so it now uses the machine's byte
+order, and `useVX` comes out true under emulation as it should.
+
 **2026-10-08 (v0.7.0).** New riscv64 RVV kernels. They are held to the same
 guards: each wrapper checks the reach of its kernel (`TestWrappersRefuseShortOperands`),
 and the PROT_NONE fence tests pass on a real SpacemiT X60 with V (cfarm95), as

docs/perf.md

@@ -875,7 +875,19 @@ stay dynamic, so its slower efficiency cores are handled as before.
   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.
-- The remaining 64-bit Go target — **s390x** — keeps
+- **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.
+  max/min and the GEMM micro-kernel stay scalar (the vector max with NumPy's
+  NaN rule needs the z14, beyond the z13 Go supports). 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). ⚠ **Validated under emulation
+  only** (Docker's qemu-s390x locally, qemu-s390x in CI): the whole suite, the
+  fence tests, and two sabotaged kernels that both failed. **Not measured**:
+  the LinuxONE host did not answer. AT_HWCAP is read in the machine's byte
+  order, which on big-endian s390x is not the little-endian the parser first
+  assumed.
+- 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
@@ -888,8 +900,8 @@ stay dynamic, so its slower efficiency cores are handled as before.
   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. So s390x kernels are
-  work not yet done, not a toolchain wall. (s390x additionally exercises the big-endian path in CI.)
+  (`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

docs/plan-ndarray.md

@@ -126,7 +126,9 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
     VSX encoders, and since v0.2.5 loong64 has LASX kernels (used when the CPU
     reports LASX), and since v0.7.0 riscv64 has RVV kernels (sum, dot, sqrt,
     elementwise, and since v0.7.1 a 6×8 GEMM micro-kernel; used when the CPU
-    reports V). The other targets (s390x, 32-bit) keep the
+    reports V), and since v0.8.0 s390x has vector-facility kernels (sum, dot,
+    sqrt, elementwise; validated under emulation, not yet measured). The
+    32-bit targets keep the
     validated scalar oracles plus a scalar 4×4 GEMM micro-kernel over the
     packed panels, with packing + cache-blocking + multicore; they have not
     been measured against NumPy. Split-CI: the pure-Go

internal/kernels/asmgen/s390x/gen.go

@@ -0,0 +1,140 @@
+// Command gen produces sum_s390x.s, the vector-facility float64 kernels for
+// s390x, via go-asmgen. Run with: go run . (or `go generate` from the kernels
+// package).
+//
+// Go's s390x assembler names the instructions (VL/VLM/VSTM, VFADB, VFSDB,
+// VFMDB, VFDDB, VFSQDB, VFMADB, VREPG, WFADB), so nothing is a WORD. A vector
+// register holds two float64; each loop does eight (four registers) and a
+// scalar tail does the rest. The Go side checks the vector facility in
+// AT_HWCAP.
+//
+// Operand order, measured by running each instruction on known inputs under
+// s390x emulation (the IBM manuals give the other order): VFSDB V1, V2, V3 is
+// V3 = V2 - V1, VFDDB V1, V2, V3 is V3 = V2 / V1, VFMADB V1, V2, V3, V4 is
+// V4 = V1*V2 + V3, and FSUB F1, F0 is F0 = F0 - F1.
+//
+// R1..R4 hold arguments and pointers (R0 cannot be a base register); V16..V19
+// are the accumulators; F0..F15 alias the first element of V0..V15.
+package main
+
+import (
+	"fmt"
+	"os"
+
+	"github.com/go-asmgen/asmgen/emit"
+	"github.com/go-asmgen/asmgen/s390x"
+)
+
+func main() {
+	f := emit.NewFile("s390x")
+	f.Add(reduceKernel("sumVX", false))
+	f.Add(reduceKernel("dotVX", true))
+	f.Add(sqrtKernel())
+	for _, k := range []struct{ name, vec, sca string }{
+		{"addVX", "VFADB", "FADD"}, {"subVX", "VFSDB", "FSUB"},
+		{"mulVX", "VFMDB", "FMUL"}, {"divVX", "VFDDB", "FDIV"},
+	} {
+		f.Add(binKernel(k.name, k.vec, k.sca))
+	}
+	if err := os.WriteFile("sum_s390x.s", []byte(f.String()), 0o644); err != nil {
+		fmt.Fprintln(os.Stderr, err)
+		os.Exit(1)
+	}
+	fmt.Println("wrote sum_s390x.s")
+}
+
+// reduceKernel builds NAME(a[, b] *float64, n int) float64: the sum of a, or
+// with dot the sum of a[i]*b[i]. Four two-lane accumulators, folded at the end
+// into F0, then the tail added one element at a time.
+func reduceKernel(name string, dot bool) *emit.Function {
+	args, types := []string{"a", "n"}, []s390x.Type{s390x.Ptr, s390x.Int64}
+	if dot {
+		args, types = []string{"a", "b", "n"}, []s390x.Type{s390x.Ptr, s390x.Ptr, s390x.Int64}
+	}
+	sig := s390x.Layout(args, types, []string{"ret"}, []s390x.Type{s390x.Float64})
+	b := s390x.NewFunc(name, sig, 0)
+	b.LoadArg("a", "R1")
+	if dot {
+		b.LoadArg("b", "R3")
+	}
+	b.LoadArg("n", "R2")
+	for r := 16; r < 20; r++ {
+		b.Raw("VZERO V%d", r)
+	}
+	b.Raw("loop:").Raw("CMP R2, $8").Raw("BLT fold").
+		Raw("VLM (R1), V0, V3")
+	if dot {
+		b.Raw("VLM (R3), V4, V7")
+		for r := 0; r < 4; r++ {
+			b.Raw("VFMADB V%d, V%d, V%d, V%d", r, 4+r, 16+r, 16+r) // acc += a*b
+		}
+		b.Raw("ADD $64, R3")
+	} else {
+		for r := 0; r < 4; r++ {
+			b.Raw("VFADB V%d, V%d, V%d", r, 16+r, 16+r) // acc += a
+		}
+	}
+	b.Raw("ADD $64, R1").Raw("SUB $8, R2").Raw("BR loop").
+		Raw("fold:").
+		Raw("VFADB V17, V16, V16").
+		Raw("VFADB V19, V18, V18").
+		Raw("VFADB V18, V16, V0").
+		Raw("VREPG $1, V0, V1").
+		Raw("WFADB V0, V1, V0"). // F0 = lane 0 + lane 1
+		Raw("tail:").Raw("CMP R2, $0").Raw("BEQ done").
+		Raw("FMOVD (R1), F1")
+	if dot {
+		b.Raw("FMOVD (R3), F2").Raw("FMUL F2, F1").Raw("ADD $8, R3")
+	}
+	b.Raw("FADD F1, F0").
+		Raw("ADD $8, R1").Raw("SUB $1, R2").Raw("BR tail").
+		Raw("done:").
+		StoreRet("F0", "ret").
+		Ret()
+	return b.Func()
+}
+
+// sqrtKernel builds sqrtVX(dst, src *float64, n int).
+func sqrtKernel() *emit.Function {
+	sig := s390x.Layout([]string{"dst", "src", "n"},
+		[]s390x.Type{s390x.Ptr, s390x.Ptr, s390x.Int64}, nil, nil)
+	b := s390x.NewFunc("sqrtVX", sig, 0)
+	b.LoadArg("dst", "R1").LoadArg("src", "R2").LoadArg("n", "R4")
+	b.Raw("loop:").Raw("CMP R4, $8").Raw("BLT tail").
+		Raw("VLM (R2), V0, V3")
+	for r := 0; r < 4; r++ {
+		b.Raw("VFSQDB V%d, V%d", r, r)
+	}
+	b.Raw("VSTM V0, V3, (R1)").
+		Raw("ADD $64, R1").Raw("ADD $64, R2").Raw("SUB $8, R4").Raw("BR loop").
+		Raw("tail:").Raw("CMP R4, $0").Raw("BEQ done").
+		Raw("FMOVD (R2), F0").Raw("FSQRT F0, F0").Raw("FMOVD F0, (R1)").
+		Raw("ADD $8, R1").Raw("ADD $8, R2").Raw("SUB $1, R4").Raw("BR tail").
+		Raw("done:").Ret()
+	return b.Func()
+}
+
+// binKernel builds NAME(dst, a, b *float64, n int): dst = a OP b. vec computes
+// Vd = Vn OP Vm from "vec Vm, Vn, Vd"; sca computes F0 = F0 OP F1 from
+// "sca F1, F0".
+func binKernel(name, vec, sca string) *emit.Function {
+	sig := s390x.Layout([]string{"dst", "a", "b", "n"},
+		[]s390x.Type{s390x.Ptr, s390x.Ptr, s390x.Ptr, s390x.Int64}, nil, nil)
+	b := s390x.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "R1").LoadArg("a", "R2").LoadArg("b", "R3").LoadArg("n", "R4")
+	b.Raw("loop:").Raw("CMP R4, $8").Raw("BLT tail").
+		Raw("VLM (R2), V0, V3"). // a
+		Raw("VLM (R3), V4, V7")  // b
+	for r := 0; r < 4; r++ {
+		b.Raw("%s V%d, V%d, V%d", vec, 4+r, r, r) // V_r = a OP b
+	}
+	b.Raw("VSTM V0, V3, (R1)").
+		Raw("ADD $64, R1").Raw("ADD $64, R2").Raw("ADD $64, R3").Raw("SUB $8, R4").Raw("BR loop").
+		Raw("tail:").Raw("CMP R4, $0").Raw("BEQ done").
+		Raw("FMOVD (R2), F0").Raw("FMOVD (R3), F1").
+		Raw(sca + " F1, F0").
+		Raw("FMOVD F0, (R1)").
+		Raw("ADD $8, R1").Raw("ADD $8, R2").Raw("ADD $8, R3").Raw("SUB $1, R4").Raw("BR tail").
+		Raw("done:").Ret()
+	return b.Func()
+}

internal/kernels/dot_generic.go

@@ -1,4 +1,4 @@
-//go:build !arm64 && !amd64 && !ppc64le && !loong64 && !riscv64
+//go:build !arm64 && !amd64 && !ppc64le && !loong64 && !riscv64 && !s390x
 
 package kernels
 

internal/kernels/fence_s390x_test.go

@@ -0,0 +1,9 @@
+package kernels
+
+// forEachFallback runs f again with the scalar kernels used without the vector facility.
+func forEachFallback(f func()) {
+	saved := useVX
+	useVX = false
+	defer func() { useVX = saved }()
+	f()
+}

internal/kernels/fence_test.go

@@ -1,4 +1,4 @@
-//go:build (amd64 || arm64 || ppc64le || loong64 || riscv64) && (linux || darwin)
+//go:build (amd64 || arm64 || ppc64le || loong64 || riscv64 || s390x) && (linux || darwin)
 
 package kernels
 

internal/kernels/guard_test.go

@@ -1,4 +1,4 @@
-//go:build amd64 || arm64 || ppc64le || loong64 || riscv64
+//go:build amd64 || arm64 || ppc64le || loong64 || riscv64 || s390x
 
 package kernels
 

internal/kernels/hwcap.go

@@ -7,12 +7,13 @@ import "encoding/binary"
 const atHWCAP = 16
 
 // hwcapFrom returns AT_HWCAP from the contents of /proc/self/auxv: pairs of
-// little-endian 64-bit words (key, value), ended by key 0. It returns 0 if the
+// 64-bit words in the machine's byte order (key, value), ended by key 0 --
+// big-endian on s390x, so the order must be native. It returns 0 if the
 // key is absent or the data is short. A pure function, so the parsing is
-// tested on every target even though only loong64 uses it today.
+// tested on every target even though only loong64, riscv64 and s390x use it.
 func hwcapFrom(auxv []byte) uint64 {
 	for len(auxv) >= 16 {
-		k, v := binary.LittleEndian.Uint64(auxv), binary.LittleEndian.Uint64(auxv[8:])
+		k, v := binary.NativeEndian.Uint64(auxv), binary.NativeEndian.Uint64(auxv[8:])
 		if k == 0 {
 			break
 		}

internal/kernels/hwcap_test.go

@@ -8,7 +8,7 @@ import (
 func auxvOf(pairs ...uint64) []byte {
 	b := make([]byte, 8*len(pairs))
 	for i, p := range pairs {
-		binary.LittleEndian.PutUint64(b[8*i:], p)
+		binary.NativeEndian.PutUint64(b[8*i:], p)
 	}
 	return b
 }

internal/kernels/reduce_generic.go

@@ -1,4 +1,4 @@
-//go:build !arm64 && !amd64 && !ppc64le && !loong64 && !riscv64
+//go:build !arm64 && !amd64 && !ppc64le && !loong64 && !riscv64 && !s390x
 
 package kernels
 
@@ -7,16 +7,15 @@ package kernels
 // SIMD-vs-scalar tests (which then trivially hold, scalar == scalar) still
 // exercise this path and the package builds uniformly across all six targets.
 //
-// amd64 (SSE2), arm64 (NEON), ppc64le (VSX), loong64 (LASX) and riscv64 (RVV),
-// the last two when the CPU has them, ship kernels (reduce_<arch>.go), so this
-// file is the s390x and 32-bit fallback. s390x's assembler has vector add/mul/
-// FMA, so its kernels are work not yet done; see docs/perf.md. The per-arch qemu jobs still
+// amd64 (SSE2), arm64 (NEON), ppc64le (VSX), loong64 (LASX), riscv64 (RVV) and
+// s390x (vector facility), the last three when the CPU has them, ship kernels
+// (reduce_<arch>.go), so this file is the 32-bit fallback. 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.
 
 // HaveReduceSIMD reports whether this build routes through a hand-vectorized
-// SIMD kernel (true on amd64/arm64/ppc64le/loong64/riscv64) rather than the
-// scalar oracle (false on s390x and the 32-bit targets). The kernels test logs it so each per-arch CI
+// SIMD kernel (true on every 64-bit target) rather than the scalar oracle
+// (false on the 32-bit targets). The kernels test logs it so each per-arch CI
 // run states which path it validated.
 const HaveReduceSIMD = false