From 2a48ffcd4275f99ea7a6727c320ab20dbe877570 Mon Sep 17 00:00:00 2001 From: tannevaled Date: Sun, 4 Oct 2026 22:14:19 +0200 Subject: [PATCH 1/3] perf(gemm): persistent workers per call with a spinning barrier; double-buffered B Co-Authored-By: Claude Opus 5.5 --- internal/kernels/barrier_test.go | 42 ++++++++++++ internal/kernels/gemm.go | 113 +++++++++++++++++-------------- internal/kernels/gemm_test.go | 13 ++-- internal/kernels/parallel.go | 33 +++++++++ 4 files changed, 142 insertions(+), 59 deletions(-) create mode 100644 internal/kernels/barrier_test.go diff --git a/internal/kernels/barrier_test.go b/internal/kernels/barrier_test.go new file mode 100644 index 0000000..3a0902d --- /dev/null +++ b/internal/kernels/barrier_test.go @@ -0,0 +1,42 @@ +package kernels + +import ( + "sync" + "sync/atomic" + "testing" +) + +// TestSpinBarrier: no goroutine leaves a generation before all have arrived, +// over many back-to-back generations, both spinning and yielding. +func TestSpinBarrier(t *testing.T) { + for _, iters := range []int{spinIters, 0} { + saved := spinIters + spinIters = iters + const n, rounds = 4, 2000 + bar := newSpinBarrier(n) + var arrived atomic.Int64 + var bad atomic.Int64 + var wg sync.WaitGroup + for g := 0; g < n; g++ { + wg.Add(1) + go func() { + defer wg.Done() + for r := int64(1); r <= rounds; r++ { + arrived.Add(1) + bar.wait() + // Everyone has arrived for round r, and nobody can have + // arrived for round r+1 before this goroutine has. + if got := arrived.Load(); got < r*n || got > r*n+n-1 { + bad.Add(1) + } + bar.wait() + } + }() + } + wg.Wait() + spinIters = saved + if bad.Load() != 0 { + t.Fatalf("spinIters=%d: %d early departures", iters, bad.Load()) + } + } +} diff --git a/internal/kernels/gemm.go b/internal/kernels/gemm.go index 7caad0e..293c06e 100644 --- a/internal/kernels/gemm.go +++ b/internal/kernels/gemm.go @@ -190,30 +190,6 @@ func packB(b Mat, pc, kc, jc, nc int, pb []float64) { } } -// packBParMin is the B-block size, in elements, from which packBP splits the -// packing across workers (64 Ki doubles = 512 KiB). A var so tests can pin it. -var packBParMin = 1 << 16 - -// packBP is packB with the NR-wide panels split across w goroutines. MatMulP -// packs the shared B block once per (jc, pc) round, and done by one goroutine -// it is a serial step while every other worker waits at the round's start: -// a 1 MiB strided copy, eight times in a 1024^2 product. -// -// Below packBParMin elements the goroutine launch costs more than the copy -// (measured on M4: a 128x128 block packed in parallel made 96^3..128^2 -// products 8-15% slower), so a small block is packed by the caller alone. -func packBP(b Mat, pc, kc, jc, nc int, pb []float64, w int) { - if kc*nc < packBParMin { - packB(b, pc, kc, jc, nc, pb) - return - } - panels := (nc + NR - 1) / NR - parallelFor(panels, min(w, panels), func(lo, hi int) { - j0 := lo * NR - packB(b, pc, kc, jc+j0, min(hi*NR, nc)-j0, pb[j0*kc:]) - }) -} - // packA copies A[ic:ic+mc][pc:pc+kc] into MR-tall row panels: for each panel ir // the layout is pa[ir*kc + p*MR + r] = A(ic+ir+r, pc+p), with the r >= mr rows // of an edge panel zero-filled so the micro-kernel reads a full MR. @@ -339,41 +315,76 @@ func MatMulStridedP(dst []float64, a, b Mat, m, k, n int) { w = nBands // never spawn more goroutines than there is work for } - // One shared B-pack buffer (packed once per jc/pc block) plus one private - // A-pack buffer per worker. The shared pb is sized like a normal pooled B - // buffer; getPackBuf grows it for any test-raised block sizes. - shared := getPackBuf() - pb := shared.pb + // The workers are started once per call and stay for every (jc, pc) block, + // meeting at a spinning barrier between blocks instead of being forked and + // joined per block. A fork/join wakes the parked worker threads one after + // another; measured on a 16-vCPU Zen 3 guest it costs 26 µs at 2 workers, + // 165 µs at 8 and 330 µs at 16, against ~0.8 ms of work per block of a + // 1024² product. With two fork/joins per block (pack B, then the bands) + // that was ~40% of the wall time, and why the CPUs were busy only 7 of 16. + // + // pb is double-buffered so one barrier per block suffices: block i packs + // into pbs[i%2] while nobody can still be reading it, because reading it + // for block i-2 ended before everyone passed block i-1's barrier. + shared, shared2 := getPackBuf(), getPackBuf() defer packPool.Put(shared) + defer packPool.Put(shared2) + pbs := [2][]float64{shared.pb, shared2.pb} + type block struct{ jc, nc, pc, kc int } + var blocks []block for jc := 0; jc < n; jc += blockNC { // L3 column block of B - nc := min(blockNC, n-jc) for pc := 0; pc < k; pc += blockKC { // L2 contraction block - kc := min(blockKC, k-pc) - packBP(b, pc, kc, jc, nc, pb, w) // pack B once, by all workers + blocks = append(blocks, block{jc, min(blockNC, n-jc), pc, min(blockKC, k-pc)}) + } + } + // Each block has its own band cursor (dynamic banding within a block) and + // its own B-panel cursor (packing shared by whoever is free). + bandNext := make([]atomic.Int64, len(blocks)) + packNext := make([]atomic.Int64, len(blocks)) + bar := newSpinBarrier(w) - var next atomic.Int64 // shared band cursor - var wg sync.WaitGroup - for g := 0; g < w; g++ { - wg.Add(1) - go func() { - defer wg.Done() - ab := getPackBuf() - for { - bi := int(next.Add(1)) - 1 - if bi >= nBands { - break - } - r0 := bi * bandRows - r1 := min(r0+bandRows, m) - gemmBand(r0, r1, jc, nc, pc, kc, dst, a, n, ab.pa, pb) - } - packPool.Put(ab) - }() + work := func() { + ab := getPackBuf() + defer packPool.Put(ab) + for i, bl := range blocks { + pb := pbs[i%2] + panels := (bl.nc + NR - 1) / NR + for { + // B panels in chunks, so packing is shared without one atomic + // per 8-column panel. + const chunk = 4 + p0 := int(packNext[i].Add(chunk)) - chunk + if p0 >= panels { + break + } + j0, j1 := p0*NR, min((p0+chunk)*NR, bl.nc) + packB(b, bl.pc, bl.kc, bl.jc+j0, j1-j0, pb[j0*bl.kc:]) } - wg.Wait() + bar.wait() // pb for block i is complete + for { + bi := int(bandNext[i].Add(1)) - 1 + if bi >= nBands { + break + } + r0 := bi * bandRows + gemmBand(r0, min(r0+bandRows, m), bl.jc, bl.nc, bl.pc, bl.kc, dst, a, n, ab.pa, pb) + } + // No barrier here: the next block packs into the other buffer, and + // the barrier after that packing orders this block's reads of pb + // before the packing of block i+2 into it. Tiles of C are still + // accumulated in pc order: a band of block i+1 starts only after + // every worker has passed block i+1's barrier, i.e. finished its + // bands of block i. } } + var wg sync.WaitGroup + for g := 1; g < w; g++ { + wg.Add(1) + go func() { defer wg.Done(); work() }() + } + work() // the caller is worker 0 + wg.Wait() } // MatMul computes dst = a(m x k) * b(k x n) serially with the packed GEMM (the diff --git a/internal/kernels/gemm_test.go b/internal/kernels/gemm_test.go index 014c9ea..60774d9 100644 --- a/internal/kernels/gemm_test.go +++ b/internal/kernels/gemm_test.go @@ -140,14 +140,11 @@ func TestGemmBandSizing(t *testing.T) { }) } -// TestGemmParallelPackB forces packBP's split path (normally only for blocks of -// 64 Ki elements or more) on small products, so the NR-panel ranges handed to -// each goroutine — including a ragged last panel and more workers than panels -// — must reassemble the exact packed block the serial packB writes. -func TestGemmParallelPackB(t *testing.T) { - o := packBParMin - packBParMin = 1 - defer func() { packBParMin = o }() +// TestGemmSharedPackB runs MatMulP's parallel path on small products, where +// the workers share the packing of B in chunks of NR-wide panels: the chunks — +// including a ragged last panel and more workers than panels — must reassemble +// the exact packed block the serial packB writes, in both pb buffers. +func TestGemmSharedPackB(t *testing.T) { withMaxProcs(4, func() { withThresholds(1<<14, 1, func() { // force MatMulP's parallel path for _, s := range []struct{ m, k, n int }{ diff --git a/internal/kernels/parallel.go b/internal/kernels/parallel.go index 835406b..11dcef5 100644 --- a/internal/kernels/parallel.go +++ b/internal/kernels/parallel.go @@ -3,6 +3,7 @@ package kernels import ( "runtime" "sync" + "sync/atomic" ) // Multicore kernels. @@ -265,3 +266,35 @@ func MaxP(a []float64) float64 { return reduceP(a, maxSIMD) } // ParThreshold, via the SIMD min kernel; result identical to the serial Min, // NaN-propagating like numpy.min. func MinP(a []float64) float64 { return reduceP(a, minSIMD) } + +// spinBarrier is a reusable barrier for n goroutines that are all running, as +// the GEMM workers are between their blocks. Waiting is a spin, then a spin +// that yields (runtime.Gosched) so an oversubscribed machine still makes +// progress; it never parks, because parking is exactly the cost it replaces: +// waking n parked threads one after another. +// +// The arrival count only ever grows: the c-th arrival (from 1) belongs to +// generation (c-1)/n and leaves once the count reaches the end of it. With no +// reset there is no window in which a released goroutine could arrive before +// the count is cleared. +type spinBarrier struct { + n int64 + count atomic.Int64 +} + +func newSpinBarrier(n int) *spinBarrier { return &spinBarrier{n: int64(n)} } + +// spinIters is how many times a waiter polls before it starts yielding. A var +// so tests can make every wait yield. +var spinIters = 1 << 10 + +// wait blocks until n goroutines have called it for the current generation. +func (b *spinBarrier) wait() { + c := b.count.Add(1) + end := ((c-1)/b.n + 1) * b.n + for i := 0; b.count.Load() < end; i++ { + if i >= spinIters { + runtime.Gosched() + } + } +} From 0ef1ad3a083dd79b3dd4f67bcfec8cd8d1ad1bca Mon Sep 17 00:00:00 2001 From: tannevaled Date: Sun, 4 Oct 2026 22:34:25 +0200 Subject: [PATCH 2/3] perf(gemm): synchronise on work done, not workers arrived; split pa/pb pools Co-Authored-By: Claude Opus 5.5 --- internal/kernels/barrier_test.go | 57 +++++++------ internal/kernels/gemm.go | 134 ++++++++++++++++--------------- internal/kernels/parallel.go | 37 +++------ 3 files changed, 109 insertions(+), 119 deletions(-) diff --git a/internal/kernels/barrier_test.go b/internal/kernels/barrier_test.go index 3a0902d..6378cbd 100644 --- a/internal/kernels/barrier_test.go +++ b/internal/kernels/barrier_test.go @@ -1,42 +1,41 @@ package kernels import ( - "sync" "sync/atomic" "testing" ) -// TestSpinBarrier: no goroutine leaves a generation before all have arrived, -// over many back-to-back generations, both spinning and yielding. -func TestSpinBarrier(t *testing.T) { +// TestSpinUntil waits on a condition another goroutine makes true, both +// spinning and yielding from the first poll. +func TestSpinUntil(t *testing.T) { for _, iters := range []int{spinIters, 0} { saved := spinIters spinIters = iters - const n, rounds = 4, 2000 - bar := newSpinBarrier(n) - var arrived atomic.Int64 - var bad atomic.Int64 - var wg sync.WaitGroup - for g := 0; g < n; g++ { - wg.Add(1) - go func() { - defer wg.Done() - for r := int64(1); r <= rounds; r++ { - arrived.Add(1) - bar.wait() - // Everyone has arrived for round r, and nobody can have - // arrived for round r+1 before this goroutine has. - if got := arrived.Load(); got < r*n || got > r*n+n-1 { - bad.Add(1) - } - bar.wait() - } - }() - } - wg.Wait() + var flag atomic.Bool + go func() { + for i := 0; i < 1000; i++ { + _ = i * i + } + flag.Store(true) + }() + spinUntil(flag.Load) spinIters = saved - if bad.Load() != 0 { - t.Fatalf("spinIters=%d: %d early departures", iters, bad.Load()) - } } } + +// TestGemmParallelYielding runs the parallel GEMM with every wait yielding at +// once, as on an oversubscribed machine, and with more workers than CPUs. +func TestGemmParallelYielding(t *testing.T) { + saved := spinIters + spinIters = 0 + defer func() { spinIters = saved }() + withMaxProcs(8, func() { + withThresholds(1<<14, 1, func() { + for _, s := range []struct{ m, k, n int }{{50, 70, 90}, {13, 300, 7}, {97, 31, 600}} { + a, b := intMat(s.m, s.k, 3), intMat(s.k, s.n, 4) + assertGemm(t, s.m, s.k, s.n, a, b) + withBlocks(MR, 8, 2*NR, func() { assertGemm(t, s.m, s.k, s.n, a, b) }) + } + }) + }) +} diff --git a/internal/kernels/gemm.go b/internal/kernels/gemm.go index 293c06e..8eadf4b 100644 --- a/internal/kernels/gemm.go +++ b/internal/kernels/gemm.go @@ -62,9 +62,11 @@ type Mat struct { // rowMajor views a contiguous row-major matrix with ld elements per row. func rowMajor(data []float64, ld int) Mat { return Mat{Data: data, RS: ld, CS: 1} } -// packBuf is one worker's reusable pair of pack buffers: paBuf for the MC*KC A -// panel, pbBuf for the KC*NC B panel. Pooling them keeps the per-call GEMM -// allocation-free above the goroutine launch, which matters at small n. +// packBuf is a pair of pack buffers: pa for the MC*KC A panel, pb for the KC*NC +// B panel. The two come from separate pools, because the parallel GEMM needs +// one pa per worker but only two pb in all: drawing a pair per worker held +// 16 MiB of unused B buffers at 16 workers, which the garbage collector empties +// out of the pool and the next call zeroes again. type packBuf struct { pa []float64 pb []float64 @@ -83,27 +85,29 @@ func roundUp(x, m int) int { return (x + m - 1) / m * m } func paCap() int { return roundUp(blockMC, MR) * blockKC } func pbCap() int { return blockKC * roundUp(blockNC, NR) } -var packPool = sync.Pool{New: func() any { - return &packBuf{ - pa: make([]float64, paCap()), - pb: make([]float64, pbCap()), - } -}} +var paPool, pbPool sync.Pool -// getPackBuf returns a pooled buffer pair sized for the current block params. If -// the pooled buffers are too small for (possibly test-tuned) larger blocks it -// reallocates them, so a test that raises blockMC/KC/NC still gets valid scratch. -func getPackBuf() *packBuf { - b := packPool.Get().(*packBuf) - if cap(b.pa) < paCap() { - b.pa = make([]float64, paCap()) - } - if cap(b.pb) < pbCap() { - b.pb = make([]float64, pbCap()) +// getPA and getPB return a pooled pack buffer sized for the current block +// parameters, allocating one when the pool is empty or its buffer is too small +// (a test may raise blockMC/KC/NC). putPA and putPB return it. +func getPA() []float64 { return getScratch(&paPool, paCap()) } +func getPB() []float64 { return getScratch(&pbPool, pbCap()) } +func putPA(b []float64) { paPool.Put(&b) } +func putPB(b []float64) { pbPool.Put(&b) } + +func getScratch(p *sync.Pool, n int) []float64 { + if b, ok := p.Get().(*[]float64); ok && cap(*b) >= n { + return (*b)[:n] } - b.pa = b.pa[:paCap()] - b.pb = b.pb[:pbCap()] - return b + return make([]float64, n) +} + +// getPackBuf returns a pooled pair, for the single-worker paths. +func getPackBuf() *packBuf { return &packBuf{pa: getPA(), pb: getPB()} } + +func putPackBuf(b *packBuf) { + putPA(b.pa) + putPB(b.pb) } // packGemmRows computes the row band [r0,r1) of dst = a(m x k) * b(k x n) with @@ -288,7 +292,7 @@ func MatMulStridedP(dst []float64, a, b Mat, m, k, n int) { if m*n < GemmThreshold { buf := getPackBuf() packGemmRows(0, m, dst, a, b, k, n, buf) - packPool.Put(buf) + putPackBuf(buf) return } w := runtime.GOMAXPROCS(0) @@ -316,66 +320,70 @@ func MatMulStridedP(dst []float64, a, b Mat, m, k, n int) { } // The workers are started once per call and stay for every (jc, pc) block, - // meeting at a spinning barrier between blocks instead of being forked and - // joined per block. A fork/join wakes the parked worker threads one after - // another; measured on a 16-vCPU Zen 3 guest it costs 26 µs at 2 workers, - // 165 µs at 8 and 330 µs at 16, against ~0.8 ms of work per block of a - // 1024² product. With two fork/joins per block (pack B, then the bands) - // that was ~40% of the wall time, and why the CPUs were busy only 7 of 16. + // instead of being forked and joined per block. A fork/join wakes the + // parked worker threads one after another; measured on a 16-vCPU Zen 3 + // guest it costs 26 µs at 2 workers, 165 µs at 8 and 330 µs at 16, against + // ~0.8 ms of work per block of a 1024² product, and there were two per + // block (pack B, then the bands): ~40% of the wall time. + // + // The workers synchronise on work done, never on workers arrived: a block's + // bands start once all of its B panels are packed and all bands of the + // previous block are done. A worker still waking up holds nobody back, so + // the wake-up chain stays off the critical path; a barrier that waited for + // every worker put it right back on it for small products (96³ on 8 POWER9 + // cores: 0.57x). Claimed work never waits, so the scheme cannot deadlock. // - // pb is double-buffered so one barrier per block suffices: block i packs - // into pbs[i%2] while nobody can still be reading it, because reading it - // for block i-2 ended before everyone passed block i-1's barrier. - shared, shared2 := getPackBuf(), getPackBuf() - defer packPool.Put(shared) - defer packPool.Put(shared2) - pbs := [2][]float64{shared.pb, shared2.pb} + // pb is double-buffered, so block i can be packed while block i-1 is being + // multiplied. Tiles of C are still accumulated in pc order, since block i's + // bands wait for block i-1's. + pbs := [2][]float64{getPB(), getPB()} + defer putPB(pbs[0]) + defer putPB(pbs[1]) - type block struct{ jc, nc, pc, kc int } + type block struct{ jc, nc, pc, kc, chunks int } + const chunk = 4 // B panels packed per claim var blocks []block for jc := 0; jc < n; jc += blockNC { // L3 column block of B + nc := min(blockNC, n-jc) + panels := (nc + NR - 1) / NR for pc := 0; pc < k; pc += blockKC { // L2 contraction block - blocks = append(blocks, block{jc, min(blockNC, n-jc), pc, min(blockKC, k-pc)}) + blocks = append(blocks, block{jc, nc, pc, min(blockKC, k-pc), (panels + chunk - 1) / chunk}) } } - // Each block has its own band cursor (dynamic banding within a block) and - // its own B-panel cursor (packing shared by whoever is free). - bandNext := make([]atomic.Int64, len(blocks)) - packNext := make([]atomic.Int64, len(blocks)) - bar := newSpinBarrier(w) + type progress struct { + packNext, packDone, bandNext, bandDone atomic.Int64 + _ [32]byte // keep blocks' counters off one cache line + } + prog := make([]progress, len(blocks)) + bandsDone := func(i int) bool { return i < 0 || prog[i].bandDone.Load() == int64(nBands) } work := func() { - ab := getPackBuf() - defer packPool.Put(ab) + pa := getPA() + defer putPA(pa) for i, bl := range blocks { - pb := pbs[i%2] - panels := (bl.nc + NR - 1) / NR + // pbs[i%2] is free: its last reader was block i-2, and to get here + // this worker passed block i-1's wait, which requires every band of + // block i-2 done. + pr, pb := &prog[i], pbs[i%2] for { - // B panels in chunks, so packing is shared without one atomic - // per 8-column panel. - const chunk = 4 - p0 := int(packNext[i].Add(chunk)) - chunk - if p0 >= panels { + c := int(pr.packNext.Add(1)) - 1 + if c >= bl.chunks { break } - j0, j1 := p0*NR, min((p0+chunk)*NR, bl.nc) + j0, j1 := c*chunk*NR, min((c+1)*chunk*NR, bl.nc) packB(b, bl.pc, bl.kc, bl.jc+j0, j1-j0, pb[j0*bl.kc:]) + pr.packDone.Add(1) } - bar.wait() // pb for block i is complete + spinUntil(func() bool { return pr.packDone.Load() == int64(bl.chunks) && bandsDone(i-1) }) for { - bi := int(bandNext[i].Add(1)) - 1 + bi := int(pr.bandNext.Add(1)) - 1 if bi >= nBands { break } r0 := bi * bandRows - gemmBand(r0, min(r0+bandRows, m), bl.jc, bl.nc, bl.pc, bl.kc, dst, a, n, ab.pa, pb) + gemmBand(r0, min(r0+bandRows, m), bl.jc, bl.nc, bl.pc, bl.kc, dst, a, n, pa, pb) + pr.bandDone.Add(1) } - // No barrier here: the next block packs into the other buffer, and - // the barrier after that packing orders this block's reads of pb - // before the packing of block i+2 into it. Tiles of C are still - // accumulated in pc order: a band of block i+1 starts only after - // every worker has passed block i+1's barrier, i.e. finished its - // bands of block i. } } var wg sync.WaitGroup @@ -394,7 +402,7 @@ func MatMulStridedP(dst []float64, a, b Mat, m, k, n int) { func MatMul(dst, a, b []float64, m, k, n int) { buf := getPackBuf() packGemmRows(0, m, dst, rowMajor(a, k), rowMajor(b, n), k, n, buf) - packPool.Put(buf) + putPackBuf(buf) } // dotRange returns sum(a[i]*b[i]) over equal-length slices, using four diff --git a/internal/kernels/parallel.go b/internal/kernels/parallel.go index 11dcef5..f13ae51 100644 --- a/internal/kernels/parallel.go +++ b/internal/kernels/parallel.go @@ -3,7 +3,6 @@ package kernels import ( "runtime" "sync" - "sync/atomic" ) // Multicore kernels. @@ -267,34 +266,18 @@ func MaxP(a []float64) float64 { return reduceP(a, maxSIMD) } // NaN-propagating like numpy.min. func MinP(a []float64) float64 { return reduceP(a, minSIMD) } -// spinBarrier is a reusable barrier for n goroutines that are all running, as -// the GEMM workers are between their blocks. Waiting is a spin, then a spin -// that yields (runtime.Gosched) so an oversubscribed machine still makes -// progress; it never parks, because parking is exactly the cost it replaces: -// waking n parked threads one after another. -// -// The arrival count only ever grows: the c-th arrival (from 1) belongs to -// generation (c-1)/n and leaves once the count reaches the end of it. With no -// reset there is no window in which a released goroutine could arrive before -// the count is cleared. -type spinBarrier struct { - n int64 - count atomic.Int64 -} - -func newSpinBarrier(n int) *spinBarrier { return &spinBarrier{n: int64(n)} } - -// spinIters is how many times a waiter polls before it starts yielding. A var -// so tests can make every wait yield. -var spinIters = 1 << 10 - -// wait blocks until n goroutines have called it for the current generation. -func (b *spinBarrier) wait() { - c := b.count.Add(1) - end := ((c-1)/b.n + 1) * b.n - for i := 0; b.count.Load() < end; i++ { +// spinUntil waits for done to become true, which other running goroutines +// will make so: a spin, then a spin that yields (runtime.Gosched) so an +// oversubscribed machine still makes progress. It never parks, because parking +// is exactly the cost the GEMM workers avoid: waking threads one after another. +func spinUntil(done func() bool) { + for i := 0; !done(); i++ { if i >= spinIters { runtime.Gosched() } } } + +// spinIters is how many times spinUntil polls before it starts yielding. A var +// so tests can make every wait yield. +var spinIters = 1 << 10 From a8d77e8685efaf03bd22d92855f9b310fada0aed Mon Sep 17 00:00:00 2001 From: tannevaled Date: Sun, 4 Oct 2026 22:49:44 +0200 Subject: [PATCH 3/3] docs(perf): parallel GEMM diagnosis and the work-based driver; correct the profile-based claim Co-Authored-By: Claude Opus 5.5 --- README.md | 5 +-- docs/perf.md | 95 ++++++++++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 92 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index d8467b0..1d6b723 100644 --- a/README.md +++ b/README.md @@ -86,8 +86,9 @@ is wrong on subnormals beat NumPy at every size measured. - **At parity:** `SumAxis(0)`, `Inner` against 16-thread OpenBLAS, `Log` on one core. -- **Slower:** `MatMul` against 16-thread OpenBLAS (0.55× at 1024², 0.3× at - 256² and below); `Dot`/mat·vec against NumPy's threaded BLAS (0.16×); `Exp` +- **Slower:** `MatMul` against 16-thread OpenBLAS at 512² and up (0.42× at + 512², 0.51× at 1024², measured in one run with OpenBLAS; at 256² ours is + 1.44× faster); `Dot`/mat·vec against NumPy's threaded BLAS (0.16×); `Exp` below 256 Ki elements (0.82×); the *allocating* forms of elementwise ops up to 256 Ki elements, `Concatenate`/`Stack` and slice copies outside a `Workspace` (0.2–0.5×: Go's allocator and garbage collector, which diff --git a/docs/perf.md b/docs/perf.md index c5299b8..10cfb6a 100644 --- a/docs/perf.md +++ b/docs/perf.md @@ -411,9 +411,13 @@ below 256²). Where the parallel loss is, measured on the Zen 3: each ran at 41 GFLOP/s (656 in total), while the parallel 1024² GEMM ran its kernel at 15–17 GFLOP/s per core. Per-core efficiency already drops from 40 to ~30 GFLOP/s at 2 cores. -- It is not scheduling overhead: 83% of the parallel profile is the - micro-kernel itself, so the kernel runs slower when cores share work. Why is - **not established**. The shared packed B panel is the suspect (on this part +- ~~It is not scheduling overhead: 83% of the parallel profile is the + micro-kernel itself, so the kernel runs slower when cores share work.~~ + **Wrong, corrected the same day:** a profile of cycles samples only the cores + that run, never the ones that wait. `perf stat` shows 7.1 of 16 CPUs busy + during a 1024² product: the loss was mostly waiting, and the fix is in + [Parallel GEMM: start once, wait on work](#parallel-gemm-start-once-wait-on-work-2026-10-04). + Why the kernel itself also runs slower in parallel is still **not established**. The shared packed B panel is the suspect (on this part the last-level cache is split between core complexes, so the panel is partly remote), but the drop already shows at 2 cores of one complex, which that explanation does not cover. @@ -426,9 +430,8 @@ below 256²). Where the parallel loss is, measured on the Zen 3: rows·cols·k multiply-adds, and a strided pack costs far more per element than an FMA, so tiles under ~256×256 lose. -The next step is to find where the kernel stalls at 2 cores (hardware -counters: cache misses per FMA, serial against parallel) before choosing a -lever; then an AVX-512 kernel for the hosts that have it. +The hardware counters were read the same day; see the next-but-one +section. An AVX-512 kernel for the hosts that have it remains open. Correctness was checked on real hardware for five of the six 64-bit targets: amd64 (Zen 3, FMA path), arm64 (cfarm185), ppc64le (cfarm120), riscv64 (cfarm94) @@ -637,6 +640,86 @@ 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. +## Parallel GEMM: start once, wait on work (2026-10-04) + +**Diagnosis, from hardware counters** (`perf stat`, Zen 3 guest cfarm421, +steady state, 1024², 3 s per configuration): + +| workers | CPUs busy | IPC | instructions | +|---|---|---|---| +| 1 | 1.00 | 3.35 | 37.5 G | +| 2 | 1.91 | 2.00 | 37.6 G | +| 8 | 5.58 | 2.16 | 172 G (192 products) | +| 16 | 7.08 | 2.28 | 193 G (226 products) | + +The instruction count does not grow, so nobody busy-waits. L1 misses do not +grow either. What changes is that at 16 workers the CPUs are busy only 44% of +the time. The driver forked and joined goroutines twice per (jc, pc) block of +the product, once to pack B and once for the row bands, 16 times in a 1024² +product. A fork/join wakes parked threads one after another. A micro-benchmark +on the same guest (16 goroutines, 700 µs of work each) measures its cost at +**26 µs with 2 workers, 165 µs with 8, 330 µs with 16**. 16 × 330 µs is ~40% +of a 12.5 ms product. + +Two hypotheses were tested and **refuted** first: + +- SMT siblings. The guest has 128 vCPUs, as many as the host has threads. + Two serial GEMMs pinned on vCPU 16 and any other vCPU slow each other by + 0–10%, never by half. +- A block of A larger than L2. With MC = 256, the A block is 512 KiB, all of + Zen 3's L2, against BLIS's MC = 72. Sweeping MC over {72, 120, 144, 192, 256} + changed nothing beyond the noise at 2 workers, and gave +17% at 16 workers + for 1024² but −8% for 512². + +**The change.** The workers start once per product. They synchronise on +**work done**, never on workers arrived: a block's bands start once all of its +B panels are packed and all bands of the previous block are done, and the +packing is shared in claims of four NR-wide panels. A first version used a +barrier that waited for every worker, and it put the wake-up chain straight +back on the critical path for small products: 96³ on 8 POWER9 cores ran at +0.57×. B is double-buffered, so block i+1 is packed while block i is still +being multiplied. Each C tile is still accumulated in pc order, so the result +is bit-identical to the serial product. The pack buffers moved to two pools, +one for A and one for B: a pair per worker had held 16 MiB of unused B buffers +at 16 workers, which the collector empties and the next call zeroes again +(`memclr` page faults, 20% of the cycles in a KVM spinlock in one profile). + +**Measured.** Every figure is the median of interleaved runs against the +previous driver. + +| | 2 workers | 8 workers | 16 workers | +|---|---|---|---| +| Zen 3 guest, 1024² | 1.06× | 1.33× | 1.39× | +| Zen 3 guest, 512² | 1.20× | 1.20× | 1.33× | +| Zen 3 guest, 256² | 0.92× | 1.19× | 1.24× | +| Zen 3 guest, 128² | **0.89×** | 1.50× | 1.16× | +| Zen 3 guest, (1024×256)·(256×1024) | **0.82×** | 1.07× | 0.89× (p25 1.00×) | +| POWER9 bare metal, 512² | 1.34× | 1.37× | | +| POWER9 bare metal, 256² | 1.61× | 1.78× | | +| POWER9 bare metal, (1024×256)·(256×1024) | 1.13× | 1.18× | | + +The Zen 3 runs were 12 rounds and the POWER9 runs 6, except 24 for 128² and +256² at 8 workers. Those two are bimodal on that machine, in both drivers, and +6 rounds had shown 256² at 0.72×. POWER9 used one thread per core (CPUs 0, 4, +8, …). + +**Still behind.** On the Zen 3 guest, two workers lose 8–18% on mid-size +products, where a fork/join costs only 26 µs. Bare-metal POWER9 gains at two +workers, so this is not a property of the scheme, and no second code path was +added for it. Against OpenBLAS with 16 threads, measured three ways in one run: + +| n | OpenBLAS-16 | previous driver | new driver | +|---|---|---|---| +| 256 | 708 µs | 1.06× | **1.44×** | +| 512 | 789 µs | 0.29× | 0.42× | +| 1024 | 4.60 ms | 0.43× | 0.51× | + +(× = OpenBLAS time / ours, so higher is faster). On this guest, absolute times +drift through the day: the previous driver measured 0.54–0.66× of OpenBLAS in +the morning. Only ratios taken in the same run are comparable. Not measured on +Apple M-series, whose machine was shared with other jobs at the time. The bands +stay dynamic, so its slower efficiency cores are handled as before. + ## SIMD coverage - **amd64 (SSE2)** ships hand-vectorized `sum` (4-accumulator `ADDPD`), `sqrt`