diff --git a/.gitignore b/.gitignore index f4709d9..0ed9083 100644 --- a/.gitignore +++ b/.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 diff --git a/docs/perf-table.md b/docs/perf-table.md index 46f22d9..baab391 100644 --- a/docs/perf-table.md +++ b/docs/perf-table.md @@ -2,7 +2,7 @@ -- 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 | diff --git a/docs/perf.json b/docs/perf.json index e36f64f..d061ffa 100644 --- a/docs/perf.json +++ b/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, diff --git a/internal/kernels/gemm.go b/internal/kernels/gemm.go index 8ad3e0e..1358f0d 100644 --- a/internal/kernels/gemm.go +++ b/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 diff --git a/internal/kernels/gemm_oneworker_test.go b/internal/kernels/gemm_oneworker_test.go new file mode 100644 index 0000000..cb94786 --- /dev/null +++ b/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]) + } + } + } +} diff --git a/linalg/blas.go b/linalg/blas.go index 8b7688a..d35fcd6 100644 --- a/linalg/blas.go +++ b/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) diff --git a/linalg/linalg.test b/linalg/linalg.test deleted file mode 100755 index 4a5edad..0000000 Binary files a/linalg/linalg.test and /dev/null differ diff --git a/linalg/scratch.go b/linalg/scratch.go new file mode 100644 index 0000000..23aba2d --- /dev/null +++ b/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) + } +}