diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index b504e5b..6acf52d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -15,8 +15,9 @@ permissions:
 # 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), 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
+# -> 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
 # more than one, so no single run covers them all.
 #
 # The honest split:
@@ -37,9 +38,10 @@ 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/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.
+# 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.
 jobs:
   test:
     name: test (${{ matrix.os }})
@@ -124,12 +126,12 @@ jobs:
           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'
+        if: matrix.goarch == 'ppc64le' || matrix.goarch == 'loong64'
         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; }
+          ( cd internal/kernels/asmgen/${{ matrix.goarch }} && go run . )
+          mv internal/kernels/asmgen/${{ matrix.goarch }}/sum_${{ matrix.goarch }}.s internal/kernels/sum_${{ matrix.goarch }}.s
+          git diff --exit-code -- internal/kernels/sum_${{ matrix.goarch }}.s || {
+            echo "::error::generated ${{ matrix.goarch }} 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 4853170..5f4794b 100644
--- a/README.md
+++ b/README.md
@@ -59,14 +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 and ppc64le** the hot loops are
+`GOMAXPROCS`). On **amd64, arm64, ppc64le and loong64** 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), 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
+arm64, VSX on ppc64le, LASX on loong64 when the CPU has it), 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; 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
 (≤ 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 9be5b01..7e51ce0 100644
--- a/docs/perf.md
+++ b/docs/perf.md
@@ -603,6 +603,37 @@ 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.
 
+## loong64: LASX kernels (2026-10-04)
+
+Go's loong64 assembler has LASX vector float64 add/sub/mul/div, square root and
+256-bit loads, but no fused multiply-add and no broadcast load; go-asmgen
+v0.14.0 encodes `xvfmadd.d` and `xvldrepl.d` as `WORD`s (pinned against GNU as
+2.43, bit-identical to `math.FMA` on a Loongson 3C5000L; transitional until
+cmd/asm names them). go-ndarray uses them for sum, dot, sqrt, add/sub/mul/div
+and an 8×8 GEMM micro-kernel (16 four-lane accumulators, A broadcast by
+`xvldrepl.d`). LASX is not on every LoongArch CPU, so the kernels run only when
+the kernel reports it in AT_HWCAP (`/proc/self/auxv`), with the scalar code as
+the fallback; the parsing is a pure function tested on every target.
+
+Loongson 3C5000L (cfarm401, 32 cores), against the scalar code, best of 3
+interleaved runs (the host's load average of ~146 is processes stuck in
+uninterruptible sleep; its CPUs measured 100% idle):
+
+| op | before | after | gain |
+|----|--:|--:|:--:|
+| micro-kernel, 1 core | 7.0 GFLOP/s | 27.5 GFLOP/s | ×3.95 |
+| serial GEMM 512², 1 core | 5.3 GFLOP/s | 20.0 GFLOP/s | ×3.75 |
+| MatMul 64² / 512² / 1024², 32 threads | | | ×1.9 / ×1.3 / ×1.1 |
+| Sum 256 Ki / 4 Mi, 1 core | 658 / 11 750 µs | 65 / 4 615 µs | ×10.2 / ×2.5 |
+| Sum 1 Ki / 16 Ki (below the parallel threshold: one core) | | | ×16 / ×13 |
+| Sqrt 1 Ki–16 Ki | | | ×3.8 |
+| MatVec 1024² | 364 µs | 224 µs | ×1.6 |
+
+At 4 Mi elements in parallel, sum and sqrt are memory-bound (parity); the
+parallel 4 Mi sum is bimodal on this host (1.6 ms or 7 ms in either version),
+so only the minima are compared. No NumPy or BLAS is installable there (the
+host has no route to package mirrors), so there is no external reference.
+
 ## SIMD coverage
 
 - **amd64 (SSE2)** ships hand-vectorized `sum` (4-accumulator `ADDPD`), `sqrt`
@@ -629,7 +660,10 @@ baseline is POWER8.
 - **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
+- **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
   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
@@ -638,10 +672,12 @@ baseline is POWER8.
   **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;
+  `vfadd.d`, `VMULD` to `vfmul.d`, `XVADDD` to `xvfadd.d`) but no vector FMA
+  or broadcast load, 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 loong64, s390x and riscv64 kernels are work not yet done, not a
+  check. So 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
diff --git a/docs/plan-ndarray.md b/docs/plan-ndarray.md
index 0b3b4f2..2b4b312 100644
--- a/docs/plan-ndarray.md
+++ b/docs/plan-ndarray.md
@@ -123,7 +123,8 @@ suite (plus the per-arch CI jobs already wired in `.github/workflows/ci.yml`).
     lane-parallel grouping is a valid reordering — the same trade-off as
     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. The other targets (riscv64, loong64, s390x, 32-bit) keep the
+    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
     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
