PERF: Use a plain loop in _pivoting! for small tableaus

[oyamad]

This PR speeds up the pivoting step shared by lcp_lemke and (via GameTheory.jl) lemke_howson, following the performance review of pivoting.jl/lcp_lemke.jl (#368/#371).

Changes

• _pivoting! dispatches between two kernels on total tableau size (PIVOTING_BLAS_CUTOFF = 2^14 elements; revised from a row-count criterion after review — the crossover tracks the memory traffic of the update, which row count misrepresents for the rectangular tableaus that arise in GameTheory.jl): at most 2^14 elements, a single fused SIMD loop in column-major order (dividing directly for non-BLAS eltypes such as Float16, whose reciprocals can overflow); above, the existing BLAS.scal! + BLAS.ger! path, unchanged. The function signature is unchanged, so lcp_lemke and GameTheory.jl's lemke_howson (which imports _pivoting! and passes col_buf) need no changes and inherit the speedup. Non-BLAS eltypes (e.g. BigFloat) now take the loop path instead of the broadcast fallback.
• A benchmark/lcp_lemke.jl suite is added, with end-to-end cases on both sides of the kernel dispatch and a preallocated lcp_lemke! case for the repeated-solve regime.
• New tests: the two kernels agree up to roundoff, _pivoting! selects the intended kernel on both sides of the cutoff, and the pivot column is reduced to the exact unit vector.

Rationale for the cutoff

A rank-1 update has no data reuse for BLAS kernels to exploit, so a plain SIMD loop matches single-threaded BLAS throughput exactly (measured: 268 μs for both at n=1000); the only structural advantage of BLAS is multithreading, which pays only where the work amortizes the thread-pool synchronization latency. The crossover size depends on the machine only through a square root of the hardware constants (sync latency, cache-level bandwidth, bandwidth scaling across cores), so a fixed constant is defensible. The crossover tracks the total tableau size (the memory traffic of the update) rather than the row count — confirmed by rectangular measurements: BLAS is 2.2x faster at (64, 1000) and 3.0x at (64, 10000), shapes that a row-count criterion would route to the loop. 2^14 elements errs low, which is the cheap direction: using the loop slightly below the true crossover costs a few percent, while the BLAS side is byte-identical to the current implementation.

Kernel crossover measured on three platforms (loop/BLAS time ratio; < 1 means the loop is faster):

n	macOS arm64 (M-series)	Linux x86-64 (EPYC 7763)	Windows x86-64 (EPYC 9V74)
16	0.25	0.54	0.25
48	0.72	0.67	0.43
64	0.72	0.66	0.07
96	1.40	1.41	0.13
128	2.02	1.21	0.15
192	3.12	1.80	0.37

macOS and Linux agree that the crossover lies between 64 and 96 rows on these square (n, 2n+2) tableaus, i.e. between ~8k and ~19k elements, bracketing the 2^14 cutoff. On the Windows CI runner, OpenBLAS's threaded ger! has pathological synchronization overhead (15 μs at n=64 against 1 μs for the loop) and the loop wins through at least n=256; below the cutoff this PR shields that platform from the pathology, above it the behavior is the status quo. An OS-conditional cutoff is left as a possible follow-up pending measurements on non-virtualized Windows hardware.

Benchmarks

lcp_lemke on random positive definite LCPs (minimum times, this machine): 40% faster at n=10 (0.9 → 0.54 μs), 9% at n=50, unchanged at n ≥ 96 (BLAS path). New suite baseline: dense_n10 0.40 μs, dense_n10_prealloc 0.35 μs (12 allocations — removing these is planned as a follow-up PR on lcp_lemke! workspace), dense_n50 16.9 μs, dense_n200 456 μs.

🤖 Generated with Claude Code (Claude Fable 5)

[Copilot]

Pull request overview

This PR improves performance of the pivoting step used by Lemke-style algorithms by introducing a fast plain-loop kernel for small tableaus while preserving the existing BLAS-based update for larger ones.

Changes:

• Added PIVOTING_BLAS_CUTOFF = 64 and dispatch in _pivoting! between a new SIMD loop kernel and the existing BLAS (scal! + ger!) kernel.
• Added tests ensuring the two kernels agree up to roundoff, _pivoting! selects the expected kernel across the cutoff, and the pivot column is exactly reduced to a unit vector.
• Added a new benchmark/lcp_lemke.jl suite and hooked it into benchmark/benchmarks.jl, with README documentation for the new suite.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated no comments.

Show a summary per file

File	Description
src/pivoting.jl	Adds loop-vs-BLAS kernel dispatch for _pivoting! and implements the new loop kernel.
test/test_pivoting.jl	Adds coverage comparing loop and BLAS kernels and verifying kernel selection and exact pivot-column cleanup.
benchmark/lcp_lemke.jl	Introduces an end-to-end benchmark suite for lcp_lemke / lcp_lemke! around the pivoting cutoff.
benchmark/benchmarks.jl	Registers the new lcp_lemke benchmark group in the top-level suite.
benchmark/README.md	Documents the new lcp_lemke benchmark group and its cases.

[oyamad]

Applied the findings of an external review by ChatGPT, in the two commits above:

FIX: direct division for non-BLAS eltypes (97f1f59). The reviewer found that the loop kernel's reciprocal normalization overflows for eltypes whose reciprocal range the pivot tolerances do not protect, and supplied a Float16 repro, which I confirmed on the previous head ([1.0 Inf NaN; 0.0 -Inf NaN]): inv overflows Float16 for pivots below 1/floatmax(Float16) ≈ 1.5e-5, and 0 * Inf then produces NaN where the rdiv!-based fallback that the kernel replaced stayed finite. The loop kernel now multiplies by the reciprocal only for BlasFloat eltypes (where any pivot passing the tolerance checks is orders of magnitude above the overflow threshold, and where this matches the BLAS kernel) and divides directly otherwise; the type-level branch folds at compile time. The reviewer's repro is added as a regression test.

PERF: dispatch on total tableau size (8411475). The reviewer questioned the row-count dispatch criterion for rectangular tableaus, such as the (n, m+n+1) Lemke-Howson tableaus of GameTheory.jl. Measurements confirmed the criterion misclassifies wide shapes — BLAS is 2.2x faster at (64, 1000) and 3.0x at (64, 10000), both previously routed to the loop — and that the crossover tracks the total tableau size (the memory traffic of the update), as the theory says it should. Dispatch is now on length(tableau) > 2^14, which classifies every measured point correctly: all six rectangular shapes on this machine plus the square series here and on the Linux CI runner.

Minor points, all applied: the dispatch boundary is tested exactly at 2^14 and 2^14 + 1 elements (with rectangular shapes), the benchmark fixtures assert num_iter > 0 so the suite cannot silently benchmark the trivial q >= 0 path, and lcp_lemke.jl is listed in the benchmark README overview.

Validation: test_pivoting.jl (13/13, including the new Float16 and boundary tests), test_lcp_lemke.jl (44/44), the benchmark suite builds, the four manual validation scenarios pass, and the full test suite is green.

🤖 Generated with Claude Code (Claude Fable 5)