simd(riscv64): RVV kernels for sum, dot, sqrt and add/sub/mul/div

.github/workflows/ci.yml

@@ -16,8 +16,9 @@ permissions:
 # amd64 SSE2 wrapper (reduce_amd64.go -> sum_amd64.s), the arm64 NEON wrapper
 # (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), or the scalar-alias wrapper (reduce_generic.go, on
-# riscv64/s390x and 32-bit). No single build sees
+# 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
+# single build sees
 # more than one, so no single run covers them all.
 #
 # The honest split:
@@ -41,7 +42,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/s390x keep the validated scalar path for now.
+# riscv64 ships RVV kernels (used when the CPU reports V); s390x keeps the
+# validated scalar path for now.
 jobs:
   test:
     name: test (${{ matrix.os }})
@@ -126,7 +128,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'
+        if: matrix.goarch == 'ppc64le' || matrix.goarch == 'loong64' || matrix.goarch == 'riscv64'
         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,15 @@ 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 and loong64** the hot loops are
-[go-asmgen](https://github.com/go-asmgen)-generated SIMD kernels: sum,
+`GOMAXPROCS`). On **amd64, arm64, ppc64le, loong64 and riscv64** 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 when the CPU has it), max/min too on
+arm64, VSX on ppc64le, LASX on loong64 and RVV on riscv64 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), the OpenBLAS/BLIS structure. The other 64-bit targets
-(riscv64, 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
+8×8; LASX 8×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
 (≤ 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 +119,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 and loong64 today, behind the same API and held to the
-same tests; riscv64 and s390x are next.
+on amd64, arm64, ppc64le, loong64 and riscv64 today, behind the same API and
+held to the same tests; s390x is next.
 
 ## 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 and loong64. Those kernels are the only code Go's
+for amd64, arm64, ppc64le, loong64 and riscv64. 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.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
+does the whole suite. On a riscv64 CPU without V (cfarm94, SiFive U74) the same
+build takes the scalar path, since V is gated on AT_HWCAP, and passes too. As a
+positive control for the fence on this hardware, `sumRVV` reading one element
+past the end faults (4 KiB pages). The
+tests were checked against two sabotaged kernels, a subtraction with swapped
+operands and a dot that does not advance `b`; both failed on the real hardware.
+
 **2026-10-07 (after v0.6.1, tests only).** `TestParallelStress` runs the
 parallel kernels from 6 goroutines at once while another keeps changing
 `GOMAXPROCS`, with the thresholds lowered so that tiny arrays go parallel. The

docs/perf.md

@@ -838,7 +838,26 @@ stay dynamic, so its slower efficiency cores are handled as before.
 - **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.
-- The other two 64-bit Go targets — **riscv64, s390x** — keep
+- **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. The GEMM micro-kernel is still the scalar 4×4. 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.
+- The remaining 64-bit Go target — **s390x** — keeps
   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
@@ -851,9 +870,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`), but the V extension is optional and needs a run-time
-  check. So s390x and riscv64 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. So s390x 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

docs/plan-ndarray.md

@@ -124,7 +124,8 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
     NumPy's pairwise sum). Since v0.2.3 ppc64le has VSX kernels too (sum,
     dot, sqrt, elementwise, an 8×8 GEMM micro-kernel), through go-asmgen's
     VSX encoders, and since v0.2.5 loong64 has LASX kernels (used when the CPU
-    reports LASX). The other targets (riscv64, s390x, 32-bit) keep the
+    reports LASX), and since v0.7.0 riscv64 has RVV kernels (sum, dot, sqrt,
+    elementwise; used when the CPU reports V). The other targets (s390x, 32-bit) 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/riscv64/gen.go

@@ -0,0 +1,127 @@
+// Command gen produces sum_riscv64.s, the RVV (vector extension) float64
+// kernels for riscv64, via go-asmgen. Run with: go run . (or `go generate`
+// from the kernels package).
+//
+// Go's riscv64 assembler names the vector instructions these kernels need
+// (VSETVLI, VLE64V/VSE64V, VFADDVV, VFSUBVV, VFMULVV, VFDIVVV, VFSQRTV,
+// VFMACCVV, VFREDUSUMVS, VFMVSF/VFMVFS), so nothing is emitted as a WORD. The
+// loops are vector-length agnostic: VSETVLI hands out as many elements as the
+// hardware vector holds (LMUL=8 groups eight registers), including a short last
+// strip, so there is no scalar tail. The V extension is optional, so the Go
+// side calls these only when the kernel reports it (AT_HWCAP).
+//
+// Operand order, decoded from Go's own testdata encodings: VFADDVV V1, V2, V3
+// computes V3 = V2 + V1 (likewise SUB and DIV: V3 = V2 - V1, V2 / V1);
+// VFMACCVV V1, V2, V3 computes V3 += V2*V1; VFREDUSUMVS Vs, Vv, Vd computes
+// Vd[0] = Vs[0] + sum(Vv).
+package main
+
+import (
+	"fmt"
+	"os"
+
+	"github.com/go-asmgen/asmgen/emit"
+	"github.com/go-asmgen/asmgen/riscv64"
+)
+
+func main() {
+	f := emit.NewFile("riscv64")
+	f.Add(reduceKernel("sumRVV", false))
+	f.Add(reduceKernel("dotRVV", true))
+	f.Add(unaryKernel("sqrtRVV", "VFSQRTV V0, V16"))
+	for _, k := range []struct{ name, op string }{
+		{"addRVV", "VFADDVV"}, {"subRVV", "VFSUBVV"},
+		{"mulRVV", "VFMULVV"}, {"divRVV", "VFDIVVV"},
+	} {
+		f.Add(binKernel(k.name, k.op))
+	}
+	if err := os.WriteFile("sum_riscv64.s", []byte(f.String()), 0o644); err != nil {
+		fmt.Fprintln(os.Stderr, err)
+		os.Exit(1)
+	}
+	fmt.Println("wrote sum_riscv64.s")
+}
+
+// advance moves each pointer register past vl (in X) float64s.
+func advance(b *riscv64.Builder, vl string, ptrs ...string) {
+	b.Raw("SLLI $3, %s, X17", vl)
+	for _, p := range ptrs {
+		b.Raw("ADD X17, %s, %s", p, p)
+	}
+}
+
+// reduceKernel builds NAME(a[, b] *float64, n int) float64 for n >= 1: the sum
+// of a, or with dot the sum of a[i]*b[i]. A vector accumulator of VLMAX lanes
+// is updated tail-undisturbed, so the lanes past a short last strip keep their
+// partial sums, then reduced to one scalar.
+func reduceKernel(name string, dot bool) *emit.Function {
+	args, types := []string{"a", "n"}, []riscv64.Type{riscv64.Ptr, riscv64.Int64}
+	if dot {
+		args, types = []string{"a", "b", "n"}, []riscv64.Type{riscv64.Ptr, riscv64.Ptr, riscv64.Int64}
+	}
+	sig := riscv64.Layout(args, types, []string{"ret"}, []riscv64.Type{riscv64.Float64})
+	b := riscv64.NewFunc(name, sig, 0)
+	b.LoadArg("a", "X10")
+	n := "X11"
+	if dot {
+		b.LoadArg("b", "X11")
+		n = "X12"
+	}
+	b.LoadArg("n", n)
+	b.Raw("VSETVLI X0, E64, M8, TA, MA, X13"). // X13 = VLMAX
+							Raw("VMVVI $0, V16") // accumulator = 0.0 in every lane
+	b.Raw("loop:").
+		Raw("VSETVLI %s, E64, M8, TU, MA, X13", n).
+		Raw("VLE64V (X10), V0")
+	if dot {
+		b.Raw("VLE64V (X11), V8").Raw("VFMACCVV V8, V0, V16") // V16 += a*b
+		advance(b, "X13", "X10", "X11")
+	} else {
+		b.Raw("VFADDVV V0, V16, V16") // V16 += a
+		advance(b, "X13", "X10")
+	}
+	b.Raw("SUB X13, %s, %s", n, n).
+		Raw("BNEZ %s, loop", n).
+		Raw("VSETVLI X0, E64, M8, TA, MA, X13").
+		Raw("FMVDX X0, F0").
+		Raw("VFMVSF F0, V24").
+		Raw("VFREDUSUMVS V24, V16, V24").
+		Raw("VFMVFS V24, F0").
+		StoreRet("F0", "ret").
+		Ret()
+	return b.Func()
+}
+
+// unaryKernel builds NAME(dst, src *float64, n int) for n >= 1: op reads V0 and
+// writes V16.
+func unaryKernel(name, op string) *emit.Function {
+	sig := riscv64.Layout([]string{"dst", "src", "n"},
+		[]riscv64.Type{riscv64.Ptr, riscv64.Ptr, riscv64.Int64}, nil, nil)
+	b := riscv64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "X10").LoadArg("src", "X11").LoadArg("n", "X13")
+	b.Raw("loop:").
+		Raw("VSETVLI X13, E64, M8, TA, MA, X14").
+		Raw("VLE64V (X11), V0").
+		Raw(op).
+		Raw("VSE64V V16, (X10)")
+	advance(b, "X14", "X10", "X11")
+	b.Raw("SUB X14, X13, X13").Raw("BNEZ X13, loop").Ret()
+	return b.Func()
+}
+
+// binKernel builds NAME(dst, a, b *float64, n int) for n >= 1: dst = a OP b.
+func binKernel(name, op string) *emit.Function {
+	sig := riscv64.Layout([]string{"dst", "a", "b", "n"},
+		[]riscv64.Type{riscv64.Ptr, riscv64.Ptr, riscv64.Ptr, riscv64.Int64}, nil, nil)
+	b := riscv64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "X10").LoadArg("a", "X11").LoadArg("b", "X12").LoadArg("n", "X13")
+	b.Raw("loop:").
+		Raw("VSETVLI X13, E64, M8, TA, MA, X14").
+		Raw("VLE64V (X11), V0").
+		Raw("VLE64V (X12), V8").
+		Raw("%s V8, V0, V16", op). // V16 = a OP b
+		Raw("VSE64V V16, (X10)")
+	advance(b, "X14", "X10", "X11", "X12")
+	b.Raw("SUB X14, X13, X13").Raw("BNEZ X13, loop").Ret()
+	return b.Func()
+}

internal/kernels/dot_generic.go

@@ -0,0 +1,6 @@
+//go:build !arm64 && !amd64 && !ppc64le && !loong64 && !riscv64
+
+package kernels
+
+// dotSIMD is dotRange here: no vector-double kernel on these targets.
+func dotSIMD(a, b []float64) float64 { return dotRange(a, b) }

internal/kernels/fence_riscv64_test.go

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

internal/kernels/fence_test.go

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

internal/kernels/gemm_generic.go

@@ -2,19 +2,16 @@
 
 package kernels
 
-// Generic GEMM micro-kernel for the targets without a vector-double SIMD kernel
-// (riscv64, s390x and the 32-bit ones): a portable scalar 4x4 register
-// tile over the packed panels. There is no hand-vectorized .s here yet (Go's
-// assembler has vector-double arithmetic for s390x and riscv64, see
-// docs/perf.md) — but the kernel still gets the two structural wins of the packed
+// Generic GEMM micro-kernel for the targets without a vector-double SIMD
+// micro-kernel (riscv64, s390x and the 32-bit ones): a portable scalar 4x4
+// register tile over the packed panels. There is no hand-vectorized tile here
+// yet (riscv64 has RVV kernels for the other loops; Go's assembler has
+// vector-double arithmetic for s390x too, 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
 // single-threaded numpy on large products.
 
-// dotSIMD is dotRange here: no vector-double kernel on these targets.
-func dotSIMD(a, b []float64) float64 { return dotRange(a, b) }
-
 // MR, NR are the register-block dimensions; 4x4 keeps 16 + a few scalars live in
 // the GP/FP register file on these targets without spilling.
 const (

internal/kernels/gemm_riscv64_test.go

@@ -0,0 +1,16 @@
+package kernels
+
+import "testing"
+
+// TestRVVDetected logs whether this machine runs the RVV kernels, and runs the
+// SIMD-vs-scalar checks again with them forced off.
+func TestRVVDetected(t *testing.T) {
+	t.Logf("useRVV = %v (AT_HWCAP)", useRVV)
+	saved := useRVV
+	useRVV = false
+	defer func() { useRVV = saved }()
+	t.Run("scalar", func(t *testing.T) {
+		TestBinSIMD(t)
+		TestSumSIMD(t)
+	})
+}

internal/kernels/guard_test.go

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

internal/kernels/reduce_generic.go

@@ -1,4 +1,4 @@
-//go:build !arm64 && !amd64 && !ppc64le && !loong64
+//go:build !arm64 && !amd64 && !ppc64le && !loong64 && !riscv64
 
 package kernels
 
@@ -7,17 +7,16 @@ 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) and loong64 (LASX, when the CPU
-// has it) ship kernels (reduce_<arch>.go), so this file is the riscv64/s390x
-// and 32-bit fallback. For those no hand-vectorized kernel is shipped yet:
-// s390x's assembler has vector add/mul/FMA and riscv64's V extension is
-// optional (a run-time check); see docs/perf.md. The per-arch qemu jobs still
+// 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
 // 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) rather than the scalar
-// oracle (false on riscv64/s390x and the 32-bit targets). The kernels test logs it so each per-arch CI
+// 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
 // run states which path it validated.
 const HaveReduceSIMD = false