diff --git a/internal/kernels/asmgen/go.mod b/internal/kernels/asmgen/go.mod
index 156d43c..e2dc623 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.13.0
+require github.com/go-asmgen/asmgen v0.14.0
diff --git a/internal/kernels/asmgen/go.sum b/internal/kernels/asmgen/go.sum
index 1f76dfc..0a54b2f 100644
--- a/internal/kernels/asmgen/go.sum
+++ b/internal/kernels/asmgen/go.sum
@@ -1,2 +1,2 @@
-github.com/go-asmgen/asmgen v0.13.0 h1:O3DIZtgTgwHJ2uof6NhMCO1AnTZqub24r3DS83IjOb8=
-github.com/go-asmgen/asmgen v0.13.0/go.mod h1:Hj29fXioLRe9oSnQ18SuDBiOcCbu8Qwxba6/5TMr4+Y=
+github.com/go-asmgen/asmgen v0.14.0 h1:3f1L0X0RK1Vf1J9nb0Pei87vGRT5yhfNLNlecZ/SNFc=
+github.com/go-asmgen/asmgen v0.14.0/go.mod h1:Hj29fXioLRe9oSnQ18SuDBiOcCbu8Qwxba6/5TMr4+Y=
diff --git a/internal/kernels/asmgen/loong64/gen.go b/internal/kernels/asmgen/loong64/gen.go
new file mode 100644
index 0000000..c2c3611
--- /dev/null
+++ b/internal/kernels/asmgen/loong64/gen.go
@@ -0,0 +1,207 @@
+// Command gen produces sum_loong64.s, the LASX float64 kernels for loong64, via
+// go-asmgen. Run with: go run . (or `go generate` from the kernels package).
+//
+// Go's loong64 assembler has the LASX vector float64 add/sub/mul/div
+// (XVADDD, XVSUBD, XVMULD, XVDIVD), the square root (XVFSQRTD) and the
+// 256-bit loads and stores (XVMOVQ, with immediate offsets), but no fused
+// multiply-add and no broadcast load; those come from go-asmgen's encoders
+// (Builder.XVFMADDD, XVLDREPLD), which emit WORDs. LASX is not on every
+// LoongArch CPU, so the Go side calls these kernels only when the kernel
+// reports it (AT_HWCAP).
+//
+// Registers: R4.. hold the arguments; R12..R15 are scratch; X0..X31 are the
+// 256-bit LASX registers (four float64 lanes), F0..F31 their low lanes.
+package main
+
+import (
+	"fmt"
+	"os"
+
+	"github.com/go-asmgen/asmgen/emit"
+	"github.com/go-asmgen/asmgen/loong64"
+)
+
+const (
+	gemmMR = 8
+	gemmNR = 8
+)
+
+func main() {
+	f := emit.NewFile("loong64")
+	f.Add(sumKernel())
+	f.Add(dotKernel())
+	f.Add(sqrtKernel())
+	for _, k := range []struct{ name, vec, sca string }{
+		{"addLASX", "XVADDD", "ADDD"}, {"subLASX", "XVSUBD", "SUBD"},
+		{"mulLASX", "XVMULD", "MULD"}, {"divLASX", "XVDIVD", "DIVD"},
+	} {
+		f.Add(binKernel(k.name, k.vec, k.sca))
+	}
+	f.Add(gemmKernel())
+	if err := os.WriteFile("sum_loong64.s", []byte(f.String()), 0o644); err != nil {
+		fmt.Fprintln(os.Stderr, err)
+		os.Exit(1)
+	}
+	fmt.Println("wrote sum_loong64.s")
+}
+
+// fold adds X0..X3 into X0, then its four lanes into F0.
+func fold(b *loong64.Builder) {
+	b.Raw("XVADDD X1, X0, X0").Raw("XVADDD X3, X2, X2").Raw("XVADDD X2, X0, X0")
+	for l := 0; l < 4; l++ {
+		b.Raw("XVMOVQ X0.V[%d], R%d", l, 12+l)
+	}
+	b.Raw("MOVV R12, F0").Raw("MOVV R13, F1").Raw("MOVV R14, F2").Raw("MOVV R15, F3")
+	b.Raw("ADDD F1, F0").Raw("ADDD F3, F2").Raw("ADDD F2, F0")
+}
+
+func zero(b *loong64.Builder, n int) {
+	for r := 0; r < n; r++ {
+		b.Raw("XVXORV X%d, X%d, X%d", r, r, r)
+	}
+}
+
+// sumKernel builds sumLASX(a *float64, n int) float64: four LASX accumulators
+// (16 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 {
+	sig := loong64.Layout([]string{"a", "n"}, []loong64.Type{loong64.Ptr, loong64.Int64},
+		[]string{"ret"}, []loong64.Type{loong64.Float64})
+	b := loong64.NewFunc("sumLASX", sig, 0)
+	b.LoadArg("a", "R4").LoadArg("n", "R5")
+	zero(b, 4)
+	b.Label("sloop").Raw("MOVV $16, R12").Raw("BLT R5, R12, sfold")
+	for k := 0; k < 4; k++ {
+		b.Raw("XVMOVQ %d(R4), X%d", 32*k, 4+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("XVADDD X%d, X%d, X%d", 4+k, k, k)
+	}
+	b.Raw("ADDV $128, R4").Raw("ADDV $-16, R5").Raw("JMP sloop")
+	b.Label("sfold")
+	fold(b)
+	b.Label("stail").Raw("BEQ R5, sdone").
+		Raw("MOVD 0(R4), F1").Raw("ADDD F1, F0").
+		Raw("ADDV $8, R4").Raw("ADDV $-1, R5").Raw("JMP stail")
+	b.Label("sdone").StoreRet("F0", "ret").Ret()
+	return b.Func()
+}
+
+// dotKernel builds dotLASX(a, b *float64, n int) float64: four fused
+// multiply-add accumulators (xvfmadd.d), folded, plus a fused scalar tail.
+func dotKernel() *emit.Function {
+	sig := loong64.Layout([]string{"a", "b", "n"}, []loong64.Type{loong64.Ptr, loong64.Ptr, loong64.Int64},
+		[]string{"ret"}, []loong64.Type{loong64.Float64})
+	b := loong64.NewFunc("dotLASX", sig, 0)
+	b.LoadArg("a", "R4").LoadArg("b", "R6").LoadArg("n", "R5")
+	zero(b, 4)
+	b.Label("dloop").Raw("MOVV $16, R12").Raw("BLT R5, R12, dfold")
+	for k := 0; k < 4; k++ {
+		b.Raw("XVMOVQ %d(R4), X%d", 32*k, 4+k)
+		b.Raw("XVMOVQ %d(R6), X%d", 32*k, 8+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.XVFMADDD(k, 4+k, 8+k, k) // X(k) = a*b + X(k)
+	}
+	b.Raw("ADDV $128, R4").Raw("ADDV $128, R6").Raw("ADDV $-16, R5").Raw("JMP dloop")
+	b.Label("dfold")
+	fold(b)
+	b.Label("dtail").Raw("BEQ R5, ddone").
+		Raw("MOVD 0(R4), F1").Raw("MOVD 0(R6), F2").Raw("FMADDD F0, F2, F1, F0"). // F0 = F1*F2 + F0; Go's order is Fa, Fk, Fj, Fd
+		Raw("ADDV $8, R4").Raw("ADDV $8, R6").Raw("ADDV $-1, R5").Raw("JMP dtail")
+	b.Label("ddone").StoreRet("F0", "ret").Ret()
+	return b.Func()
+}
+
+// sqrtKernel builds sqrtLASX(dst, src *float64, n int): XVFSQRTD over 16
+// doubles per iteration, SQRTD 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 := loong64.Layout([]string{"dst", "src", "n"}, []loong64.Type{loong64.Ptr, loong64.Ptr, loong64.Int64}, nil, nil)
+	b := loong64.NewFunc("sqrtLASX", sig, 0)
+	b.LoadArg("dst", "R4").LoadArg("src", "R6").LoadArg("n", "R5")
+	b.Label("qloop").Raw("MOVV $16, R12").Raw("BLT R5, R12, qtail")
+	for k := 0; k < 4; k++ {
+		b.Raw("XVMOVQ %d(R6), X%d", 32*k, k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("XVFSQRTD X%d, X%d", k, 4+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("XVMOVQ X%d, %d(R4)", 4+k, 32*k)
+	}
+	b.Raw("ADDV $128, R4").Raw("ADDV $128, R6").Raw("ADDV $-16, R5").Raw("JMP qloop")
+	b.Label("qtail").Raw("BEQ R5, qdone").
+		Raw("MOVD 0(R6), F1").Raw("SQRTD F1, F1").Raw("MOVD F1, 0(R4)").
+		Raw("ADDV $8, R4").Raw("ADDV $8, R6").Raw("ADDV $-1, R5").Raw("JMP qtail")
+	b.Label("qdone").Ret()
+	return b.Func()
+}
+
+// binKernel builds name(dst, a, b *float64, n int): dst[i] = a[i] OP b[i], 16
+// doubles per iteration with the LASX op, a scalar tail with sca. Each lane is
+// one correctly rounded IEEE operation, bit-identical to the scalar oracle.
+func binKernel(name, vec, sca string) *emit.Function {
+	sig := loong64.Layout([]string{"dst", "a", "b", "n"},
+		[]loong64.Type{loong64.Ptr, loong64.Ptr, loong64.Ptr, loong64.Int64}, nil, nil)
+	b := loong64.NewFunc(name, sig, 0)
+	b.LoadArg("dst", "R4").LoadArg("a", "R6").LoadArg("b", "R7").LoadArg("n", "R5")
+	l := name[:3]
+	b.Label(l+"loop").Raw("MOVV $16, R12").Raw("BLT R5, R12, %stail", l)
+	for k := 0; k < 4; k++ {
+		b.Raw("XVMOVQ %d(R6), X%d", 32*k, k)
+		b.Raw("XVMOVQ %d(R7), X%d", 32*k, 4+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("%s X%d, X%d, X%d", vec, 4+k, k, 8+k) // X(8+k) = X(k) OP X(4+k)
+	}
+	for k := 0; k < 4; k++ {
+		b.Raw("XVMOVQ X%d, %d(R4)", 8+k, 32*k)
+	}
+	b.Raw("ADDV $128, R4").Raw("ADDV $128, R6").Raw("ADDV $128, R7").Raw("ADDV $-16, R5").Raw("JMP %sloop", l)
+	b.Label(l+"tail").Raw("BEQ R5, %sdone", l).
+		Raw("MOVD 0(R6), F1").Raw("MOVD 0(R7), F2").Raw("%s F2, F1", sca).Raw("MOVD F1, 0(R4)").
+		Raw("ADDV $8, R4").Raw("ADDV $8, R6").Raw("ADDV $8, R7").Raw("ADDV $-1, R5").Raw("JMP %stail", l)
+	b.Label(l + "done").Ret()
+	return b.Func()
+}
+
+// gemmKernel builds gemmMicro8x8LASX(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)
+//
+// 16 accumulators X0..X15 (row r in X(2r), X(2r+1), four columns each), the B
+// row in X16, X17, and each A value broadcast by xvldrepl.d into X18..X25,
+// then xvfmadd.d: 16 fused multiply-adds per k step against 10 loads. Fused,
+// like the other FMA kernels; on integer-valued data it equals the scalar
+// ikj oracle exactly.
+func gemmKernel() *emit.Function {
+	sig := loong64.Layout([]string{"kc", "pa", "pb", "c", "ldc"},
+		[]loong64.Type{loong64.Int64, loong64.Ptr, loong64.Ptr, loong64.Ptr, loong64.Int64}, nil, nil)
+	b := loong64.NewFunc("gemmMicro8x8LASX", sig, 0)
+	b.LoadArg("kc", "R5").LoadArg("pa", "R6").LoadArg("pb", "R7").LoadArg("c", "R8").LoadArg("ldc", "R9")
+	b.Raw("SLLV $3, R9") // ldc -> bytes
+	zero(b, 16)
+	b.Raw("BEQ R5, gstore")
+	b.Label("gloop")
+	b.Raw("XVMOVQ 0(R7), X16").Raw("XVMOVQ 32(R7), X17")
+	for r := 0; r < gemmMR; r++ {
+		b.XVLDREPLD(18+r, 6, 8*r)
+	}
+	for r := 0; r < gemmMR; r++ {
+		b.XVFMADDD(2*r, 18+r, 16, 2*r)
+		b.XVFMADDD(2*r+1, 18+r, 17, 2*r+1)
+	}
+	b.Raw("ADDV $%d, R6", 8*gemmMR).Raw("ADDV $%d, R7", 8*gemmNR).Raw("ADDV $-1, R5").Raw("BNE R5, gloop")
+	b.Label("gstore")
+	for r := 0; r < gemmMR; r++ {
+		b.Raw("XVMOVQ 0(R8), X16").Raw("XVMOVQ 32(R8), X17")
+		b.Raw("XVADDD X%d, X16, X16", 2*r).Raw("XVADDD X%d, X17, X17", 2*r+1)
+		b.Raw("XVMOVQ X16, 0(R8)").Raw("XVMOVQ X17, 32(R8)")
+		if r < gemmMR-1 {
+			b.Raw("ADDV R9, R8")
+		}
+	}
+	b.Ret()
+	return b.Func()
+}
diff --git a/internal/kernels/gemm_generic.go b/internal/kernels/gemm_generic.go
index e91b246..87fdbef 100644
--- a/internal/kernels/gemm_generic.go
+++ b/internal/kernels/gemm_generic.go
@@ -1,11 +1,11 @@
-//go:build !arm64 && !amd64 && !ppc64le
+//go:build !arm64 && !amd64 && !ppc64le && !loong64
 
 package kernels
 
 // 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
+// (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
+// 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
 // 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
diff --git a/internal/kernels/gemm_loong64.go b/internal/kernels/gemm_loong64.go
new file mode 100644
index 0000000..592c3c0
--- /dev/null
+++ b/internal/kernels/gemm_loong64.go
@@ -0,0 +1,56 @@
+package kernels
+
+// loong64 GEMM micro-kernel: an 8x8 LASX tile (gemmMicro8x8LASX in
+// sum_loong64.s, generated by go-asmgen): 16 four-lane accumulators, the B row
+// in two vectors, each A value broadcast by xvldrepl.d and a fused xvfmadd.d.
+// Without LASX, gemmMicroGo computes the same tile in scalar code.
+
+// MR, NR are the micro-kernel's register-block dimensions.
+const (
+	MR = 8
+	NR = 8
+)
+
+func gemmMicro(kc int, pa, pb, dst []float64, ldc int) {
+	if useLASX {
+		gemmMicro8x8LASX(kc, &pa[0], &pb[0], &dst[0], ldc)
+		return
+	}
+	gemmMicroGo(kc, pa, pb, dst, ldc)
+}
+
+// gemmMicroGo is the scalar MR x NR tile: dst[r*ldc+c] += sum_p
+// pa[p*MR+r]*pb[p*NR+c], in p order per element.
+func gemmMicroGo(kc int, pa, pb, dst []float64, ldc int) {
+	var acc [MR * NR]float64
+	for p := 0; p < kc; p++ {
+		a := pa[p*MR : p*MR+MR : p*MR+MR]
+		b := pb[p*NR : p*NR+NR : p*NR+NR]
+		for r, av := range a {
+			row := acc[r*NR : r*NR+NR : r*NR+NR]
+			for c, bv := range b {
+				row[c] += av * bv
+			}
+		}
+	}
+	for r := 0; r < MR; r++ {
+		d := dst[r*ldc : r*ldc+NR : r*ldc+NR]
+		for c := range d {
+			d[c] += acc[r*NR+c]
+		}
+	}
+}
+
+// dotSIMD is the inner product (LASX kernel, four fused accumulators).
+func dotSIMD(a, b []float64) float64 {
+	if !useLASX || len(a) == 0 {
+		return dotRange(a, b)
+	}
+	return dotLASX(&a[0], &b[0], len(a))
+}
+
+//go:noescape
+func gemmMicro8x8LASX(kc int, pa, pb, c *float64, ldc int)
+
+//go:noescape
+func dotLASX(a, b *float64, n int) float64
diff --git a/internal/kernels/gemm_loong64_test.go b/internal/kernels/gemm_loong64_test.go
new file mode 100644
index 0000000..38623c7
--- /dev/null
+++ b/internal/kernels/gemm_loong64_test.go
@@ -0,0 +1,19 @@
+package kernels
+
+import "testing"
+
+// TestScalarFallbackLoong64 runs the GEMM, dot and elementwise suites with
+// LASX forced off, so the path a CPU without LASX takes is exercised on one
+// that has it.
+func TestScalarFallbackLoong64(t *testing.T) {
+	saved := useLASX
+	useLASX = false
+	defer func() { useLASX = saved }()
+	t.Run("gemm", TestGemmCorrectness)
+	t.Run("dot", TestDotSIMD)
+	t.Run("bin", TestBinSIMD)
+}
+
+func TestLASXDetected(t *testing.T) {
+	t.Logf("useLASX = %v (AT_HWCAP)", useLASX)
+}
diff --git a/internal/kernels/hwcap.go b/internal/kernels/hwcap.go
new file mode 100644
index 0000000..ca64163
--- /dev/null
+++ b/internal/kernels/hwcap.go
@@ -0,0 +1,25 @@
+package kernels
+
+import "encoding/binary"
+
+// atHWCAP is the auxiliary-vector key of the hardware-capability bits the
+// kernel exports (AT_HWCAP).
+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
+// 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.
+func hwcapFrom(auxv []byte) uint64 {
+	for len(auxv) >= 16 {
+		k, v := binary.LittleEndian.Uint64(auxv), binary.LittleEndian.Uint64(auxv[8:])
+		if k == 0 {
+			break
+		}
+		if k == atHWCAP {
+			return v
+		}
+		auxv = auxv[16:]
+	}
+	return 0
+}
diff --git a/internal/kernels/hwcap_test.go b/internal/kernels/hwcap_test.go
new file mode 100644
index 0000000..b3c202f
--- /dev/null
+++ b/internal/kernels/hwcap_test.go
@@ -0,0 +1,32 @@
+package kernels
+
+import (
+	"encoding/binary"
+	"testing"
+)
+
+func auxvOf(pairs ...uint64) []byte {
+	b := make([]byte, 8*len(pairs))
+	for i, p := range pairs {
+		binary.LittleEndian.PutUint64(b[8*i:], p)
+	}
+	return b
+}
+
+func TestHWCAPFrom(t *testing.T) {
+	for _, c := range []struct {
+		name string
+		auxv []byte
+		want uint64
+	}{
+		{"found after others", auxvOf(6, 4096, atHWCAP, 0x3f, 0, 0), 0x3f},
+		{"absent", auxvOf(6, 4096, 0, 0), 0},
+		{"after the terminator", auxvOf(0, 0, atHWCAP, 1), 0},
+		{"empty", nil, 0},
+		{"short tail", append(auxvOf(6, 1), 1, 2, 3), 0},
+	} {
+		if got := hwcapFrom(c.auxv); got != c.want {
+			t.Errorf("%s: %#x, want %#x", c.name, got, c.want)
+		}
+	}
+}
diff --git a/internal/kernels/kernels.go b/internal/kernels/kernels.go
index 984ae6b..91e3c1e 100644
--- a/internal/kernels/kernels.go
+++ b/internal/kernels/kernels.go
@@ -3,7 +3,7 @@
 // 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 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
+// on amd64, arm64, ppc64le and loong64 (the *_<arch> files), while these
 // scalar versions remain the reference and the fallback elsewhere.
 //
 // Kernels operate on flat, contiguous []float64 slices. Shape, stride and
diff --git a/internal/kernels/reduce_generic.go b/internal/kernels/reduce_generic.go
index f204162..d3d3b55 100644
--- a/internal/kernels/reduce_generic.go
+++ b/internal/kernels/reduce_generic.go
@@ -1,4 +1,4 @@
-//go:build !arm64 && !amd64 && !ppc64le
+//go:build !arm64 && !amd64 && !ppc64le && !loong64
 
 package kernels
 
@@ -7,20 +7,17 @@ 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) 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
+// 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
-// 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.
+// 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) 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
+// 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
 // run states which path it validated.
 const HaveReduceSIMD = false
 
diff --git a/internal/kernels/reduce_loong64.go b/internal/kernels/reduce_loong64.go
new file mode 100644
index 0000000..5fce39d
--- /dev/null
+++ b/internal/kernels/reduce_loong64.go
@@ -0,0 +1,95 @@
+package kernels
+
+import "os"
+
+// loong64 LASX kernels generated by go-asmgen (sum_loong64.s). LASX (256-bit,
+// four float64 lanes) is not on every LoongArch CPU, so each wrapper uses it
+// only when the kernel reports it in AT_HWCAP, and runs the scalar oracle
+// otherwise.
+//
+//go:generate sh -c "cd asmgen/loong64 && go run . && mv sum_loong64.s ../../sum_loong64.s"
+
+// HaveReduceSIMD reports that this build ships SIMD kernels; they run when
+// useLASX is true.
+const HaveReduceSIMD = true
+
+// hwcapLASX is HWCAP_LOONGARCH_LASX (Linux arch/loongarch/include/uapi/asm/hwcap.h).
+const hwcapLASX = 1 << 5
+
+// useLASX selects the LASX kernels. A var so tests can force the scalar path
+// on hardware that has LASX (never the reverse: that would fault).
+var useLASX = lasxFromAuxv()
+
+func lasxFromAuxv() bool {
+	auxv, err := os.ReadFile("/proc/self/auxv")
+	return err == nil && hwcapFrom(auxv)&hwcapLASX != 0
+}
+
+func sumSIMD(a []float64) float64 {
+	if !useLASX || len(a) == 0 {
+		return Sum(a)
+	}
+	return sumLASX(&a[0], len(a))
+}
+
+func sqrtSIMD(dst, src []float64) {
+	if !useLASX || len(dst) == 0 {
+		sqrtScalar(dst, src)
+		return
+	}
+	sqrtLASX(&dst[0], &src[0], len(dst))
+}
+
+// max/min stay on the four-accumulator scalar reducers.
+func maxSIMD(a []float64) float64 { return maxUnrolled(a) }
+func minSIMD(a []float64) float64 { return minUnrolled(a) }
+
+func addBin(dst, a, b []float64) {
+	if !useLASX || len(dst) == 0 {
+		Add(dst, a, b)
+		return
+	}
+	addLASX(&dst[0], &a[0], &b[0], len(dst))
+}
+
+func subBin(dst, a, b []float64) {
+	if !useLASX || len(dst) == 0 {
+		Sub(dst, a, b)
+		return
+	}
+	subLASX(&dst[0], &a[0], &b[0], len(dst))
+}
+
+func mulBin(dst, a, b []float64) {
+	if !useLASX || len(dst) == 0 {
+		Mul(dst, a, b)
+		return
+	}
+	mulLASX(&dst[0], &a[0], &b[0], len(dst))
+}
+
+func divBin(dst, a, b []float64) {
+	if !useLASX || len(dst) == 0 {
+		Div(dst, a, b)
+		return
+	}
+	divLASX(&dst[0], &a[0], &b[0], len(dst))
+}
+
+//go:noescape
+func sumLASX(a *float64, n int) float64
+
+//go:noescape
+func sqrtLASX(dst, src *float64, n int)
+
+//go:noescape
+func addLASX(dst, a, b *float64, n int)
+
+//go:noescape
+func subLASX(dst, a, b *float64, n int)
+
+//go:noescape
+func mulLASX(dst, a, b *float64, n int)
+
+//go:noescape
+func divLASX(dst, a, b *float64, n int)
diff --git a/internal/kernels/sum_loong64.s b/internal/kernels/sum_loong64.s
new file mode 100644
index 0000000..55d17ba
--- /dev/null
+++ b/internal/kernels/sum_loong64.s
@@ -0,0 +1,428 @@
+// Code generated by go-asmgen. DO NOT EDIT.
+
+//go:build loong64
+
+#include "textflag.h"
+
+TEXT ·sumLASX(SB), NOSPLIT, $0-24
+	MOVV a+0(FP), R4
+	MOVV n+8(FP), R5
+	XVXORV X0, X0, X0
+	XVXORV X1, X1, X1
+	XVXORV X2, X2, X2
+	XVXORV X3, X3, X3
+sloop:
+	MOVV $16, R12
+	BLT R5, R12, sfold
+	XVMOVQ 0(R4), X4
+	XVMOVQ 32(R4), X5
+	XVMOVQ 64(R4), X6
+	XVMOVQ 96(R4), X7
+	XVADDD X4, X0, X0
+	XVADDD X5, X1, X1
+	XVADDD X6, X2, X2
+	XVADDD X7, X3, X3
+	ADDV $128, R4
+	ADDV $-16, R5
+	JMP sloop
+sfold:
+	XVADDD X1, X0, X0
+	XVADDD X3, X2, X2
+	XVADDD X2, X0, X0
+	XVMOVQ X0.V[0], R12
+	XVMOVQ X0.V[1], R13
+	XVMOVQ X0.V[2], R14
+	XVMOVQ X0.V[3], R15
+	MOVV R12, F0
+	MOVV R13, F1
+	MOVV R14, F2
+	MOVV R15, F3
+	ADDD F1, F0
+	ADDD F3, F2
+	ADDD F2, F0
+stail:
+	BEQ R5, sdone
+	MOVD 0(R4), F1
+	ADDD F1, F0
+	ADDV $8, R4
+	ADDV $-1, R5
+	JMP stail
+sdone:
+	MOVD F0, ret+16(FP)
+	RET
+
+TEXT ·dotLASX(SB), NOSPLIT, $0-32
+	MOVV a+0(FP), R4
+	MOVV b+8(FP), R6
+	MOVV n+16(FP), R5
+	XVXORV X0, X0, X0
+	XVXORV X1, X1, X1
+	XVXORV X2, X2, X2
+	XVXORV X3, X3, X3
+dloop:
+	MOVV $16, R12
+	BLT R5, R12, dfold
+	XVMOVQ 0(R4), X4
+	XVMOVQ 0(R6), X8
+	XVMOVQ 32(R4), X5
+	XVMOVQ 32(R6), X9
+	XVMOVQ 64(R4), X6
+	XVMOVQ 64(R6), X10
+	XVMOVQ 96(R4), X7
+	XVMOVQ 96(R6), X11
+	WORD $0x0a202080 // xvfmadd.d $xr0, $xr4, $xr8, $xr0
+	WORD $0x0a20a4a1 // xvfmadd.d $xr1, $xr5, $xr9, $xr1
+	WORD $0x0a2128c2 // xvfmadd.d $xr2, $xr6, $xr10, $xr2
+	WORD $0x0a21ace3 // xvfmadd.d $xr3, $xr7, $xr11, $xr3
+	ADDV $128, R4
+	ADDV $128, R6
+	ADDV $-16, R5
+	JMP dloop
+dfold:
+	XVADDD X1, X0, X0
+	XVADDD X3, X2, X2
+	XVADDD X2, X0, X0
+	XVMOVQ X0.V[0], R12
+	XVMOVQ X0.V[1], R13
+	XVMOVQ X0.V[2], R14
+	XVMOVQ X0.V[3], R15
+	MOVV R12, F0
+	MOVV R13, F1
+	MOVV R14, F2
+	MOVV R15, F3
+	ADDD F1, F0
+	ADDD F3, F2
+	ADDD F2, F0
+dtail:
+	BEQ R5, ddone
+	MOVD 0(R4), F1
+	MOVD 0(R6), F2
+	FMADDD F0, F2, F1, F0
+	ADDV $8, R4
+	ADDV $8, R6
+	ADDV $-1, R5
+	JMP dtail
+ddone:
+	MOVD F0, ret+24(FP)
+	RET
+
+TEXT ·sqrtLASX(SB), NOSPLIT, $0-24
+	MOVV dst+0(FP), R4
+	MOVV src+8(FP), R6
+	MOVV n+16(FP), R5
+qloop:
+	MOVV $16, R12
+	BLT R5, R12, qtail
+	XVMOVQ 0(R6), X0
+	XVMOVQ 32(R6), X1
+	XVMOVQ 64(R6), X2
+	XVMOVQ 96(R6), X3
+	XVFSQRTD X0, X4
+	XVFSQRTD X1, X5
+	XVFSQRTD X2, X6
+	XVFSQRTD X3, X7
+	XVMOVQ X4, 0(R4)
+	XVMOVQ X5, 32(R4)
+	XVMOVQ X6, 64(R4)
+	XVMOVQ X7, 96(R4)
+	ADDV $128, R4
+	ADDV $128, R6
+	ADDV $-16, R5
+	JMP qloop
+qtail:
+	BEQ R5, qdone
+	MOVD 0(R6), F1
+	SQRTD F1, F1
+	MOVD F1, 0(R4)
+	ADDV $8, R4
+	ADDV $8, R6
+	ADDV $-1, R5
+	JMP qtail
+qdone:
+	RET
+
+TEXT ·addLASX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+addloop:
+	MOVV $16, R12
+	BLT R5, R12, addtail
+	XVMOVQ 0(R6), X0
+	XVMOVQ 0(R7), X4
+	XVMOVQ 32(R6), X1
+	XVMOVQ 32(R7), X5
+	XVMOVQ 64(R6), X2
+	XVMOVQ 64(R7), X6
+	XVMOVQ 96(R6), X3
+	XVMOVQ 96(R7), X7
+	XVADDD X4, X0, X8
+	XVADDD X5, X1, X9
+	XVADDD X6, X2, X10
+	XVADDD X7, X3, X11
+	XVMOVQ X8, 0(R4)
+	XVMOVQ X9, 32(R4)
+	XVMOVQ X10, 64(R4)
+	XVMOVQ X11, 96(R4)
+	ADDV $128, R4
+	ADDV $128, R6
+	ADDV $128, R7
+	ADDV $-16, R5
+	JMP addloop
+addtail:
+	BEQ R5, adddone
+	MOVD 0(R6), F1
+	MOVD 0(R7), F2
+	ADDD F2, F1
+	MOVD F1, 0(R4)
+	ADDV $8, R4
+	ADDV $8, R6
+	ADDV $8, R7
+	ADDV $-1, R5
+	JMP addtail
+adddone:
+	RET
+
+TEXT ·subLASX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+subloop:
+	MOVV $16, R12
+	BLT R5, R12, subtail
+	XVMOVQ 0(R6), X0
+	XVMOVQ 0(R7), X4
+	XVMOVQ 32(R6), X1
+	XVMOVQ 32(R7), X5
+	XVMOVQ 64(R6), X2
+	XVMOVQ 64(R7), X6
+	XVMOVQ 96(R6), X3
+	XVMOVQ 96(R7), X7
+	XVSUBD X4, X0, X8
+	XVSUBD X5, X1, X9
+	XVSUBD X6, X2, X10
+	XVSUBD X7, X3, X11
+	XVMOVQ X8, 0(R4)
+	XVMOVQ X9, 32(R4)
+	XVMOVQ X10, 64(R4)
+	XVMOVQ X11, 96(R4)
+	ADDV $128, R4
+	ADDV $128, R6
+	ADDV $128, R7
+	ADDV $-16, R5
+	JMP subloop
+subtail:
+	BEQ R5, subdone
+	MOVD 0(R6), F1
+	MOVD 0(R7), F2
+	SUBD F2, F1
+	MOVD F1, 0(R4)
+	ADDV $8, R4
+	ADDV $8, R6
+	ADDV $8, R7
+	ADDV $-1, R5
+	JMP subtail
+subdone:
+	RET
+
+TEXT ·mulLASX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+mulloop:
+	MOVV $16, R12
+	BLT R5, R12, multail
+	XVMOVQ 0(R6), X0
+	XVMOVQ 0(R7), X4
+	XVMOVQ 32(R6), X1
+	XVMOVQ 32(R7), X5
+	XVMOVQ 64(R6), X2
+	XVMOVQ 64(R7), X6
+	XVMOVQ 96(R6), X3
+	XVMOVQ 96(R7), X7
+	XVMULD X4, X0, X8
+	XVMULD X5, X1, X9
+	XVMULD X6, X2, X10
+	XVMULD X7, X3, X11
+	XVMOVQ X8, 0(R4)
+	XVMOVQ X9, 32(R4)
+	XVMOVQ X10, 64(R4)
+	XVMOVQ X11, 96(R4)
+	ADDV $128, R4
+	ADDV $128, R6
+	ADDV $128, R7
+	ADDV $-16, R5
+	JMP mulloop
+multail:
+	BEQ R5, muldone
+	MOVD 0(R6), F1
+	MOVD 0(R7), F2
+	MULD F2, F1
+	MOVD F1, 0(R4)
+	ADDV $8, R4
+	ADDV $8, R6
+	ADDV $8, R7
+	ADDV $-1, R5
+	JMP multail
+muldone:
+	RET
+
+TEXT ·divLASX(SB), NOSPLIT, $0-32
+	MOVV dst+0(FP), R4
+	MOVV a+8(FP), R6
+	MOVV b+16(FP), R7
+	MOVV n+24(FP), R5
+divloop:
+	MOVV $16, R12
+	BLT R5, R12, divtail
+	XVMOVQ 0(R6), X0
+	XVMOVQ 0(R7), X4
+	XVMOVQ 32(R6), X1
+	XVMOVQ 32(R7), X5
+	XVMOVQ 64(R6), X2
+	XVMOVQ 64(R7), X6
+	XVMOVQ 96(R6), X3
+	XVMOVQ 96(R7), X7
+	XVDIVD X4, X0, X8
+	XVDIVD X5, X1, X9
+	XVDIVD X6, X2, X10
+	XVDIVD X7, X3, X11
+	XVMOVQ X8, 0(R4)
+	XVMOVQ X9, 32(R4)
+	XVMOVQ X10, 64(R4)
+	XVMOVQ X11, 96(R4)
+	ADDV $128, R4
+	ADDV $128, R6
+	ADDV $128, R7
+	ADDV $-16, R5
+	JMP divloop
+divtail:
+	BEQ R5, divdone
+	MOVD 0(R6), F1
+	MOVD 0(R7), F2
+	DIVD F2, F1
+	MOVD F1, 0(R4)
+	ADDV $8, R4
+	ADDV $8, R6
+	ADDV $8, R7
+	ADDV $-1, R5
+	JMP divtail
+divdone:
+	RET
+
+TEXT ·gemmMicro8x8LASX(SB), NOSPLIT, $0-40
+	MOVV kc+0(FP), R5
+	MOVV pa+8(FP), R6
+	MOVV pb+16(FP), R7
+	MOVV c+24(FP), R8
+	MOVV ldc+32(FP), R9
+	SLLV $3, R9
+	XVXORV X0, X0, X0
+	XVXORV X1, X1, X1
+	XVXORV X2, X2, X2
+	XVXORV X3, X3, X3
+	XVXORV X4, X4, X4
+	XVXORV X5, X5, X5
+	XVXORV X6, X6, X6
+	XVXORV X7, X7, X7
+	XVXORV X8, X8, X8
+	XVXORV X9, X9, X9
+	XVXORV X10, X10, X10
+	XVXORV X11, X11, X11
+	XVXORV X12, X12, X12
+	XVXORV X13, X13, X13
+	XVXORV X14, X14, X14
+	XVXORV X15, X15, X15
+	BEQ R5, gstore
+gloop:
+	XVMOVQ 0(R7), X16
+	XVMOVQ 32(R7), X17
+	WORD $0x321000d2 // xvldrepl.d $xr18, $r6, 0
+	WORD $0x321004d3 // xvldrepl.d $xr19, $r6, 8
+	WORD $0x321008d4 // xvldrepl.d $xr20, $r6, 16
+	WORD $0x32100cd5 // xvldrepl.d $xr21, $r6, 24
+	WORD $0x321010d6 // xvldrepl.d $xr22, $r6, 32
+	WORD $0x321014d7 // xvldrepl.d $xr23, $r6, 40
+	WORD $0x321018d8 // xvldrepl.d $xr24, $r6, 48
+	WORD $0x32101cd9 // xvldrepl.d $xr25, $r6, 56
+	WORD $0x0a204240 // xvfmadd.d $xr0, $xr18, $xr16, $xr0
+	WORD $0x0a20c641 // xvfmadd.d $xr1, $xr18, $xr17, $xr1
+	WORD $0x0a214262 // xvfmadd.d $xr2, $xr19, $xr16, $xr2
+	WORD $0x0a21c663 // xvfmadd.d $xr3, $xr19, $xr17, $xr3
+	WORD $0x0a224284 // xvfmadd.d $xr4, $xr20, $xr16, $xr4
+	WORD $0x0a22c685 // xvfmadd.d $xr5, $xr20, $xr17, $xr5
+	WORD $0x0a2342a6 // xvfmadd.d $xr6, $xr21, $xr16, $xr6
+	WORD $0x0a23c6a7 // xvfmadd.d $xr7, $xr21, $xr17, $xr7
+	WORD $0x0a2442c8 // xvfmadd.d $xr8, $xr22, $xr16, $xr8
+	WORD $0x0a24c6c9 // xvfmadd.d $xr9, $xr22, $xr17, $xr9
+	WORD $0x0a2542ea // xvfmadd.d $xr10, $xr23, $xr16, $xr10
+	WORD $0x0a25c6eb // xvfmadd.d $xr11, $xr23, $xr17, $xr11
+	WORD $0x0a26430c // xvfmadd.d $xr12, $xr24, $xr16, $xr12
+	WORD $0x0a26c70d // xvfmadd.d $xr13, $xr24, $xr17, $xr13
+	WORD $0x0a27432e // xvfmadd.d $xr14, $xr25, $xr16, $xr14
+	WORD $0x0a27c72f // xvfmadd.d $xr15, $xr25, $xr17, $xr15
+	ADDV $64, R6
+	ADDV $64, R7
+	ADDV $-1, R5
+	BNE R5, gloop
+gstore:
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X0, X16, X16
+	XVADDD X1, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X2, X16, X16
+	XVADDD X3, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X4, X16, X16
+	XVADDD X5, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X6, X16, X16
+	XVADDD X7, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X8, X16, X16
+	XVADDD X9, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X10, X16, X16
+	XVADDD X11, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X12, X16, X16
+	XVADDD X13, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	ADDV R9, R8
+	XVMOVQ 0(R8), X16
+	XVMOVQ 32(R8), X17
+	XVADDD X14, X16, X16
+	XVADDD X15, X17, X17
+	XVMOVQ X16, 0(R8)
+	XVMOVQ X17, 32(R8)
+	RET
+