perf(solvers): davidson-native few-k thread split — measured non-win, shipped as opt-in

[wladerer]

davidson-native for few-k large-nb — MEASURED NEGATIVE, shipped as opt-in

Campaign: exploit the cores that the native Davidson kernel strands when
nk < nthreads (e.g. Si-64 at nk=8 on a 16/22-core box). The kernel runs
#pragma omp parallel over k with openblas_set_num_threads(1), so at nk=8 it
uses 8 cores and leaves the rest idle while each k's dense subspace algebra
(QR / Rayleigh-Ritz GEMMs) runs single-threaded. Target ceiling was ~1.7×
(the ~43% dense-algebra share of the Si-64 SCF wall).

Verdict: the lever is a non-win. The per-k dense algebra is
memory-bandwidth-bound, so the k-parallel path already saturates the memory
system; handing the stranded cores to OpenBLAS does not accelerate it, and
splitting k into fewer parallel waves to give each more BLAS threads makes it
monotonically worse. Shipped anyway as an explicit opt-in — auto always
resolves to kpar, so no solve regresses relative to the pre-split kernel.

What changed

• C kernel (davidson_native.c): the three k-loops take an outer_threads /
  blas_threads split. if(outer_threads>1) means outer=1 runs the k loop
  serially in the master thread, keeping the multithreaded BLAS a plain level-1
  OMP region (no nesting). kpar (outer=nthreads, blas=1) is byte-identical to
  the old single-argument kernel.
• Adapter (native_davidson.py): _resolve_thread_split + the
  GRADWAVE_NATIVE_KMODE (auto|kpar|fewk) / GRADWAVE_NATIVE_FEWK_OUTER gate.
  auto → kpar (measured-safe default). fewk is opt-in for A/B and future
  higher-bandwidth hardware.
• Unit tests: mode-gate resolution + a fewk-vs-eager SCF correctness test.
• Docs: env knobs + the measured verdict in the module docstring.

Measurements (asus, all under the bench-lock protocol)

OpenBLAS scaling at Si-64 shapes (USE_OPENMP build, m≈15k, dim=308):

threads	zgemm RR-build	QR (geqrf+ungqr)
1	184 ms	673 ms
4	60 ms	288 ms
8	58 ms	290 ms
11	60 ms	297 ms
22	61 ms	533 ms

OpenBLAS plateaus at ~4 threads and QR regresses past ~11 — there is no
single-call speedup to harvest from the idle cores.

Nested OMP dense-wall microbench (nk=8, one iteration-equivalent):

outer×blas (=cores)	dense wall
8×1 (=8)	2608 ms (best)
8×2 (=16)	2772 ms
4×4 (=16)	3085 ms
2×11 (=22)	5204 ms
1×22 (=22)	16547 ms

The 8-way k-parallel kpar path (8×1) is already optimal; every reallocation
of cores toward BLAS is neutral-to-catastrophic — the classic signature of a
bandwidth-bound workload.

Isolated real-Hamiltonian solve — Si-64, nk=8, nb=154 (16 threads):
captured mid-SCF from the actual BatchedHamiltonian. Run at ecut=15 Ry
(m=8480) because the box has only 14 GB RAM and the 30 Ry config (m≈24k)
swaps at ~5 min/iter — see the hardware note below.

arm	outer	blas	ms/solve	niter	napply	max|ΔEv| vs kpar
kpar	16	1	44507	29	18495	0
fewk-o8	8	2	44062	29	18495	0
fewk-o4	4	4	47079	29	18495	0
fewk-o2	2	8	71820	29	18495	0
fewk-o1	1	16	≫ kpar (aborted — all 16 threads on a QR that already regresses past 11)

Monotonic: giving each k more BLAS threads (fewer k in parallel) is strictly
slower. kpar wins. Eigenvalues bit-identical across every split (ΔEv = 0),
so the split is a pure threading change — no numerical effect.

QE reference for the same Si-64 box/pseudos: 87 s full SCF.

Regression guard

auto resolves to kpar for every nk (few-k and many-k alike), so the
al-4^3-mesh (nk=36) and al-8^3-mesh (nk≥60) paths are byte-identical to the
pre-split kernel: _resolve_thread_split(36,16) == (16,1) etc. — covered by the
new test_resolve_thread_split_gate.

Engagement verification (mandatory)

Every timed arm wrapped registry.get("davidson-native") in a counting spy and
asserted fallbacks == 0 (via diagnostics["fallback_reason"]) and the actual
(outer_threads, blas_threads) used per solve — the split shown in each row is
the kernel's own report, not the requested value. GRADWAVE_EIGENSOLVER pinned
in the shell (read call-time since #469). Native engaged on every arm, zero
eager fallbacks. Smoke-verified separately on al4 (nk=36).

What remains of the 1.7× ceiling

Nothing reachable by adding cores. The 43% dense-algebra share is
memory-bandwidth-bound (and the few-k large-nb regime is memory-capacity-
bound on a 14 GB box). The ceiling would only move by cutting memory traffic —
lower-precision subspace storage, or a communication-avoiding/blocked QR that
touches the tall-skinny factors fewer times — not by threading. The fewk
lever is retained for hardware with materially more memory bandwidth (H100 HBM,
or a many-channel server), where the OpenBLAS plateau may sit higher.

Hardware caveat

The full-SCF Si-64 arm (eager baseline ≈2258 s from the campaign premise, and
current-native full SCF) could not be timed cleanly: asus has 14 GB RAM and the
30 Ry Si-64 working set (~10 GB native, more eager) swaps hard (~5 min/iter).
The isolated-solve table is the direct measurement of the lever; the per-solve
kpar number is current-native's real behaviour on the true Hamiltonian.

🤖 Generated with Claude Code

https://claude.ai/code/session_018j6BBnPtBSjkrnASSc52mT

[wladerer]

Parked per the no-complexity-for-no-wins rule: the few-k thread split is a measured monotonic non-win (kpar 16×1 = 44507 ms/solve; 8×2 tie; 4×4 0.95×; 2×8 0.62×; 1×16 aborted-worse; eigenvalues bit-identical, zero fallbacks). Mechanism: the Si-64 dense-subspace wall is memory-bandwidth-bound, not core-starved — nixpkgs OpenBLAS TSQR regresses past ~11 threads and kpar already ran at ~3.4 effective cores of 8 active. Branch perf/native-largenb kept as the ready-made A/B harness if high-bandwidth hardware (H100 HBM) materializes. Full record in the session memory.