From f650cb8bb7198a80ef6df0da84ce740188707cb4 Mon Sep 17 00:00:00 2001 From: andig Date: Fri, 3 Jul 2026 13:01:03 +0200 Subject: [PATCH 1/4] simplex: extract axpy, add experimental Go 1.27 SIMD backend MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit pivot's Binv row-update and duals' accumulate loop are ~76% of solver CPU time (profiled via pprof on a 300x300 random LP). Both are the same dense AXPY pattern (dst[k] += factor*src[k]), now factored into a shared axpy() with two implementations selected by build tag: - axpy_generic.go (default): the original scalar loop, unchanged behavior. - axpy_simd.go (//go:build goexperiment.simd): vectorized via the new portable `simd` package, gated behind GOEXPERIMENT=simd (Go 1.27rc1+). Benchmarked scalar vs SIMD (benchstat, n=10, Apple M4 Pro/arm64 NEON, 128-bit = 2 float64 lanes): ColdSolve_50x50 +5.7% slower (call/bounds-check overhead dominates) ColdSolve_150x150 -10.6% faster ColdSolve_300x300 -2.8% faster geomean -2.8% Mixed, modest result — not a clean win. Not recommending this for the default build: it regresses at small problem sizes (likely the common case), and the whole simd package is still experimental on an unreleased toolchain (1.27rc1). Landing it behind the build tag costs nothing (zero effect on default `go build`/`go test`) and gives a starting point to revisit once Go 1.27 ships and/or wider (256/512-bit) SIMD is targeted. Co-Authored-By: Claude Sonnet 5 --- simplex/axpy_generic.go | 11 +++++++++++ simplex/axpy_simd.go | 22 ++++++++++++++++++++++ simplex/simplex.go | 10 ++-------- 3 files changed, 35 insertions(+), 8 deletions(-) create mode 100644 simplex/axpy_generic.go create mode 100644 simplex/axpy_simd.go diff --git a/simplex/axpy_generic.go b/simplex/axpy_generic.go new file mode 100644 index 0000000..73cefef --- /dev/null +++ b/simplex/axpy_generic.go @@ -0,0 +1,11 @@ +//go:build !goexperiment.simd + +package simplex + +// axpy computes dst[k] += factor*src[k] for all k. Scalar fallback used by +// default builds (Go 1.27+ with GOEXPERIMENT=simd gets a vectorized version). +func axpy(dst, src []float64, factor float64) { + for k := range dst { + dst[k] += factor * src[k] + } +} diff --git a/simplex/axpy_simd.go b/simplex/axpy_simd.go new file mode 100644 index 0000000..1544ac7 --- /dev/null +++ b/simplex/axpy_simd.go @@ -0,0 +1,22 @@ +//go:build goexperiment.simd + +package simplex + +import "simd" + +// axpy computes dst[k] += factor*src[k] for all k, vectorized via the +// experimental simd package (Go 1.27+, GOEXPERIMENT=simd). +func axpy(dst, src []float64, factor float64) { + fv := simd.BroadcastFloat64s(factor) + width := fv.Len() + n := len(dst) + k := 0 + for ; k+width <= n; k += width { + d := simd.LoadFloat64s(dst[k : k+width]) + s := simd.LoadFloat64s(src[k : k+width]) + s.MulAdd(fv, d).Store(dst[k : k+width]) + } + for ; k < n; k++ { + dst[k] += factor * src[k] + } +} diff --git a/simplex/simplex.go b/simplex/simplex.go index a25eae2..bb18b00 100644 --- a/simplex/simplex.go +++ b/simplex/simplex.go @@ -221,9 +221,7 @@ func duals(st *State, cost []float64, m int) []float64 { if cb == 0 { continue } - for k := 0; k < m; k++ { - y[k] += cb * st.binv[i][k] - } + axpy(y, st.binv[i], cb) } return y } @@ -472,11 +470,7 @@ func (lp *LP) pivot(st *State, q int, dir float64, a []float64, t float64, leave if i == leaveRow || a[i] == 0 { continue } - factor := a[i] - row := st.binv[i] - for k := 0; k < lp.m; k++ { - row[k] -= factor * rowR[k] - } + axpy(st.binv[i], rowR, -a[i]) } st.basicOf[leaveRow] = q st.status[q] = basic From 0ce9f7eaaff230471b91d285a79d9f21b485a76b Mon Sep 17 00:00:00 2001 From: andig Date: Fri, 3 Jul 2026 13:32:43 +0200 Subject: [PATCH 2/4] simplex: use simd/archsimd directly on arm64, drop portable dispatch cost The portable simd package's per-call CPU-feature dispatch stub was eating most of the 2-lane NEON win at small m (previous PR update: +5.7% regression at 50x50). Switching arm64 to simd/archsimd's Float64x2 directly (same 128-bit/2-lane hardware, no dispatch stub) turns this into a clean win across every size: ColdSolve_50x50 -11.6% ColdSolve_150x150 -22.8% ColdSolve_300x300 -19.1% geomean -18.0% vs the prior portable-simd geomean of -2.8%. simd/archsimd is AMD64/ arm64 only for now (per its doc), so axpy_simd.go (portable) stays as the fallback for other GOEXPERIMENT=simd architectures. Co-Authored-By: Claude Sonnet 5 --- simplex/axpy_archsimd_arm64.go | 21 +++++++++++++++++++++ simplex/axpy_simd.go | 6 +++--- 2 files changed, 24 insertions(+), 3 deletions(-) create mode 100644 simplex/axpy_archsimd_arm64.go diff --git a/simplex/axpy_archsimd_arm64.go b/simplex/axpy_archsimd_arm64.go new file mode 100644 index 0000000..de617ec --- /dev/null +++ b/simplex/axpy_archsimd_arm64.go @@ -0,0 +1,21 @@ +//go:build goexperiment.simd && arm64 + +package simplex + +import "simd/archsimd" + +// axpy computes dst[k] += factor*src[k], using NEON directly (no portable +// dispatch stub) via simd/archsimd's Float64x2 (128-bit, 2 lanes). +func axpy(dst, src []float64, factor float64) { + fv := archsimd.BroadcastFloat64x2(factor) + n := len(dst) + k := 0 + for ; k+2 <= n; k += 2 { + d := archsimd.LoadFloat64x2(dst[k : k+2]) + s := archsimd.LoadFloat64x2(src[k : k+2]) + s.MulAdd(fv, d).Store(dst[k : k+2]) + } + for ; k < n; k++ { + dst[k] += factor * src[k] + } +} diff --git a/simplex/axpy_simd.go b/simplex/axpy_simd.go index 1544ac7..3c23110 100644 --- a/simplex/axpy_simd.go +++ b/simplex/axpy_simd.go @@ -1,11 +1,11 @@ -//go:build goexperiment.simd +//go:build goexperiment.simd && !arm64 package simplex import "simd" -// axpy computes dst[k] += factor*src[k] for all k, vectorized via the -// experimental simd package (Go 1.27+, GOEXPERIMENT=simd). +// axpy computes dst[k] += factor*src[k], vectorized via the portable simd +// package. arm64 uses simd/archsimd directly (axpy_archsimd_arm64.go). func axpy(dst, src []float64, factor float64) { fv := simd.BroadcastFloat64s(factor) width := fv.Len() From 55e511bc936f8b4262f067c00b873b37d4e53a28 Mon Sep 17 00:00:00 2001 From: andig Date: Fri, 3 Jul 2026 13:37:37 +0200 Subject: [PATCH 3/4] simplex: reuse pivot-loop scratch buffers instead of allocating per pivot MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Allocation profile (300x300 random LP) showed phaseCost, duals, and alpha each allocating a fresh slice on every single pivot iteration: 6858 allocs/op, 21.2MB/op. None of that data needs to survive past the iteration, so run() now allocates the three scratch buffers once and clears+reuses them across pivots. Result: 916 allocs/op, 927KB/op (7.5x/23x less). Timing-wise this is a real but modest win, independent of the SIMD work in this branch — applies to every build, no GOEXPERIMENT flag needed (benchstat, n=10): ColdSolve_50x50 -5.36% (p=0.000) ColdSolve_150x150 -0.92% (p=0.007) ColdSolve_300x300 -0.24% (p=0.002) geomean -2.20% Co-Authored-By: Claude Sonnet 5 --- simplex/simplex.go | 45 +++++++++++++++++++++++++-------------------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/simplex/simplex.go b/simplex/simplex.go index bb18b00..39c7ea0 100644 --- a/simplex/simplex.go +++ b/simplex/simplex.go @@ -198,32 +198,30 @@ func (lp *LP) recomputeBasics(st *State) { } } -func (lp *LP) alpha(st *State, j int) []float64 { +// alpha fills dst (length lp.m, assumed zeroed) with column j's entries +// against the current basis: Binv * column(j). +func (lp *LP) alpha(st *State, j int, dst []float64) { rows, vals := lp.column(j) - a := make([]float64, lp.m) for k, r := range rows { c := vals[k] if c == 0 { continue } for i := 0; i < lp.m; i++ { - a[i] += st.binv[i][r] * c + dst[i] += st.binv[i][r] * c } } - return a } -// duals computes y = cost_B^T * Binv (as a length-m row vector). -func duals(st *State, cost []float64, m int) []float64 { - y := make([]float64, m) +// duals fills dst (length m, assumed zeroed) with y = cost_B^T * Binv. +func duals(st *State, cost []float64, m int, dst []float64) { for i := 0; i < m; i++ { cb := cost[st.basicOf[i]] if cb == 0 { continue } - axpy(y, st.binv[i], cb) + axpy(dst, st.binv[i], cb) } - return y } func (lp *LP) reducedCost(y []float64, cost []float64, j int) float64 { @@ -273,38 +271,43 @@ func (lp *LP) SetBound(j int, lb, ub float64) { lp.lb[j] = lb; lp.ub[j] = ub } // NumCols returns the number of structural columns. func (lp *LP) NumCols() int { return lp.n } -// phaseCost recomputes the Phase-1 cost vector: -1 for a basic var below -// its lower bound, +1 above its upper bound; false once all are feasible. -func (lp *LP) phaseCost(st *State) ([]float64, bool) { - cost := make([]float64, lp.nTotal()) +// phaseCost fills dst (length nTotal, zeroed) with the Phase-1 cost +// vector; returns false once all basic variables are feasible. +func (lp *LP) phaseCost(st *State, dst []float64) bool { inPhase1 := false for i := 0; i < lp.m; i++ { bv := st.basicOf[i] v := st.value[bv] switch { case v < lp.lb[bv]-eps: - cost[bv] = -1 + dst[bv] = -1 inPhase1 = true case v > lp.ub[bv]+eps: - cost[bv] = 1 + dst[bv] = 1 inPhase1 = true } } - return cost, inPhase1 + return inPhase1 } // run recomputes the active cost vector fresh before every pivot, so a // variable that becomes feasible mid-sequence stops influencing pricing. func (lp *LP) run(st *State) Status { + phase1Cost := make([]float64, lp.nTotal()) + y := make([]float64, lp.m) + a := make([]float64, lp.m) for iter := 0; ; iter++ { if iter > maxIter { return IterLimit } - cost, inPhase1 := lp.phaseCost(st) + clear(phase1Cost) + inPhase1 := lp.phaseCost(st, phase1Cost) + cost := phase1Cost if !inPhase1 { cost = lp.cost } - y := duals(st, cost, lp.m) + clear(y) + duals(st, cost, lp.m, y) q, dir := lp.chooseEntering(st, y, cost) if q < 0 { if inPhase1 { @@ -312,7 +315,8 @@ func (lp *LP) run(st *State) Status { } return Optimal } - a := lp.alpha(st, q) + clear(a) + lp.alpha(st, q, a) t, row, isFlip := lp.ratioTest(st, a, q, dir, inPhase1) if row < 0 && !isFlip { if inPhase1 { @@ -484,7 +488,8 @@ func (lp *LP) Solution(st *State) (x, rowActivity, reducedCost, rowPrice []float for i := 0; i < lp.m; i++ { rowActivity[i] = -st.value[lp.n+i] * -1 // logical var value == row activity (y_i = Ax_i) } - y := duals(st, lp.cost, lp.m) + y := make([]float64, lp.m) + duals(st, lp.cost, lp.m, y) rowPrice = make([]float64, lp.m) for i := range rowPrice { rowPrice[i] = y[i] * lp.objSign From aa8388970fc8b80da7258f7575bcbcfa3aed8750 Mon Sep 17 00:00:00 2001 From: andig Date: Fri, 3 Jul 2026 14:57:50 +0200 Subject: [PATCH 4/4] wip --- simplex/axpy_archsimd_amd64.go | 34 ++++++++++++++++++++++++++++++++++ simplex/axpy_archsimd_arm64.go | 13 +++++++++++++ simplex/axpy_generic.go | 7 +++++++ simplex/axpy_simd.go | 19 +++++++++++++++++-- 4 files changed, 71 insertions(+), 2 deletions(-) create mode 100644 simplex/axpy_archsimd_amd64.go diff --git a/simplex/axpy_archsimd_amd64.go b/simplex/axpy_archsimd_amd64.go new file mode 100644 index 0000000..6bd64e5 --- /dev/null +++ b/simplex/axpy_archsimd_amd64.go @@ -0,0 +1,34 @@ +//go:build goexperiment.simd && amd64 + +package simplex + +import "simd/archsimd" + +// axpy computes dst[k] += factor*src[k], using SSE2 directly (no portable +// dispatch stub) via simd/archsimd's Float64x2 — part of the amd64 baseline. +func axpy(dst, src []float64, factor float64) { + fv := archsimd.BroadcastFloat64x2(factor) + n := len(dst) + k := 0 + for ; k+2 <= n; k += 2 { + d := archsimd.LoadFloat64x2(dst[k : k+2]) + s := archsimd.LoadFloat64x2(src[k : k+2]) + s.MulAdd(fv, d).Store(dst[k : k+2]) + } + for ; k < n; k++ { + dst[k] += factor * src[k] + } +} + +// scale computes dst[k] /= divisor for all k, using SSE2 directly. +func scale(dst []float64, divisor float64) { + dv := archsimd.BroadcastFloat64x2(divisor) + n := len(dst) + k := 0 + for ; k+2 <= n; k += 2 { + archsimd.LoadFloat64x2(dst[k : k+2]).Div(dv).Store(dst[k : k+2]) + } + for ; k < n; k++ { + dst[k] /= divisor + } +} diff --git a/simplex/axpy_archsimd_arm64.go b/simplex/axpy_archsimd_arm64.go index de617ec..cb22ea6 100644 --- a/simplex/axpy_archsimd_arm64.go +++ b/simplex/axpy_archsimd_arm64.go @@ -19,3 +19,16 @@ func axpy(dst, src []float64, factor float64) { dst[k] += factor * src[k] } } + +// scale computes dst[k] /= divisor for all k, using NEON directly. +func scale(dst []float64, divisor float64) { + dv := archsimd.BroadcastFloat64x2(divisor) + n := len(dst) + k := 0 + for ; k+2 <= n; k += 2 { + archsimd.LoadFloat64x2(dst[k : k+2]).Div(dv).Store(dst[k : k+2]) + } + for ; k < n; k++ { + dst[k] /= divisor + } +} diff --git a/simplex/axpy_generic.go b/simplex/axpy_generic.go index 73cefef..747e2a8 100644 --- a/simplex/axpy_generic.go +++ b/simplex/axpy_generic.go @@ -9,3 +9,10 @@ func axpy(dst, src []float64, factor float64) { dst[k] += factor * src[k] } } + +// scale computes dst[k] /= divisor for all k. +func scale(dst []float64, divisor float64) { + for k := range dst { + dst[k] /= divisor + } +} diff --git a/simplex/axpy_simd.go b/simplex/axpy_simd.go index 3c23110..0a809ba 100644 --- a/simplex/axpy_simd.go +++ b/simplex/axpy_simd.go @@ -1,11 +1,11 @@ -//go:build goexperiment.simd && !arm64 +//go:build goexperiment.simd && !arm64 && !amd64 package simplex import "simd" // axpy computes dst[k] += factor*src[k], vectorized via the portable simd -// package. arm64 uses simd/archsimd directly (axpy_archsimd_arm64.go). +// package. arm64 and amd64 use simd/archsimd directly (axpy_archsimd_*.go). func axpy(dst, src []float64, factor float64) { fv := simd.BroadcastFloat64s(factor) width := fv.Len() @@ -20,3 +20,18 @@ func axpy(dst, src []float64, factor float64) { dst[k] += factor * src[k] } } + +// scale computes dst[k] /= divisor for all k, vectorized via the portable +// simd package. +func scale(dst []float64, divisor float64) { + dv := simd.BroadcastFloat64s(divisor) + width := dv.Len() + n := len(dst) + k := 0 + for ; k+width <= n; k += width { + simd.LoadFloat64s(dst[k : k+width]).Div(dv).Store(dst[k : k+width]) + } + for ; k < n; k++ { + dst[k] /= divisor + } +}