perf(linalg): fewer allocations per factorization (first step of #104)

.gitignore

@@ -1,2 +1,5 @@
 /cover.out
 *.out
+# go test -c / -memprofile leave the package's test binary beside its sources;
+# a 9.5 MB linalg.test was committed with #46.
+*.test

docs/perf-table.md

@@ -2,7 +2,7 @@
 
 <!-- Generated by refbench from docs/perf.json. Do not edit. -->
 
-- Date: 2026-11-12T12:08:39+01:00
+- Date: 2026-11-12T21:55:30+01:00
 - Commit: 
 - Go: go1.27.4
 - Machine: darwin/arm64, 16 CPUs
@@ -65,8 +65,8 @@ ns/call, median [first quartile, third quartile]; ratio = Go / reference.
 | Array.DivScalar | `Array.DivScalar` | numpy x / 1.5 | — | — | — | deferred | none | 8 |
 | Array.SqrtInto | `Array.SqrtInto` | numpy np.sqrt (out= preallocated) | — | — | — | deferred | none | 1 |
 | Array.Clip | `Array.Clip` | numpy np.clip | — | — | — | deferred | none | 4 |
-| Array.Dot | `Array.Dot` | numpy np.dot, 128 x 128 (OpenBLAS, one thread) | — | — | — | deferred | none | 13 |
-| Array.MatMul | `Array.MatMul` | numpy np.matmul, 128 x 128 (OpenBLAS, one thread) | — | — | — | deferred | none | 13 |
+| Array.Dot | `Array.Dot` | numpy np.dot, 128 x 128 (OpenBLAS, one thread) | — | — | — | deferred | none | 8 |
+| Array.MatMul | `Array.MatMul` | numpy np.matmul, 128 x 128 (OpenBLAS, one thread) | — | — | — | deferred | none | 8 |
 | Array.Outer | `Array.Outer` | numpy np.outer, 1000 x 1000 | — | — | — | deferred | none | 4 |
 | Array.Transpose | `Array.Transpose` | numpy x.T (a view; Go may copy) | — | — | — | deferred | none | 3 |
 | Array.Reshape | `Array.Reshape` | numpy x.reshape (a view) | — | — | — | deferred | none | 5 |
@@ -156,27 +156,27 @@ ns/call, median [first quartile, third quartile]; ratio = Go / reference.
 | fft.IFFTShift | `fft.IFFTShift` | numpy np.fft.ifftshift | — | — | — | deferred | none | 29 |
 | fft.FFTFreq | `fft.FFTFreq` | numpy np.fft.fftfreq | — | — | — | deferred | none | 7 |
 | fft.RFFTFreq | `fft.RFFTFreq` | numpy np.fft.rfftfreq | — | — | — | deferred | none | 7 |
-| linalg.Det | `linalg.Det` | numpy np.linalg.det(a) | — | — | — | deferred | none | 25 |
-| linalg.Inv | `linalg.Inv` | numpy np.linalg.inv(a) | — | — | — | deferred | none | 30 |
+| linalg.Det | `linalg.Det` | numpy np.linalg.det(a) | — | — | — | deferred | none | 24 |
+| linalg.Inv | `linalg.Inv` | numpy np.linalg.inv(a) | — | — | — | deferred | none | 27 |
 | linalg.Eigvals | `linalg.Eigvals` | numpy np.linalg.eigvals(a) | — | — | — | deferred | none | 25 |
-| linalg.Pinv | `linalg.Pinv` | numpy np.linalg.pinv(a) | — | — | — | deferred | none | 255 |
+| linalg.Pinv | `linalg.Pinv` | numpy np.linalg.pinv(a) | — | — | — | deferred | none | 230 |
 | linalg.SVDVals | `linalg.SVDVals` | numpy np.linalg.svd(a, compute_uv=False) | — | — | — | deferred | none | 23 |
 | linalg.MatrixRank | `linalg.MatrixRank` | numpy np.linalg.matrix_rank(a) | — | — | — | deferred | none | 34 |
-| linalg.Cholesky | `linalg.Cholesky` | numpy np.linalg.cholesky(a) | — | — | — | deferred | none | 20 |
+| linalg.Cholesky | `linalg.Cholesky` | numpy np.linalg.cholesky(a) | — | — | — | deferred | none | 19 |
 | linalg.Eigvalsh | `linalg.Eigvalsh` | numpy np.linalg.eigvalsh(a) | — | — | — | deferred | none | 22 |
 | linalg.Cond | `linalg.Cond` | numpy np.linalg.cond(a) | — | — | — | deferred | none | 39 |
 | linalg.Norm | `linalg.Norm` | numpy np.linalg.norm(a) | — | — | — | deferred | none | 11 |
 | linalg.MatrixRankTol | `linalg.MatrixRankTol` | numpy np.linalg.matrix_rank(a, tol=1e-8) | — | — | — | deferred | none | 34 |
-| linalg.PinvRcond | `linalg.PinvRcond` | numpy np.linalg.pinv(a, rcond=1e-10) | — | — | — | deferred | none | 255 |
+| linalg.PinvRcond | `linalg.PinvRcond` | numpy np.linalg.pinv(a, rcond=1e-10) | — | — | — | deferred | none | 230 |
 | linalg.MatrixPower | `linalg.MatrixPower` | numpy np.linalg.matrix_power(a, 5) | — | — | — | deferred | none | 28 |
-| linalg.Eig | `linalg.Eig` | numpy np.linalg.eig(a) | — | — | — | deferred | none | 89 |
-| linalg.Eigh | `linalg.Eigh` | numpy np.linalg.eigh(a) | — | — | — | deferred | none | 166 |
-| linalg.QR | `linalg.QR` | numpy np.linalg.qr(a) | — | — | — | deferred | none | 93 |
-| linalg.SVD | `linalg.SVD` | numpy np.linalg.svd(a, full_matrices=False) | — | — | — | deferred | none | 265 |
-| linalg.Slogdet | `linalg.Slogdet` | numpy np.linalg.slogdet(a) | — | — | — | deferred | none | 29 |
-| linalg.Solve | `linalg.Solve` | numpy np.linalg.solve(a, b) | — | — | — | deferred | none | 40 |
-| linalg.Lstsq | `linalg.Lstsq` | numpy np.linalg.lstsq(a, b, rcond=None) | — | — | — | deferred | none | 276 |
-| linalg.LstsqRcond | `linalg.LstsqRcond` | numpy np.linalg.lstsq(a, b, rcond=1e-10) | — | — | — | deferred | none | 276 |
+| linalg.Eig | `linalg.Eig` | numpy np.linalg.eig(a) | — | — | — | deferred | none | 80 |
+| linalg.Eigh | `linalg.Eigh` | numpy np.linalg.eigh(a) | — | — | — | deferred | none | 154 |
+| linalg.QR | `linalg.QR` | numpy np.linalg.qr(a) | — | — | — | deferred | none | 81 |
+| linalg.SVD | `linalg.SVD` | numpy np.linalg.svd(a, full_matrices=False) | — | — | — | deferred | none | 241 |
+| linalg.Slogdet | `linalg.Slogdet` | numpy np.linalg.slogdet(a) | — | — | — | deferred | none | 28 |
+| linalg.Solve | `linalg.Solve` | numpy np.linalg.solve(a, b) | — | — | — | deferred | none | 37 |
+| linalg.Lstsq | `linalg.Lstsq` | numpy np.linalg.lstsq(a, b, rcond=None) | — | — | — | deferred | none | 250 |
+| linalg.LstsqRcond | `linalg.LstsqRcond` | numpy np.linalg.lstsq(a, b, rcond=1e-10) | — | — | — | deferred | none | 250 |
 | linalg.MultiDot | `linalg.MultiDot` | numpy np.linalg.multi_dot([a, b, c]) | — | — | — | deferred | none | 25 |
 | random.Generator.Beta | `random.Generator.Beta` | numpy Generator.beta(2.0, 3.0, 100000) | — | — | — | deferred | none | 6 |
 | random.Generator.Binomial | `random.Generator.Binomial` | numpy Generator.binomial(100, 0.3, 100000) | — | — | — | deferred | none | 6 |

docs/perf.json

@@ -1,7 +1,7 @@
 {
   "schema": 1,
   "status": "deferred",
-  "date": "2026-11-12T12:08:39+01:00",
+  "date": "2026-11-12T21:55:30+01:00",
   "commit": "",
   "repo": "go-ndarray/ndarray",
   "go_version": "go1.27.4",
@@ -973,7 +973,7 @@
       "func": "Array.Dot",
       "reference": "numpy np.dot, 128 x 128 (OpenBLAS, one thread)",
       "parity": "none",
-      "allocs": 13,
+      "allocs": 8,
       "go": {
         "median": 0,
         "q1": 0,
@@ -991,7 +991,7 @@
       "func": "Array.MatMul",
       "reference": "numpy np.matmul, 128 x 128 (OpenBLAS, one thread)",
       "parity": "none",
-      "allocs": 13,
+      "allocs": 8,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2611,7 +2611,7 @@
       "func": "linalg.Det",
       "reference": "numpy np.linalg.det(a)",
       "parity": "none",
-      "allocs": 25,
+      "allocs": 24,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2629,7 +2629,7 @@
       "func": "linalg.Inv",
       "reference": "numpy np.linalg.inv(a)",
       "parity": "none",
-      "allocs": 30,
+      "allocs": 27,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2665,7 +2665,7 @@
       "func": "linalg.Pinv",
       "reference": "numpy np.linalg.pinv(a)",
       "parity": "none",
-      "allocs": 255,
+      "allocs": 230,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2719,7 +2719,7 @@
       "func": "linalg.Cholesky",
       "reference": "numpy np.linalg.cholesky(a)",
       "parity": "none",
-      "allocs": 20,
+      "allocs": 19,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2809,7 +2809,7 @@
       "func": "linalg.PinvRcond",
       "reference": "numpy np.linalg.pinv(a, rcond=1e-10)",
       "parity": "none",
-      "allocs": 255,
+      "allocs": 230,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2845,7 +2845,7 @@
       "func": "linalg.Eig",
       "reference": "numpy np.linalg.eig(a)",
       "parity": "none",
-      "allocs": 89,
+      "allocs": 80,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2863,7 +2863,7 @@
       "func": "linalg.Eigh",
       "reference": "numpy np.linalg.eigh(a)",
       "parity": "none",
-      "allocs": 166,
+      "allocs": 154,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2881,7 +2881,7 @@
       "func": "linalg.QR",
       "reference": "numpy np.linalg.qr(a)",
       "parity": "none",
-      "allocs": 93,
+      "allocs": 81,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2899,7 +2899,7 @@
       "func": "linalg.SVD",
       "reference": "numpy np.linalg.svd(a, full_matrices=False)",
       "parity": "none",
-      "allocs": 265,
+      "allocs": 241,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2917,7 +2917,7 @@
       "func": "linalg.Slogdet",
       "reference": "numpy np.linalg.slogdet(a)",
       "parity": "none",
-      "allocs": 29,
+      "allocs": 28,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2935,7 +2935,7 @@
       "func": "linalg.Solve",
       "reference": "numpy np.linalg.solve(a, b)",
       "parity": "none",
-      "allocs": 40,
+      "allocs": 37,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2953,7 +2953,7 @@
       "func": "linalg.Lstsq",
       "reference": "numpy np.linalg.lstsq(a, b, rcond=None)",
       "parity": "none",
-      "allocs": 276,
+      "allocs": 250,
       "go": {
         "median": 0,
         "q1": 0,
@@ -2971,7 +2971,7 @@
       "func": "linalg.LstsqRcond",
       "reference": "numpy np.linalg.lstsq(a, b, rcond=1e-10)",
       "parity": "none",
-      "allocs": 276,
+      "allocs": 250,
       "go": {
         "median": 0,
         "q1": 0,

internal/kernels/gemm.go

@@ -491,6 +491,15 @@ func (g *gemmKernel[T]) matMul(dst []T, a, b MatOf[T], m, k, n int) {
 	if w > m {
 		w = m
 	}
+	if w == 1 {
+		// One worker: the serial path, without the block list, the progress
+		// counters and the worker closure, which cost four allocations per call
+		// and linalg's blocked factorizations make one call per panel (#104).
+		buf := g.getPackBuf()
+		g.packGemmRows(0, m, dst, a, b, k, n, buf)
+		g.putPackBuf(buf)
+		return
+	}
 	// Dice the m rows into MR-aligned bands a few rows of micro-tiles tall, so
 	// there are several bands per worker. That oversubscription is what lets the
 	// dynamic counter balance the M-series big.LITTLE cores: fast P-cores grab

internal/kernels/gemm_oneworker_test.go

@@ -0,0 +1,63 @@
+package kernels
+
+import (
+	"math/rand/v2"
+	"runtime"
+	"testing"
+)
+
+// TestGemmOneWorkerSameBits pins the one-worker path of matMul (#104): at
+// GOMAXPROCS=1 a product above GemmThreshold takes packGemmRows instead of the
+// banded parallel path, and on non-integer data the two must agree bit for
+// bit (same tiles, same pc order), in float64 and float32, including m = 1,
+// which takes the same branch at any GOMAXPROCS.
+func TestGemmOneWorkerSameBits(t *testing.T) {
+	r := rand.New(rand.NewPCG(11, 12))
+	for _, s := range [][3]int{{97, 131, 83}, {1, 300, 7000}, {256, 64, 200}} {
+		m, k, n := s[0], s[1], s[2]
+		if m*n < GemmThreshold {
+			t.Fatalf("%v is below GemmThreshold %d: the case would not reach the branch", s, GemmThreshold)
+		}
+		a, b := make([]float64, m*k), make([]float64, k*n)
+		for i := range a {
+			a[i] = r.NormFloat64()
+		}
+		for i := range b {
+			b[i] = r.NormFloat64()
+		}
+		run := func(procs int) []float64 {
+			defer runtime.GOMAXPROCS(runtime.GOMAXPROCS(procs))
+			dst := make([]float64, m*n)
+			MatMulStridedP(dst, rowMajor(a, k), rowMajor(b, n), m, k, n)
+			return dst
+		}
+		one, many := run(1), run(4)
+		for i := range one {
+			if one[i] != many[i] {
+				t.Fatalf("float64 %v: element %d is %v on one worker, %v on four", s, i, one[i], many[i])
+			}
+		}
+		if !HasGemm32() {
+			continue
+		}
+		a32, b32 := make([]float32, m*k), make([]float32, k*n)
+		for i := range a32 {
+			a32[i] = float32(a[i])
+		}
+		for i := range b32 {
+			b32[i] = float32(b[i])
+		}
+		run32 := func(procs int) []float32 {
+			defer runtime.GOMAXPROCS(runtime.GOMAXPROCS(procs))
+			dst := make([]float32, m*n)
+			MatMul32StridedP(dst, MatOf[float32]{Data: a32, RS: k, CS: 1}, MatOf[float32]{Data: b32, RS: n, CS: 1}, m, k, n)
+			return dst
+		}
+		one32, many32 := run32(1), run32(4)
+		for i := range one32 {
+			if one32[i] != many32[i] {
+				t.Fatalf("float32 %v: element %d is %v on one worker, %v on four", s, i, one32[i], many32[i])
+			}
+		}
+	}
+}

linalg/blas.go

@@ -204,7 +204,9 @@ func gemm[T scalar](m, n, k int, alpha T, a, b mview[T], beta T, c []T, ldc int)
 // fastGemm adds alpha*A*B to C through the float64 packed GEMM: directly
 // for real operands, as four real products for complex ones.
 func fastGemm[T scalar](m, n, k int, alpha T, a, b mview[T], c []T, ldc int) {
-	tmp := make([]float64, m*n)
+	tp := getScratch(m * n)
+	defer putScratch(tp)
+	tmp := *tp
 	if cf, ok := any(c).([]float64); ok {
 		av, bv := any(a).(mview[float64]), any(b).(mview[float64])
 		kernels.MatMulStridedP(tmp, kmat(av), kmat(bv), m, k, n)

linalg/scratch.go

@@ -0,0 +1,30 @@
+package linalg
+
+import "sync"
+
+// scratchMax is the largest buffer (in float64s, 8 MiB) the scratch pool keeps:
+// a bigger one is allocated and dropped, so a single huge product does not pin
+// its temporary for the life of the process.
+const scratchMax = 1 << 20
+
+var scratchPool sync.Pool // of *[]float64
+
+// getScratch returns a zeroed float64 buffer of length n, reused across calls
+// when one is free: the blocked factorizations call fastGemm once per panel,
+// and its temporary was one allocation each time (#104).
+func getScratch(n int) *[]float64 {
+	if p, ok := scratchPool.Get().(*[]float64); ok && cap(*p) >= n {
+		*p = (*p)[:n]
+		clear(*p)
+		return p
+	}
+	s := make([]float64, n)
+	return &s
+}
+
+// putScratch returns p to the pool unless it is over scratchMax.
+func putScratch(p *[]float64) {
+	if cap(*p) <= scratchMax {
+		scratchPool.Put(p)
+	}
+}