From 106f5fc257fef4597f515031274502189f10ac3b Mon Sep 17 00:00:00 2001 From: tannevaled Date: Sun, 4 Oct 2026 18:14:01 +0200 Subject: [PATCH] perf: VSX kernels on ppc64le (sum, dot, sqrt, elementwise, 8x8 GEMM) 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 now encodes xvadddp/xvmaddadp/... as WORDs (pinned against GNU as, run bit-identical to Go's arithmetic and math.FMA on POWER8 and POWER9); this uses them for sum, dot, sqrt, add/sub/mul/div and an 8x8 GEMM micro-kernel (sum_ppc64le.s, generated). max/min stay scalar: xvmaxdp's NaN rule is not NumPy's. The GEMM kernel follows OpenBLAS's POWER9 dgemm scheme: no splats; A pairs and their swaps multiply B column pairs into diagonal/anti-diagonal accumulators, rebuilt into rows with xxpermdi at the store. Measured in one binary on a POWER9 core: 13.3 GFLOP/s with lxvdsx splats, 14.8 with xxpermdi splats, 15.4 unrolled, 16.9 with the diagonal scheme (shipped). POWER9 (cfarm29), against the scalar code, best of 3-4 interleaved runs: micro-kernel 8.3 -> 16.2 GFLOP/s; serial GEMM 512^2 7.9 -> 15.8 MatMul 1024^2 (32 threads) 39.2 -> 23.6 ms; 128^2 x1.8; 512^2 x1.5 MatVec x3.1; Sum x3.2-7.3 below 4 Mi; AddInto x2.7-4.1; Sqrt x1.05-1.18 Against OpenBLAS 0.3.34's POWER9 kernel on the same host: serial 0.69x, 32 threads 0.56x. Tests: the whole suite passes on a real POWER8 and POWER9; replacing xvsubdp's operand order or the store's xxpermdi selector makes it fail there. CI regenerates sum_ppc64le.s and fails if it is stale. Docs and comments that said ppc64le has no SIMD kernel are updated. Co-Authored-By: Claude Opus 5.5 --- .github/workflows/ci.yml | 19 +- README.md | 14 +- docs/perf.md | 46 +- internal/kernels/asmgen/go.mod | 2 +- internal/kernels/asmgen/go.sum | 4 +- internal/kernels/asmgen/ppc64le/gen.go | 284 +++++++++++ internal/kernels/gemm_generic.go | 12 +- internal/kernels/gemm_ppc64le.go | 35 ++ internal/kernels/kernels.go | 9 +- internal/kernels/reduce_generic.go | 20 +- internal/kernels/reduce_ppc64le.go | 79 +++ internal/kernels/sum_ppc64le.s | 652 +++++++++++++++++++++++++ 12 files changed, 1137 insertions(+), 39 deletions(-) create mode 100644 internal/kernels/asmgen/ppc64le/gen.go create mode 100644 internal/kernels/gemm_ppc64le.go create mode 100644 internal/kernels/reduce_ppc64le.go create mode 100644 internal/kernels/sum_ppc64le.s diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 68a6d6b..b504e5b 100644 --- a/.github/workflows/ci.yml +++ b/.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 ./... diff --git a/README.md b/README.md index b66ecc1..4853170 100644 --- a/README.md +++ b/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 diff --git a/docs/perf.md b/docs/perf.md index 93668ee..9be5b01 100644 --- a/docs/perf.md +++ b/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 diff --git a/internal/kernels/asmgen/go.mod b/internal/kernels/asmgen/go.mod index bafd420..156d43c 100644 --- a/internal/kernels/asmgen/go.mod +++ b/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 diff --git a/internal/kernels/asmgen/go.sum b/internal/kernels/asmgen/go.sum index 5813d72..1f76dfc 100644 --- a/internal/kernels/asmgen/go.sum +++ b/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= diff --git a/internal/kernels/asmgen/ppc64le/gen.go b/internal/kernels/asmgen/ppc64le/gen.go new file mode 100644 index 0000000..67e2b9a --- /dev/null +++ b/internal/kernels/asmgen/ppc64le/gen.go @@ -0,0 +1,284 @@ +// Command gen produces sum_ppc64le.s, the VSX float64 kernels for ppc64le, via +// go-asmgen. Run with: go run . (or `go generate` from the kernels package). +// +// Go's ppc64 assembler has the VSX loads, stores and splats (LXVD2X, LXVDSX, +// STXVD2X, XXPERMDI) but no vector double arithmetic; the arithmetic comes from +// go-asmgen's VSX encoders (Builder.XVADDDP, XVMADDADP, ...), which emit WORDs. +// +// Lane order: LXVD2X/STXVD2X swap the two doublewords of a vector on +// little-endian. Every vector of one kernel is loaded and stored that way, so a +// lane-wise operation pairs the same elements in every operand and the store +// puts each result back at its own index. LXVDSX (load and splat) puts the same +// double in both lanes, so it needs no care. +// +// Registers: R3.. hold the arguments; R14..R20 hold the byte offsets that +// LXVD2X's indexed form needs (R0 reads as 0); R30/R31/R13/R2/R1 are Go's and +// untouched. VS0..VS31 are accumulators, VS32..VS43 operands. +package main + +import ( + "fmt" + "os" + + "github.com/go-asmgen/asmgen/emit" + "github.com/go-asmgen/asmgen/ppc64" +) + +const ( + gemmMR = 8 + gemmNR = 8 +) + +func main() { + f := emit.NewFile("ppc64le") + f.Add(sumKernel()) + f.Add(dotKernel()) + f.Add(sqrtKernel()) + for _, k := range []struct { + name, vec, sca string + }{{"addVSX", "add", "FADD"}, {"subVSX", "sub", "FSUB"}, {"mulVSX", "mul", "FMUL"}, {"divVSX", "div", "FDIV"}} { + f.Add(binKernel(k.name, k.vec, k.sca)) + } + f.Add(gemmKernel()) + if err := os.WriteFile("sum_ppc64le.s", []byte(f.String()), 0o644); err != nil { + fmt.Fprintln(os.Stderr, err) + os.Exit(1) + } + fmt.Println("wrote sum_ppc64le.s") +} + +// offsets loads the byte offsets 16, 32, 48 into R14, R15, R16 (and, with +// splat, 8..56 into R17..R20 for the A-panel splats): LXVD2X/LXVDSX take only +// a (base)(index) address. +func offsets(b *ppc64.Builder, splat bool) { + b.Raw("MOVD $16, R14").Raw("MOVD $32, R15").Raw("MOVD $48, R16") + if splat { + b.Raw("MOVD $8, R17").Raw("MOVD $24, R18").Raw("MOVD $40, R19").Raw("MOVD $56, R20") + } +} + +// vidx is the index register for vector k (0..3) of a 64-byte block. +func vidx(k int) string { return [...]string{"R0", "R14", "R15", "R16"}[k] } + +// fold adds VS0..VS3 into VS0 and then its two lanes together, leaving the +// total in both lanes, so F0 (the scalar view of VS0) holds it. +func fold(b *ppc64.Builder) { + b.XVADDDP(0, 0, 1).XVADDDP(2, 2, 3).XVADDDP(0, 0, 2) + b.Raw("XXPERMDI VS0, VS0, $2, VS1") // swap the doublewords + b.XVADDDP(0, 0, 1) +} + +func reduceSig() ppc64.Signature { + return ppc64.Layout([]string{"a", "n"}, []ppc64.Type{ppc64.Ptr, ppc64.Int64}, + []string{"ret"}, []ppc64.Type{ppc64.Float64}) +} + +// sumKernel builds sumVSX(a *float64, n int) float64: four VSX accumulators +// (8 doubles per iteration), folded, plus a scalar tail. A regrouping of the +// sum, held to a tolerance like the other targets' SIMD sums. +func sumKernel() *emit.Function { + b := ppc64.NewFunc("sumVSX", reduceSig(), 0) + b.LoadArg("a", "R3").LoadArg("n", "R4") + offsets(b, false) + for r := 0; r < 4; r++ { + b.Raw("XXLXOR VS%d, VS%d, VS%d", r, r, r) + } + b.Label("sloop").Raw("CMP R4, $8").Raw("BLT sfold") + for k := 0; k < 4; k++ { + b.Raw("LXVD2X (R3)(%s), VS%d", vidx(k), 4+k) + } + for k := 0; k < 4; k++ { + b.XVADDDP(k, k, 4+k) + } + b.Raw("ADD $64, R3").Raw("ADD $-8, R4").Raw("BR sloop") + b.Label("sfold") + fold(b) + b.Label("stail").Raw("CMP R4, $0").Raw("BEQ sdone"). + Raw("FMOVD 0(R3), F1").Raw("FADD F1, F0"). + Raw("ADD $8, R3").Raw("ADD $-1, R4").Raw("BR stail") + b.Label("sdone").StoreRet("F0", "ret").Ret() + return b.Func() +} + +// dotKernel builds dotVSX(a, b *float64, n int) float64: four fused +// multiply-add accumulators (XVMADDADP), folded, plus a scalar tail. +func dotKernel() *emit.Function { + sig := ppc64.Layout([]string{"a", "b", "n"}, []ppc64.Type{ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, + []string{"ret"}, []ppc64.Type{ppc64.Float64}) + b := ppc64.NewFunc("dotVSX", sig, 0) + b.LoadArg("a", "R3").LoadArg("b", "R5").LoadArg("n", "R4") + offsets(b, false) + for r := 0; r < 4; r++ { + b.Raw("XXLXOR VS%d, VS%d, VS%d", r, r, r) + } + b.Label("dloop").Raw("CMP R4, $8").Raw("BLT dfold") + for k := 0; k < 4; k++ { + b.Raw("LXVD2X (R3)(%s), VS%d", vidx(k), 32+k) + b.Raw("LXVD2X (R5)(%s), VS%d", vidx(k), 36+k) + } + for k := 0; k < 4; k++ { + b.XVMADDADP(k, 32+k, 36+k) + } + b.Raw("ADD $64, R3").Raw("ADD $64, R5").Raw("ADD $-8, R4").Raw("BR dloop") + b.Label("dfold") + fold(b) + b.Label("dtail").Raw("CMP R4, $0").Raw("BEQ ddone"). + Raw("FMOVD 0(R3), F1").Raw("FMOVD 0(R5), F2").Raw("FMUL F2, F1").Raw("FADD F1, F0"). + Raw("ADD $8, R3").Raw("ADD $8, R5").Raw("ADD $-1, R4").Raw("BR dtail") + b.Label("ddone").StoreRet("F0", "ret").Ret() + return b.Func() +} + +// sqrtKernel builds sqrtVSX(dst, src *float64, n int): XVSQRTDP over 8 doubles +// per iteration, FSQRT for the tail. Both are the correctly rounded IEEE +// square root, so the result is bit-identical to math.Sqrt. +func sqrtKernel() *emit.Function { + sig := ppc64.Layout([]string{"dst", "src", "n"}, []ppc64.Type{ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, nil, nil) + b := ppc64.NewFunc("sqrtVSX", sig, 0) + b.LoadArg("dst", "R3").LoadArg("src", "R5").LoadArg("n", "R4") + offsets(b, false) + b.Label("qloop").Raw("CMP R4, $8").Raw("BLT qtail") + for k := 0; k < 4; k++ { + b.Raw("LXVD2X (R5)(%s), VS%d", vidx(k), 32+k) + } + for k := 0; k < 4; k++ { + b.XVSQRTDP(36+k, 32+k) + } + for k := 0; k < 4; k++ { + b.Raw("STXVD2X VS%d, (R3)(%s)", 36+k, vidx(k)) + } + b.Raw("ADD $64, R3").Raw("ADD $64, R5").Raw("ADD $-8, R4").Raw("BR qloop") + b.Label("qtail").Raw("CMP R4, $0").Raw("BEQ qdone"). + Raw("FMOVD 0(R5), F1").Raw("FSQRT F1, F1").Raw("FMOVD F1, 0(R3)"). + Raw("ADD $8, R3").Raw("ADD $8, R5").Raw("ADD $-1, R4").Raw("BR qtail") + b.Label("qdone").Ret() + return b.Func() +} + +// binKernel builds name(dst, a, b *float64, n int): dst[i] = a[i] OP b[i], 8 +// doubles per iteration with the VSX op and a scalar tail with sca. Each lane +// is one correctly rounded IEEE operation, bit-identical to the scalar oracle. +func binKernel(name, op, sca string) *emit.Function { + sig := ppc64.Layout([]string{"dst", "a", "b", "n"}, + []ppc64.Type{ppc64.Ptr, ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, nil, nil) + b := ppc64.NewFunc(name, sig, 0) + b.LoadArg("dst", "R3").LoadArg("a", "R5").LoadArg("b", "R6").LoadArg("n", "R4") + offsets(b, false) + l := name[:3] + b.Label(l+"loop").Raw("CMP R4, $8").Raw("BLT %stail", l) + for k := 0; k < 4; k++ { + b.Raw("LXVD2X (R5)(%s), VS%d", vidx(k), 32+k) + b.Raw("LXVD2X (R6)(%s), VS%d", vidx(k), 36+k) + } + for k := 0; k < 4; k++ { + switch op { + case "add": + b.XVADDDP(40+k, 32+k, 36+k) + case "sub": + b.XVSUBDP(40+k, 32+k, 36+k) + case "mul": + b.XVMULDP(40+k, 32+k, 36+k) + case "div": + b.XVDIVDP(40+k, 32+k, 36+k) + } + } + for k := 0; k < 4; k++ { + b.Raw("STXVD2X VS%d, (R3)(%s)", 40+k, vidx(k)) + } + b.Raw("ADD $64, R3").Raw("ADD $64, R5").Raw("ADD $64, R6").Raw("ADD $-8, R4").Raw("BR %sloop", l) + b.Label(l+"tail").Raw("CMP R4, $0").Raw("BEQ %sdone", l). + Raw("FMOVD 0(R5), F1").Raw("FMOVD 0(R6), F2").Raw("%s F2, F1", sca).Raw("FMOVD F1, 0(R3)"). + Raw("ADD $8, R3").Raw("ADD $8, R5").Raw("ADD $8, R6").Raw("ADD $-1, R4").Raw("BR %stail", l) + b.Label(l + "done").Ret() + return b.Func() +} + +// gemmKernel builds gemmMicro8x8VSX(kc int, pa, pb, c *float64, ldc int): +// +// for p in [0,kc): C[r][col] += pa[p*8+r] * pb[p*8+col] (r<8, col<8) +// +// with the scheme of OpenBLAS's POWER9 dgemm kernel: no splats. Each A pair +// (a_2p, a_2p+1) and its swap (one XXPERMDI) multiply each B column pair, so +// accumulator VS(8p+2k) holds the diagonal (C[2p][2k], C[2p+1][2k+1]) and +// VS(8p+2k+1) the anti-diagonal (C[2p+1][2k], C[2p][2k+1]); the store rebuilds +// the rows with XXPERMDI (DM=1). 32 fused multiply-adds per k step against 8 +// loads and 4 swaps; two k steps per iteration, the second step's loads issued +// between the first step's FMAs. Measured on a POWER9 core (cfarm29): 16.9 +// GFLOP/s, against 13.3 with LXVDSX splats, 14.8 with XXPERMDI splats and 15.4 +// with those unrolled; OpenBLAS's own kernel reaches 23 there. +func gemmKernel() *emit.Function { + name := "gemmMicro8x8VSX" + sig := ppc64.Layout([]string{"kc", "pa", "pb", "c", "ldc"}, + []ppc64.Type{ppc64.Int64, ppc64.Ptr, ppc64.Ptr, ppc64.Ptr, ppc64.Int64}, nil, nil) + b := ppc64.NewFunc(name, sig, 0) + b.LoadArg("kc", "R3").LoadArg("pa", "R4").LoadArg("pb", "R5").LoadArg("c", "R6").LoadArg("ldc", "R7") + b.Raw("SLD $3, R7") + offsets(b, false) + for r := 0; r < 32; r++ { + b.Raw("XXLXOR VS%d, VS%d, VS%d", r, r, r) + } + b.Raw("CMP R3, $0").Raw("BEQ %sstore", name) + // step uses B vectors bv..bv+3, A pairs av..av+3, swapped pairs sw..sw+3. + loadsOf := func(bv, av int) []string { + var l []string + for k := 0; k < 4; k++ { + l = append(l, fmt.Sprintf("LXVD2X (R5)(%s), VS%d", vidx(k), bv+k), fmt.Sprintf("LXVD2X (R4)(%s), VS%d", vidx(k), av+k)) + } + return l + } + step := func(bv, av, sw int, pending []string) { + for p := 0; p < 4; p++ { + b.Raw("XXPERMDI VS%d, VS%d, $2, VS%d", av+p, av+p, sw+p) + } + i := 0 + for p := 0; p < 4; p++ { + for k := 0; k < 4; k++ { + b.XVMADDADP(8*p+2*k, av+p, bv+k) + b.XVMADDADP(8*p+2*k+1, sw+p, bv+k) + if i%2 == 1 && len(pending) > 0 { + b.Raw("%s", pending[0]) + pending = pending[1:] + } + i++ + } + } + for _, l := range pending { // never leave a load out + b.Raw("%s", l) + } + } + b.Label(name+"loop2").Raw("CMP R3, $2").Raw("BLT %sone", name) + for _, l := range loadsOf(32, 36) { + b.Raw("%s", l) + } + b.Raw("ADD $64, R4").Raw("ADD $64, R5") + step(32, 36, 40, loadsOf(44, 48)) + b.Raw("ADD $64, R4").Raw("ADD $64, R5") + step(44, 48, 52, nil) + b.Raw("ADD $-2, R3").Raw("BR %sloop2", name) + b.Label(name+"one").Raw("CMP R3, $0").Raw("BEQ %sstore", name) + for _, l := range loadsOf(32, 36) { + b.Raw("%s", l) + } + step(32, 36, 40, nil) + b.Label(name + "store") + for p := 0; p < 4; p++ { + for half := 0; half < 2; half++ { // row 2p, then 2p+1 + for k := 0; k < 4; k++ { + d, a := 8*p+2*k, 8*p+2*k+1 + if half == 0 { + b.Raw("XXPERMDI VS%d, VS%d, $1, VS%d", d, a, 56) // [diag.dw0, anti.dw1] + } else { + b.Raw("XXPERMDI VS%d, VS%d, $1, VS%d", a, d, 56) // [anti.dw0, diag.dw1] + } + b.Raw("LXVD2X (R6)(%s), VS%d", vidx(k), 57) + b.XVADDDP(57, 57, 56) + b.Raw("STXVD2X VS%d, (R6)(%s)", 57, vidx(k)) + } + if !(p == 3 && half == 1) { + b.Raw("ADD R7, R6") + } + } + } + b.Ret() + return b.Func() +} diff --git a/internal/kernels/gemm_generic.go b/internal/kernels/gemm_generic.go index ba202f5..e91b246 100644 --- a/internal/kernels/gemm_generic.go +++ b/internal/kernels/gemm_generic.go @@ -1,12 +1,12 @@ -//go:build !arm64 && !amd64 +//go:build !arm64 && !amd64 && !ppc64le 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 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 +// Generic GEMM micro-kernel for the targets without a vector-double SIMD kernel +// (loong64, 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 loong64, s390x and riscv64, 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/gemm_ppc64le.go b/internal/kernels/gemm_ppc64le.go new file mode 100644 index 0000000..fc53ee9 --- /dev/null +++ b/internal/kernels/gemm_ppc64le.go @@ -0,0 +1,35 @@ +package kernels + +// ppc64le GEMM micro-kernel: an 8x8 VSX tile (gemmMicro8x8VSX in +// sum_ppc64le.s, generated by go-asmgen): 32 two-lane accumulators, the B row +// loaded with LXVD2X, each A value splatted with LXVDSX, and a fused +// XVMADDADP per (row, column pair). Fused like the arm64 and amd64 FMA +// kernels: on integer-valued data it equals the scalar ikj oracle exactly. + +// MR, NR are the micro-kernel's register-block dimensions: it computes an +// MR-row x NR-col tile of C per call. The packing in gemm.go matches. +const ( + MR = 8 + NR = 8 +) + +// gemmMicro adds the MR x NR tile sum_p pa[p*MR..]*pb[p*NR..] into the C block +// at dst[0:], whose rows are ldc apart. +func gemmMicro(kc int, pa, pb, dst []float64, ldc int) { + gemmMicro8x8VSX(kc, &pa[0], &pb[0], &dst[0], ldc) +} + +// dotSIMD is the inner product of two equal-length slices (VSX kernel, four +// fused multiply-add accumulators). +func dotSIMD(a, b []float64) float64 { + if len(a) == 0 { + return 0 + } + return dotVSX(&a[0], &b[0], len(a)) +} + +//go:noescape +func gemmMicro8x8VSX(kc int, pa, pb, c *float64, ldc int) + +//go:noescape +func dotVSX(a, b *float64, n int) float64 diff --git a/internal/kernels/kernels.go b/internal/kernels/kernels.go index a422ee3..984ae6b 100644 --- a/internal/kernels/kernels.go +++ b/internal/kernels/kernels.go @@ -1,11 +1,10 @@ // Package kernels holds the inner numeric loops for the ndarray package. // // Every loop here is a portable, pure-Go (CGO=0) scalar kernel. They are kept -// behind this small, contiguous-slice API so that SIMD variants can drop in -// later without changing callers or tests: Phase 1 of the roadmap replaces -// these with go-asmgen-generated kernels across all six 64-bit Go SIMD targets -// (amd64, arm64, riscv64, loong64, ppc64le, s390x), selected at runtime, while -// these scalar versions remain the reference and the fallback. +// behind this small, contiguous-slice API so that SIMD variants drop in without +// changing callers or tests: go-asmgen-generated kernels replace the hot ones +// on amd64, arm64 and ppc64le (the *_amd64/_arm64/_ppc64le files), while these +// scalar versions remain the reference and the fallback elsewhere. // // Kernels operate on flat, contiguous []float64 slices. Shape, stride and // broadcasting concerns live in the parent package; by the time a slice reaches diff --git a/internal/kernels/reduce_generic.go b/internal/kernels/reduce_generic.go index 6a3d38a..f204162 100644 --- a/internal/kernels/reduce_generic.go +++ b/internal/kernels/reduce_generic.go @@ -1,4 +1,4 @@ -//go:build !arm64 && !amd64 +//go:build !arm64 && !amd64 && !ppc64le package kernels @@ -7,20 +7,20 @@ 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) 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 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 +// amd64 (SSE2), arm64 (NEON) and ppc64le (VSX, through go-asmgen's encoders) +// ship kernels (reduce_amd64.go / reduce_arm64.go / reduce_ppc64le.go), so this +// file is the loong64/riscv64/s390x and 32-bit fallback. For those no +// hand-vectorized kernel is shipped yet: 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. // HaveReduceSIMD reports whether this build routes through a hand-vectorized -// SIMD kernel (true on amd64/arm64) rather than the scalar oracle (false on -// loong64/ppc64le/riscv64/s390x). The kernels test logs it so each per-arch CI +// SIMD kernel (true on amd64/arm64/ppc64le) rather than the scalar oracle +// (false on loong64/riscv64/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 diff --git a/internal/kernels/reduce_ppc64le.go b/internal/kernels/reduce_ppc64le.go new file mode 100644 index 0000000..c62df0c --- /dev/null +++ b/internal/kernels/reduce_ppc64le.go @@ -0,0 +1,79 @@ +package kernels + +// ppc64le VSX kernels generated by go-asmgen (sum_ppc64le.s). VSX is part of +// the ppc64le baseline (POWER8, GOPPC64=power8), so they are always callable. +// Go's ppc64 assembler has no vector double arithmetic; go-asmgen encodes it +// (xvadddp, xvmaddadp, ...) as WORDs, checked against GNU as. +// +//go:generate sh -c "cd asmgen/ppc64le && go run . && mv sum_ppc64le.s ../../sum_ppc64le.s" + +// HaveReduceSIMD reports that this build ships SIMD reduction/elementwise +// kernels. +const HaveReduceSIMD = true + +// sumSIMD returns the VSX sum of a (empty -> 0): four accumulators, a valid +// regrouping of the sum, like the other targets' SIMD sums. +func sumSIMD(a []float64) float64 { + if len(a) == 0 { + return 0 + } + return sumVSX(&a[0], len(a)) +} + +// sqrtSIMD writes sqrt(src[i]) into dst[i] with XVSQRTDP, bit-identical to +// math.Sqrt (one correctly rounded operation per lane). +func sqrtSIMD(dst, src []float64) { + if len(dst) == 0 { + return + } + sqrtVSX(&dst[0], &src[0], len(dst)) +} + +// max/min stay on the four-accumulator scalar reducers: XVMAXDP/XVMINDP do not +// propagate NaN the way Go's max/min and numpy.max do. +func maxSIMD(a []float64) float64 { return maxUnrolled(a) } +func minSIMD(a []float64) float64 { return minUnrolled(a) } + +// addBin, subBin, mulBin, divBin write dst[i] = a[i] OP b[i] with the VSX +// kernels, each lane one IEEE operation, bit-identical to the scalar oracle. +func addBin(dst, a, b []float64) { + if len(dst) > 0 { + addVSX(&dst[0], &a[0], &b[0], len(dst)) + } +} + +func subBin(dst, a, b []float64) { + if len(dst) > 0 { + subVSX(&dst[0], &a[0], &b[0], len(dst)) + } +} + +func mulBin(dst, a, b []float64) { + if len(dst) > 0 { + mulVSX(&dst[0], &a[0], &b[0], len(dst)) + } +} + +func divBin(dst, a, b []float64) { + if len(dst) > 0 { + divVSX(&dst[0], &a[0], &b[0], len(dst)) + } +} + +//go:noescape +func sumVSX(a *float64, n int) float64 + +//go:noescape +func sqrtVSX(dst, src *float64, n int) + +//go:noescape +func addVSX(dst, a, b *float64, n int) + +//go:noescape +func subVSX(dst, a, b *float64, n int) + +//go:noescape +func mulVSX(dst, a, b *float64, n int) + +//go:noescape +func divVSX(dst, a, b *float64, n int) diff --git a/internal/kernels/sum_ppc64le.s b/internal/kernels/sum_ppc64le.s new file mode 100644 index 0000000..88ddca2 --- /dev/null +++ b/internal/kernels/sum_ppc64le.s @@ -0,0 +1,652 @@ +// Code generated by go-asmgen. DO NOT EDIT. + +//go:build ppc64le + +#include "textflag.h" + +TEXT ·sumVSX(SB), NOSPLIT, $0-24 + MOVD a+0(FP), R3 + MOVD n+8(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 + XXLXOR VS0, VS0, VS0 + XXLXOR VS1, VS1, VS1 + XXLXOR VS2, VS2, VS2 + XXLXOR VS3, VS3, VS3 +sloop: + CMP R4, $8 + BLT sfold + LXVD2X (R3)(R0), VS4 + LXVD2X (R3)(R14), VS5 + LXVD2X (R3)(R15), VS6 + LXVD2X (R3)(R16), VS7 + WORD $0xf0002300 // xvadddp vs0, vs0, vs4 + WORD $0xf0212b00 // xvadddp vs1, vs1, vs5 + WORD $0xf0423300 // xvadddp vs2, vs2, vs6 + WORD $0xf0633b00 // xvadddp vs3, vs3, vs7 + ADD $64, R3 + ADD $-8, R4 + BR sloop +sfold: + WORD $0xf0000b00 // xvadddp vs0, vs0, vs1 + WORD $0xf0421b00 // xvadddp vs2, vs2, vs3 + WORD $0xf0001300 // xvadddp vs0, vs0, vs2 + XXPERMDI VS0, VS0, $2, VS1 + WORD $0xf0000b00 // xvadddp vs0, vs0, vs1 +stail: + CMP R4, $0 + BEQ sdone + FMOVD 0(R3), F1 + FADD F1, F0 + ADD $8, R3 + ADD $-1, R4 + BR stail +sdone: + FMOVD F0, ret+16(FP) + RET + +TEXT ·dotVSX(SB), NOSPLIT, $0-32 + MOVD a+0(FP), R3 + MOVD b+8(FP), R5 + MOVD n+16(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 + XXLXOR VS0, VS0, VS0 + XXLXOR VS1, VS1, VS1 + XXLXOR VS2, VS2, VS2 + XXLXOR VS3, VS3, VS3 +dloop: + CMP R4, $8 + BLT dfold + LXVD2X (R3)(R0), VS32 + LXVD2X (R5)(R0), VS36 + LXVD2X (R3)(R14), VS33 + LXVD2X (R5)(R14), VS37 + LXVD2X (R3)(R15), VS34 + LXVD2X (R5)(R15), VS38 + LXVD2X (R3)(R16), VS35 + LXVD2X (R5)(R16), VS39 + WORD $0xf000230e // xvmaddadp vs0, vs32, vs36 + WORD $0xf0212b0e // xvmaddadp vs1, vs33, vs37 + WORD $0xf042330e // xvmaddadp vs2, vs34, vs38 + WORD $0xf0633b0e // xvmaddadp vs3, vs35, vs39 + ADD $64, R3 + ADD $64, R5 + ADD $-8, R4 + BR dloop +dfold: + WORD $0xf0000b00 // xvadddp vs0, vs0, vs1 + WORD $0xf0421b00 // xvadddp vs2, vs2, vs3 + WORD $0xf0001300 // xvadddp vs0, vs0, vs2 + XXPERMDI VS0, VS0, $2, VS1 + WORD $0xf0000b00 // xvadddp vs0, vs0, vs1 +dtail: + CMP R4, $0 + BEQ ddone + FMOVD 0(R3), F1 + FMOVD 0(R5), F2 + FMUL F2, F1 + FADD F1, F0 + ADD $8, R3 + ADD $8, R5 + ADD $-1, R4 + BR dtail +ddone: + FMOVD F0, ret+24(FP) + RET + +TEXT ·sqrtVSX(SB), NOSPLIT, $0-24 + MOVD dst+0(FP), R3 + MOVD src+8(FP), R5 + MOVD n+16(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 +qloop: + CMP R4, $8 + BLT qtail + LXVD2X (R5)(R0), VS32 + LXVD2X (R5)(R14), VS33 + LXVD2X (R5)(R15), VS34 + LXVD2X (R5)(R16), VS35 + WORD $0xf080032f // xvsqrtdp vs36, vs32 + WORD $0xf0a00b2f // xvsqrtdp vs37, vs33 + WORD $0xf0c0132f // xvsqrtdp vs38, vs34 + WORD $0xf0e01b2f // xvsqrtdp vs39, vs35 + STXVD2X VS36, (R3)(R0) + STXVD2X VS37, (R3)(R14) + STXVD2X VS38, (R3)(R15) + STXVD2X VS39, (R3)(R16) + ADD $64, R3 + ADD $64, R5 + ADD $-8, R4 + BR qloop +qtail: + CMP R4, $0 + BEQ qdone + FMOVD 0(R5), F1 + FSQRT F1, F1 + FMOVD F1, 0(R3) + ADD $8, R3 + ADD $8, R5 + ADD $-1, R4 + BR qtail +qdone: + RET + +TEXT ·addVSX(SB), NOSPLIT, $0-32 + MOVD dst+0(FP), R3 + MOVD a+8(FP), R5 + MOVD b+16(FP), R6 + MOVD n+24(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 +addloop: + CMP R4, $8 + BLT addtail + LXVD2X (R5)(R0), VS32 + LXVD2X (R6)(R0), VS36 + LXVD2X (R5)(R14), VS33 + LXVD2X (R6)(R14), VS37 + LXVD2X (R5)(R15), VS34 + LXVD2X (R6)(R15), VS38 + LXVD2X (R5)(R16), VS35 + LXVD2X (R6)(R16), VS39 + WORD $0xf1002307 // xvadddp vs40, vs32, vs36 + WORD $0xf1212b07 // xvadddp vs41, vs33, vs37 + WORD $0xf1423307 // xvadddp vs42, vs34, vs38 + WORD $0xf1633b07 // xvadddp vs43, vs35, vs39 + STXVD2X VS40, (R3)(R0) + STXVD2X VS41, (R3)(R14) + STXVD2X VS42, (R3)(R15) + STXVD2X VS43, (R3)(R16) + ADD $64, R3 + ADD $64, R5 + ADD $64, R6 + ADD $-8, R4 + BR addloop +addtail: + CMP R4, $0 + BEQ adddone + FMOVD 0(R5), F1 + FMOVD 0(R6), F2 + FADD F2, F1 + FMOVD F1, 0(R3) + ADD $8, R3 + ADD $8, R5 + ADD $8, R6 + ADD $-1, R4 + BR addtail +adddone: + RET + +TEXT ·subVSX(SB), NOSPLIT, $0-32 + MOVD dst+0(FP), R3 + MOVD a+8(FP), R5 + MOVD b+16(FP), R6 + MOVD n+24(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 +subloop: + CMP R4, $8 + BLT subtail + LXVD2X (R5)(R0), VS32 + LXVD2X (R6)(R0), VS36 + LXVD2X (R5)(R14), VS33 + LXVD2X (R6)(R14), VS37 + LXVD2X (R5)(R15), VS34 + LXVD2X (R6)(R15), VS38 + LXVD2X (R5)(R16), VS35 + LXVD2X (R6)(R16), VS39 + WORD $0xf1002347 // xvsubdp vs40, vs32, vs36 + WORD $0xf1212b47 // xvsubdp vs41, vs33, vs37 + WORD $0xf1423347 // xvsubdp vs42, vs34, vs38 + WORD $0xf1633b47 // xvsubdp vs43, vs35, vs39 + STXVD2X VS40, (R3)(R0) + STXVD2X VS41, (R3)(R14) + STXVD2X VS42, (R3)(R15) + STXVD2X VS43, (R3)(R16) + ADD $64, R3 + ADD $64, R5 + ADD $64, R6 + ADD $-8, R4 + BR subloop +subtail: + CMP R4, $0 + BEQ subdone + FMOVD 0(R5), F1 + FMOVD 0(R6), F2 + FSUB F2, F1 + FMOVD F1, 0(R3) + ADD $8, R3 + ADD $8, R5 + ADD $8, R6 + ADD $-1, R4 + BR subtail +subdone: + RET + +TEXT ·mulVSX(SB), NOSPLIT, $0-32 + MOVD dst+0(FP), R3 + MOVD a+8(FP), R5 + MOVD b+16(FP), R6 + MOVD n+24(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 +mulloop: + CMP R4, $8 + BLT multail + LXVD2X (R5)(R0), VS32 + LXVD2X (R6)(R0), VS36 + LXVD2X (R5)(R14), VS33 + LXVD2X (R6)(R14), VS37 + LXVD2X (R5)(R15), VS34 + LXVD2X (R6)(R15), VS38 + LXVD2X (R5)(R16), VS35 + LXVD2X (R6)(R16), VS39 + WORD $0xf1002387 // xvmuldp vs40, vs32, vs36 + WORD $0xf1212b87 // xvmuldp vs41, vs33, vs37 + WORD $0xf1423387 // xvmuldp vs42, vs34, vs38 + WORD $0xf1633b87 // xvmuldp vs43, vs35, vs39 + STXVD2X VS40, (R3)(R0) + STXVD2X VS41, (R3)(R14) + STXVD2X VS42, (R3)(R15) + STXVD2X VS43, (R3)(R16) + ADD $64, R3 + ADD $64, R5 + ADD $64, R6 + ADD $-8, R4 + BR mulloop +multail: + CMP R4, $0 + BEQ muldone + FMOVD 0(R5), F1 + FMOVD 0(R6), F2 + FMUL F2, F1 + FMOVD F1, 0(R3) + ADD $8, R3 + ADD $8, R5 + ADD $8, R6 + ADD $-1, R4 + BR multail +muldone: + RET + +TEXT ·divVSX(SB), NOSPLIT, $0-32 + MOVD dst+0(FP), R3 + MOVD a+8(FP), R5 + MOVD b+16(FP), R6 + MOVD n+24(FP), R4 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 +divloop: + CMP R4, $8 + BLT divtail + LXVD2X (R5)(R0), VS32 + LXVD2X (R6)(R0), VS36 + LXVD2X (R5)(R14), VS33 + LXVD2X (R6)(R14), VS37 + LXVD2X (R5)(R15), VS34 + LXVD2X (R6)(R15), VS38 + LXVD2X (R5)(R16), VS35 + LXVD2X (R6)(R16), VS39 + WORD $0xf10023c7 // xvdivdp vs40, vs32, vs36 + WORD $0xf1212bc7 // xvdivdp vs41, vs33, vs37 + WORD $0xf14233c7 // xvdivdp vs42, vs34, vs38 + WORD $0xf1633bc7 // xvdivdp vs43, vs35, vs39 + STXVD2X VS40, (R3)(R0) + STXVD2X VS41, (R3)(R14) + STXVD2X VS42, (R3)(R15) + STXVD2X VS43, (R3)(R16) + ADD $64, R3 + ADD $64, R5 + ADD $64, R6 + ADD $-8, R4 + BR divloop +divtail: + CMP R4, $0 + BEQ divdone + FMOVD 0(R5), F1 + FMOVD 0(R6), F2 + FDIV F2, F1 + FMOVD F1, 0(R3) + ADD $8, R3 + ADD $8, R5 + ADD $8, R6 + ADD $-1, R4 + BR divtail +divdone: + RET + +TEXT ·gemmMicro8x8VSX(SB), NOSPLIT, $0-40 + MOVD kc+0(FP), R3 + MOVD pa+8(FP), R4 + MOVD pb+16(FP), R5 + MOVD c+24(FP), R6 + MOVD ldc+32(FP), R7 + SLD $3, R7 + MOVD $16, R14 + MOVD $32, R15 + MOVD $48, R16 + XXLXOR VS0, VS0, VS0 + XXLXOR VS1, VS1, VS1 + XXLXOR VS2, VS2, VS2 + XXLXOR VS3, VS3, VS3 + XXLXOR VS4, VS4, VS4 + XXLXOR VS5, VS5, VS5 + XXLXOR VS6, VS6, VS6 + XXLXOR VS7, VS7, VS7 + XXLXOR VS8, VS8, VS8 + XXLXOR VS9, VS9, VS9 + XXLXOR VS10, VS10, VS10 + XXLXOR VS11, VS11, VS11 + XXLXOR VS12, VS12, VS12 + XXLXOR VS13, VS13, VS13 + XXLXOR VS14, VS14, VS14 + XXLXOR VS15, VS15, VS15 + XXLXOR VS16, VS16, VS16 + XXLXOR VS17, VS17, VS17 + XXLXOR VS18, VS18, VS18 + XXLXOR VS19, VS19, VS19 + XXLXOR VS20, VS20, VS20 + XXLXOR VS21, VS21, VS21 + XXLXOR VS22, VS22, VS22 + XXLXOR VS23, VS23, VS23 + XXLXOR VS24, VS24, VS24 + XXLXOR VS25, VS25, VS25 + XXLXOR VS26, VS26, VS26 + XXLXOR VS27, VS27, VS27 + XXLXOR VS28, VS28, VS28 + XXLXOR VS29, VS29, VS29 + XXLXOR VS30, VS30, VS30 + XXLXOR VS31, VS31, VS31 + CMP R3, $0 + BEQ gemmMicro8x8VSXstore +gemmMicro8x8VSXloop2: + CMP R3, $2 + BLT gemmMicro8x8VSXone + LXVD2X (R5)(R0), VS32 + LXVD2X (R4)(R0), VS36 + LXVD2X (R5)(R14), VS33 + LXVD2X (R4)(R14), VS37 + LXVD2X (R5)(R15), VS34 + LXVD2X (R4)(R15), VS38 + LXVD2X (R5)(R16), VS35 + LXVD2X (R4)(R16), VS39 + ADD $64, R4 + ADD $64, R5 + XXPERMDI VS36, VS36, $2, VS40 + XXPERMDI VS37, VS37, $2, VS41 + XXPERMDI VS38, VS38, $2, VS42 + XXPERMDI VS39, VS39, $2, VS43 + WORD $0xf004030e // xvmaddadp vs0, vs36, vs32 + WORD $0xf028030e // xvmaddadp vs1, vs40, vs32 + WORD $0xf0440b0e // xvmaddadp vs2, vs36, vs33 + WORD $0xf0680b0e // xvmaddadp vs3, vs40, vs33 + LXVD2X (R5)(R0), VS44 + WORD $0xf084130e // xvmaddadp vs4, vs36, vs34 + WORD $0xf0a8130e // xvmaddadp vs5, vs40, vs34 + WORD $0xf0c41b0e // xvmaddadp vs6, vs36, vs35 + WORD $0xf0e81b0e // xvmaddadp vs7, vs40, vs35 + LXVD2X (R4)(R0), VS48 + WORD $0xf105030e // xvmaddadp vs8, vs37, vs32 + WORD $0xf129030e // xvmaddadp vs9, vs41, vs32 + WORD $0xf1450b0e // xvmaddadp vs10, vs37, vs33 + WORD $0xf1690b0e // xvmaddadp vs11, vs41, vs33 + LXVD2X (R5)(R14), VS45 + WORD $0xf185130e // xvmaddadp vs12, vs37, vs34 + WORD $0xf1a9130e // xvmaddadp vs13, vs41, vs34 + WORD $0xf1c51b0e // xvmaddadp vs14, vs37, vs35 + WORD $0xf1e91b0e // xvmaddadp vs15, vs41, vs35 + LXVD2X (R4)(R14), VS49 + WORD $0xf206030e // xvmaddadp vs16, vs38, vs32 + WORD $0xf22a030e // xvmaddadp vs17, vs42, vs32 + WORD $0xf2460b0e // xvmaddadp vs18, vs38, vs33 + WORD $0xf26a0b0e // xvmaddadp vs19, vs42, vs33 + LXVD2X (R5)(R15), VS46 + WORD $0xf286130e // xvmaddadp vs20, vs38, vs34 + WORD $0xf2aa130e // xvmaddadp vs21, vs42, vs34 + WORD $0xf2c61b0e // xvmaddadp vs22, vs38, vs35 + WORD $0xf2ea1b0e // xvmaddadp vs23, vs42, vs35 + LXVD2X (R4)(R15), VS50 + WORD $0xf307030e // xvmaddadp vs24, vs39, vs32 + WORD $0xf32b030e // xvmaddadp vs25, vs43, vs32 + WORD $0xf3470b0e // xvmaddadp vs26, vs39, vs33 + WORD $0xf36b0b0e // xvmaddadp vs27, vs43, vs33 + LXVD2X (R5)(R16), VS47 + WORD $0xf387130e // xvmaddadp vs28, vs39, vs34 + WORD $0xf3ab130e // xvmaddadp vs29, vs43, vs34 + WORD $0xf3c71b0e // xvmaddadp vs30, vs39, vs35 + WORD $0xf3eb1b0e // xvmaddadp vs31, vs43, vs35 + LXVD2X (R4)(R16), VS51 + ADD $64, R4 + ADD $64, R5 + XXPERMDI VS48, VS48, $2, VS52 + XXPERMDI VS49, VS49, $2, VS53 + XXPERMDI VS50, VS50, $2, VS54 + XXPERMDI VS51, VS51, $2, VS55 + WORD $0xf010630e // xvmaddadp vs0, vs48, vs44 + WORD $0xf034630e // xvmaddadp vs1, vs52, vs44 + WORD $0xf0506b0e // xvmaddadp vs2, vs48, vs45 + WORD $0xf0746b0e // xvmaddadp vs3, vs52, vs45 + WORD $0xf090730e // xvmaddadp vs4, vs48, vs46 + WORD $0xf0b4730e // xvmaddadp vs5, vs52, vs46 + WORD $0xf0d07b0e // xvmaddadp vs6, vs48, vs47 + WORD $0xf0f47b0e // xvmaddadp vs7, vs52, vs47 + WORD $0xf111630e // xvmaddadp vs8, vs49, vs44 + WORD $0xf135630e // xvmaddadp vs9, vs53, vs44 + WORD $0xf1516b0e // xvmaddadp vs10, vs49, vs45 + WORD $0xf1756b0e // xvmaddadp vs11, vs53, vs45 + WORD $0xf191730e // xvmaddadp vs12, vs49, vs46 + WORD $0xf1b5730e // xvmaddadp vs13, vs53, vs46 + WORD $0xf1d17b0e // xvmaddadp vs14, vs49, vs47 + WORD $0xf1f57b0e // xvmaddadp vs15, vs53, vs47 + WORD $0xf212630e // xvmaddadp vs16, vs50, vs44 + WORD $0xf236630e // xvmaddadp vs17, vs54, vs44 + WORD $0xf2526b0e // xvmaddadp vs18, vs50, vs45 + WORD $0xf2766b0e // xvmaddadp vs19, vs54, vs45 + WORD $0xf292730e // xvmaddadp vs20, vs50, vs46 + WORD $0xf2b6730e // xvmaddadp vs21, vs54, vs46 + WORD $0xf2d27b0e // xvmaddadp vs22, vs50, vs47 + WORD $0xf2f67b0e // xvmaddadp vs23, vs54, vs47 + WORD $0xf313630e // xvmaddadp vs24, vs51, vs44 + WORD $0xf337630e // xvmaddadp vs25, vs55, vs44 + WORD $0xf3536b0e // xvmaddadp vs26, vs51, vs45 + WORD $0xf3776b0e // xvmaddadp vs27, vs55, vs45 + WORD $0xf393730e // xvmaddadp vs28, vs51, vs46 + WORD $0xf3b7730e // xvmaddadp vs29, vs55, vs46 + WORD $0xf3d37b0e // xvmaddadp vs30, vs51, vs47 + WORD $0xf3f77b0e // xvmaddadp vs31, vs55, vs47 + ADD $-2, R3 + BR gemmMicro8x8VSXloop2 +gemmMicro8x8VSXone: + CMP R3, $0 + BEQ gemmMicro8x8VSXstore + LXVD2X (R5)(R0), VS32 + LXVD2X (R4)(R0), VS36 + LXVD2X (R5)(R14), VS33 + LXVD2X (R4)(R14), VS37 + LXVD2X (R5)(R15), VS34 + LXVD2X (R4)(R15), VS38 + LXVD2X (R5)(R16), VS35 + LXVD2X (R4)(R16), VS39 + XXPERMDI VS36, VS36, $2, VS40 + XXPERMDI VS37, VS37, $2, VS41 + XXPERMDI VS38, VS38, $2, VS42 + XXPERMDI VS39, VS39, $2, VS43 + WORD $0xf004030e // xvmaddadp vs0, vs36, vs32 + WORD $0xf028030e // xvmaddadp vs1, vs40, vs32 + WORD $0xf0440b0e // xvmaddadp vs2, vs36, vs33 + WORD $0xf0680b0e // xvmaddadp vs3, vs40, vs33 + WORD $0xf084130e // xvmaddadp vs4, vs36, vs34 + WORD $0xf0a8130e // xvmaddadp vs5, vs40, vs34 + WORD $0xf0c41b0e // xvmaddadp vs6, vs36, vs35 + WORD $0xf0e81b0e // xvmaddadp vs7, vs40, vs35 + WORD $0xf105030e // xvmaddadp vs8, vs37, vs32 + WORD $0xf129030e // xvmaddadp vs9, vs41, vs32 + WORD $0xf1450b0e // xvmaddadp vs10, vs37, vs33 + WORD $0xf1690b0e // xvmaddadp vs11, vs41, vs33 + WORD $0xf185130e // xvmaddadp vs12, vs37, vs34 + WORD $0xf1a9130e // xvmaddadp vs13, vs41, vs34 + WORD $0xf1c51b0e // xvmaddadp vs14, vs37, vs35 + WORD $0xf1e91b0e // xvmaddadp vs15, vs41, vs35 + WORD $0xf206030e // xvmaddadp vs16, vs38, vs32 + WORD $0xf22a030e // xvmaddadp vs17, vs42, vs32 + WORD $0xf2460b0e // xvmaddadp vs18, vs38, vs33 + WORD $0xf26a0b0e // xvmaddadp vs19, vs42, vs33 + WORD $0xf286130e // xvmaddadp vs20, vs38, vs34 + WORD $0xf2aa130e // xvmaddadp vs21, vs42, vs34 + WORD $0xf2c61b0e // xvmaddadp vs22, vs38, vs35 + WORD $0xf2ea1b0e // xvmaddadp vs23, vs42, vs35 + WORD $0xf307030e // xvmaddadp vs24, vs39, vs32 + WORD $0xf32b030e // xvmaddadp vs25, vs43, vs32 + WORD $0xf3470b0e // xvmaddadp vs26, vs39, vs33 + WORD $0xf36b0b0e // xvmaddadp vs27, vs43, vs33 + WORD $0xf387130e // xvmaddadp vs28, vs39, vs34 + WORD $0xf3ab130e // xvmaddadp vs29, vs43, vs34 + WORD $0xf3c71b0e // xvmaddadp vs30, vs39, vs35 + WORD $0xf3eb1b0e // xvmaddadp vs31, vs43, vs35 +gemmMicro8x8VSXstore: + XXPERMDI VS0, VS1, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS2, VS3, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS4, VS5, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS6, VS7, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS1, VS0, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS3, VS2, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS5, VS4, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS7, VS6, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS8, VS9, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS10, VS11, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS12, VS13, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS14, VS15, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS9, VS8, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS11, VS10, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS13, VS12, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS15, VS14, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS16, VS17, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS18, VS19, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS20, VS21, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS22, VS23, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS17, VS16, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS19, VS18, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS21, VS20, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS23, VS22, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS24, VS25, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS26, VS27, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS28, VS29, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS30, VS31, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + ADD R7, R6 + XXPERMDI VS25, VS24, $1, VS56 + LXVD2X (R6)(R0), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R0) + XXPERMDI VS27, VS26, $1, VS56 + LXVD2X (R6)(R14), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R14) + XXPERMDI VS29, VS28, $1, VS56 + LXVD2X (R6)(R15), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R15) + XXPERMDI VS31, VS30, $1, VS56 + LXVD2X (R6)(R16), VS57 + WORD $0xf339c307 // xvadddp vs57, vs57, vs56 + STXVD2X VS57, (R6)(R16) + RET +