simplex: experimental Go 1.27 SIMD axpy for pivot/duals (arm64: -18% via archsimd)

[andig]

Summary

• Profiling (pprof on a 300x300 random LP, forced to actually pivot) showed 76% of solver CPU time in two dense O(m²) loops: pivot's Binv row update (51%) and duals' accumulate (25%). Both are the same AXPY shape: dst[k] += factor*src[k].
• Extracted both into a shared axpy(), with build-tag-selected implementations:
  • axpy_generic.go (default, all toolchains): original scalar loop, byte-for-byte same behavior.
  • axpy_archsimd_arm64.go / axpy_archsimd_amd64.go (//go:build goexperiment.simd && arm64|amd64): vectorized directly via simd/archsimd's Float64x2 (NEON/SSE2), no portable dispatch overhead.
  • axpy_simd.go (//go:build goexperiment.simd && !arm64 && !amd64): portable simd package fallback for other SIMD-capable architectures.
  • All SIMD variants require Go 1.27rc1+ with GOEXPERIMENT=simd — 1.27 is not yet released.
• Default build/test is completely unaffected — the SIMD files only compile when the experiment flag is set, so this costs nothing to land as-is.

Benchmark results

benchstat, n=10 runs each, Apple M4 Pro (arm64 NEON, 128-bit = 2 float64 lanes), this branch vs. current main:

Scalar only (Go 1.26.4, default build, no SIMD — isolates the effect of the axpy() extraction itself):

ColdSolve_50x50-12      ~ (noise, no significant difference)
ColdSolve_150x150-12   -9.53%  (p=0.000)  faster
ColdSolve_300x300-12  -11.17%  (p=0.000)  faster

The extraction alone is a free win: a for k := range dst loop lets the compiler eliminate bounds checks that the original for k := 0; k < m; k++ loop couldn't prove statically.

With Go 1.27rc1 + GOEXPERIMENT=simd (direct archsimd NEON path):

ColdSolve_50x50-12      ~ (noise; per-call overhead dominates on 2-wide vectors)
ColdSolve_150x150-12  -35.49%  (p=0.000)  faster
ColdSolve_300x300-12  -31.05%  (p=0.000)  faster

Verdict: worth landing. Switching from the portable simd dispatch to archsimd directly (dropping the dispatch-stub overhead) turned this from a marginal, mixed-size-dependent win into a consistent 30%+ speedup at realistic problem sizes, with the scalar path also getting a smaller unconditional improvement from the refactor itself. It still regresses (within noise) at very small m, where call/bounds-check overhead on 2-wide vectors outweighs the parallelism — but that's fully mitigated by keeping this behind the GOEXPERIMENT=simd build tag until Go 1.27 ships.

🤖 Generated with Claude Code

[andig]

Verified this isn't a no-op or overhead artifact: disassembled the compiled `axpy@simd128` (`go tool objdump`) and confirmed real NEON instructions are emitted — `VDUP` (broadcast), 128-bit `FMOVQ` load/store, and a genuine vector `VFMLA` (fused multiply-add across 2 float64 lanes). `simd.Emulated()` also reports `false` on this hardware.

So the mixed benchmark result is real, not noise: NEON here is only 128-bit (2 float64 lanes), so the theoretical ceiling is 2x, and there's a runtime CPU-feature dispatch stub on every call (checked via a cached flag, branching to the vector body or a scalar fallback). That fixed per-call overhead outweighs the 2-lane gain at small m (regression at 50x50) and only pays off once the loop is long enough to amortize it (wins at 150x150/300x300). Consistent with the numbers, not an artifact.

[andig]

Split the allocation-reuse fix out into a standalone PR against main, since it doesn't depend on the SIMD work: #2. Once that merges, I'll rebase this branch — the third commit here (55e511b) will become redundant/need dropping at that point.

[andig]

@copilot resolve the merge conflicts in this pull request

[Copilot]

@copilot resolve the merge conflicts in this pull request

Done. Merged origin/main into this branch — the only conflict was in simplex/simplex.go in the duals function (inline loop vs. axpy() call); resolved by keeping the axpy() call. All tests pass.

[andig]

@copilot resolve the merge conflicts in this pull request

[Copilot]

@copilot resolve the merge conflicts in this pull request

Done — merged origin/main (through #5's binv flattening). Both conflicts were axpy() call sites in duals and pivot; resolved by updating to st.binv[i*m:i*m+m] slice notation to match the new flat layout. All tests pass.