Status: current · RTX 4090 (sm_89), CUDA 13.2 · checked 2026-10-06
Routing, tuning constants and measurements for the dense symmetric/Hermitian eigensolver syev: which provider Auto picks per scalar type, n and jobz, and the evidence behind each rule. The measurements date from 2026-08-03 to 2026-08-08. Machine: device 1 only, cuSOLVER 12.2.0.1, driver 595.84, unless a table says otherwise. Times are µs per matrix unless a table says ms.
The families live in src/ops/syev/choice.hh: cta, cta_fused and jacobi (the three n ≤ 32 sub-group solvers), blocked, two_stage and vendor. select::choose takes the first runnable entry of the nearest row of tuned/syev.<arch>.txt (arch sm_89 or sm_120), keyed on jobz, n and batch. can_run is correctness only: the native families need a non-NETLIB GPU backend and a square A with n ≥ 1, and the small three also need n ≤ 32 and a 32-wide sub-group. Every family, vendor included, needs a homogeneous batch (syev.cc:52-62). Uplo is not a key, because Upper is mirrored into Lower.
Note: the sm_89 syev table is transcribed at 424a45bc: it replays the previous router's order per grid cell, identical across batch, and the rule table below reads its rows. The sm_120 table is deep-measured (1226 rows, plus 260 transcribed rows where the lattice measured nothing near), and batch matters there (tuned/README.md, flat kernel selection).
| decision | where | rule |
|---|---|---|
| eigenvectors, n > 32 | jobz=V rows | float: blocked to 448, two-stage 449..1024, vendor above. double: blocked to 448, then vendor. complex<float>: blocked to 512, then vendor. complex<double>: blocked to 256, then vendor |
| eigenvalues only, n > 32 | jobz=N rows | blocked to 320, two-stage above, every type |
| n ≤ 32 family | rows with n ≤ 32, both modes | double: jacobi. float: jacobi to 8, cta_fused above. complex<float>: cta_fused to 8, cta above. complex<double>: cta |
| n ≤ 32 vendor handover | jobz=V, complex<double>, n = 25..32 | vendor first, cta second. Every other type stays native to 32 |
| backend scope | can_run | any GPU backend with a sub-group of 32; Cta and CtaFused not on ROCm (steqr_cta gives wrong eigenvalues on gfx1200). ROCm has no table and borrows the sm tables with a warning (untested) |
| values-only blocked solve | syev_blocked.stebz_evals | stebz bisection off the tridiagonal, no stedc |
| complex panel width | sytrd_block_size_default<T>, src/ops/syev/kernels/syev_blocked.cc | 32 for complex at 256 < n ≤ 512, harness value otherwise |
| complex trailing update | rank2k_trailing_update_supported, src/ops/syev/sytrd/sytrd_blocked.cc | float and complex<float> use the rank-2k path on the CUDA backend. double, complex<double> and other backends keep the GEMM pair |
| sb2st back-transform | sb2st_back_tile_for / sb2st_back_subs_for, src/ops/syev/sytrd/sytrd_sb2st_hh.cc | complex<float> at n = 512: tile 2, subs 4. Otherwise the shape-adaptive heuristic |
| block sizes | ORMQR_BLOCK_SIZE_*, SYTRD_BLOCK_SIZE_*, include/batchlas/tuning_params.hh | ORMQR 16/16/24/48/56, SYTRD 8/8/16/8/48 (buckets n ≤ 64/128/256/512/else) |
| sy2sb WY width | SY2SB_ORMQR_NB_LARGE = 32 | n ≤ 512 bucket only |
| stedc | STEDC_*, tuning_params.hh | leaf threshold 32, merge variant 2 (FusedCta), threads-per-root 8, wg multiplier 8 |
| stedc leaf cap | plan_stedc_levels, src/ops/syev/stedc/stedc_levels_plan.hh | leaf ≤ threshold, a hard cap |
| two-stage band width | choose_two_stage_kd, src/kernels/common/two_stage_common.hh | kd = 32 |
Pinning a family. BATCHLAS_SYEV_ROUTE=<family> (case-folded; also native, vendor, auto) or a select::ScopedPin. A pin that can_run refuses, or an unknown word, throws invalid_argument. The retired BATCHLAS_SYEV_PROVIDER, BATCHLAS_SYEV_SMALL_KERNEL and BATCHLAS_SYEV_CTA_MAX_N are no longer read.
The rule is keyed on n alone, never on batch: a tiny-batch caller pays launch overhead whatever is picked, and tuning for that regime costs saturated throughput. Auto dispatched to the winner at every shape of the grids below (Auto/expected 1.000 ± 0.002, idle GPU).
Float, eigenvectors, median of 3, build 12963a8; batch is the largest at which all three providers fit. µs per matrix.
| n | batch | blocked | vendor | two_stage | winner |
|---|---|---|---|---|---|
| 64 | 16384 | 1.64 | 2.02 | 2.21 | blocked 1.23× |
| 320 | 1302 | 74.53 | 215.23 | 111.20 | blocked 1.49× |
| 512 | 508 | 384.27 | 504.69 | 380.02 | two_stage 1.01× |
| 1024 | 254 | 4089.02 | 2706.84 | 2441.93 | two_stage 1.11× |
| 2048 | 64 | 33842.6 | 15019.1 | 24782.3 | vendor 1.65× |
Also n = 128, 256 (blocked, 1.18× and 1.01×) and 640, 768 (two_stage, 1.15× and 1.19×). The batch-keyed rule this replaced was wrong at five of nine sizes (1.11×–1.23×). For n ≤ 32 at batch 2048, cta/Jacobi (n = 16) and cta/fused (n = 32) beat every non-CTA provider by 2.8×–3.0×. Caveats: n = 2048 was measured at batch 64 (below the 128 SMs, memory-limited), so that boundary rests on the weakest row; the grid is float only.
Quoted on the documentation front page. 2026-08-07, float, eigenvectors, median of 5, default nb, fills the 320..512 gap. µs per matrix.
| n | batch | blocked | two_stage | vendor | winner |
|---|---|---|---|---|---|
| 320 | 819 | 67.84 | 113.39 | 203.00 | blocked (3.0× over vendor) |
| 448 | 585 | 195.27 | 262.79 | 400.59 | blocked (2.1× over vendor) |
| 512 | 512 | 326.71 | 332.55 | 504.40 | tie |
| 768 | 192 | 1241.90 | 1172.90 | 1348.30 | two_stage |
| 1024 | 128 | 3296.80 | 2614.70 | 2552.10 | tie (vendor/two-stage) |
(n = 384: blocked 120.99, two_stage 179.18, vendor 309.10.) The float crossover is 448, not 320 (the former rule lost 1.35×–1.48× at n = 384 and 448). Two-stage wins by at most 1.06× for float in eigenvector mode, so the 449..1024 band is neutral and was kept.
Same method, 2026-08-07; blocked at its own best nb (see the complex panel width). Two-stage is never the complex winner in eigenvector mode. µs per matrix (blocked / vendor):
Full grids: benchmarks/results/syev_retune_phaseA.jsonl.
Warning: the complex<double> boundary (256) is the most hardware-specific number on this page. The card runs FP64 at 1/64 rate, which penalises the native panel far more than cuSOLVER. Re-measure on a data-center GPU before trusting it. After the Annex G fix, complex<double> at n = 288..512 is a dead heat inside the 1.10× neutral band, so the boundary stays at 256.
Annex G fix (e4186908). Clang emitted the C99 Annex G complex multiply inside the panel symv (an isnan branch plus __mulsc3). Writing that multiply-add in real arithmetic made the panel kernel 1.22×–1.29× faster (real types unchanged: mac routes them to a plain multiply-add). At n = 512, batch 512 the complex/float cost ratio is 2.44–2.93 for blocked and two-stage against 1.38–1.41 for cuSOLVER; the rest of the gap is sytrd_blocked.panel_only (68% of the complex solve), see the open panel symv. Old Auto to new Auto (73eb5510), n = 320..1024: complex<float> 1.18×–1.69×, complex<double> 1.22×–2.88×, double 1.00×–1.97×, float 1.00×–1.50×; no cell regresses.
Float, 2026-08-04, same method (syev_benchmark gained a jobz argument for this; before that the mode had no vendor arm). µs per matrix.
| n | batch | blocked | vendor | two_stage | winner |
|---|---|---|---|---|---|
| 64 | 16384 | 1.06 | 1.12 | 1.14 | blocked 1.06× |
| 320 | 2604 | 48.16 | 200.94 | 50.35 | blocked 1.05× |
| 512 | 1017 | 298.30 | 458.95 | 158.72 | two_stage 1.88× |
| 1024 | 254 | 3537.41 | 2494.44 | 908.32 | two_stage 2.75× |
| 2048 | 64 | 29547.70 | 13908.80 | 10804.10 | two_stage 1.29× |
Blocked wins n = 64..320, two-stage 512..2048 (n = 128, 256, 640, 768 not shown). With no eigenvectors to apply, two-stage never pays for \(Q_2\). It is also the complex winner above n = 320 (1.30× over cuSOLVER at n ≥ 768), so there is no per-type split here.
These figures predate the stebz change (WP1), which made blocked 1.36×–1.87× faster up to n = 320; the 320 | 512 boundary has not been re-derived (see open debts). The rule replaced syev_prefer_two_stage_values, whose batch >= 256 term sent n = 1024 at batch 254 to the vendor (2.75× loss): a batch floor reopens that hole. At batch 1 the vendor wins by 12.5× at n = 1024 (batch starvation of both reductions).
The Lower-only blocked and two-stage pipelines get Upper input by mirroring the upper triangle into the lower one ( \(O(n^2)\) before an \(O(n^3)\) solve). Float, 2026-08-04, median of 3, n = 128..1024, both modes: Auto (Upper) beats vendor (Upper) at every shape (n = 320 vectors 2.95×, n = 1024 values 3.39×, n = 256 vectors 1.02×) and Upper/Lower is 1.003×–1.020×. A native Upper reduction could recover at most that 0.3–2.0%, so it was rejected.
mirror_upper_to_lower (src/ops/syev/uplo_mirror.hh, kernels/uplo_mirror.cc) writes \(A_{ji} := \overline{A_{ij}}\) in place (syev overwrites A); the diagonal is left alone and a complex diagonal's imaginary part is not forced to zero, as the Lower path already assumes. The header is declaration-only with explicit instantiations in uplo_mirror.cc (a SYCL kernel name needs one definition: inline in both syev_blocked.cc and syev_two_stage.cc gave "definition with same mangled name"). The obvious test is vacuous because Matrix::Random(..., symmetric=true) makes Upper and Lower interchangeable: the fixture poisons the strictly-lower triangle after taking the reference from the upper one and asserts a Lower-read solve differs by more than 1.0 ("fixture is vacuous" otherwise). Two cases in tests/syev_blocked_tests.cc cover all four types, one per call site.
The CTA-family kernels syev_cta, syev_cta_fused, syev_jacobi_cta are the families cta, cta_fused, jacobi (BATCHLAS_SYEV_ROUTE=cta|cta_fused|jacobi). The measurements below used the retired BATCHLAS_SYEV_SMALL_KERNEL knob.
Build 7911847. Median of 5 process-level repeats with IQR, each cell at its knee batch, best over cta_wg_size_multiplier. A margin below 1.10× or with overlapping IQRs is neutral. Key results on the CSV name column (each binary registers two benchmarks). µs per matrix.
| type | n | mode | winner | µs | runner-up | ratio | vs cuSOLVER |
|---|---|---|---|---|---|---|---|
| double | 4 | vectors | jacobi | 0.01072 | cta_fused | 3.05× | 15.17× |
| double | 32 | vectors | jacobi | 3.31753 | vendor | 1.37× | 1.37× |
| float | 8 | vectors | jacobi | 0.00478 | cta_fused | 1.37× | 8.42× |
| float | 16 | vectors | cta_fused | 0.02933 | jacobi | 1.22× | 4.10× |
| float | 16 | values | cta_fused | 0.02740 | jacobi | 1.10× (neutral) | - |
| float | 32 | vectors | cta_fused | 0.28426 | tridiag_ref | 1.11× | 1.79× |
Double (n = 4, 8, 16, 32, both modes): Jacobi wins every cell, 1.37×–3.59×. Full grid: benchmarks/results/syev_retune_phaseA.jsonl. syev_cta, where Auto sent every n ≤ 32, won no cell. The rule is the float and double n ≤ 32 rows of tuned/syev.sm_89.txt: double Jacobi everywhere; float Jacobi to n = 8, cta_fused above. The float rule rests on the vector cells; n = 9..15 was not swept. A later A/B found fused with its tuned multiplier beats Jacobi from n = 7; the 8 | 9 boundary stays for Jacobi's relative accuracy on graded input (4.5e-07 against 2.7e+28 for syev_cta, quoted from the uncommitted Jacobi plan, see the Jacobi design).
Warning: FP64 at 1/64 rate inflates Jacobi's margin over the tridiagonalising paths. The float column predicts a 1:2 FP64 datacenter GPU better; gate the double rule on measured FP64 throughput there. See the small-n steqr routing.
complex<float>, eigenvectors, 2026-08-07, median of 3 at each cell's saturating batch (µs per matrix):
| n | batch | cta | cta_fused | jacobi | vendor | winner |
|---|---|---|---|---|---|---|
| 8 | 32768 | 0.0132190 | 0.0105830 | 0.0143920 | 0.044342 | fused 1.25× |
| 9 | 29127 | 0.0213110 | 0.0281870 | 0.0544510 | 0.079653 | cta |
| 32 | 8192 | 0.3825800 | 0.5948600 | 2.4824000 | 0.470900 | cta |
(n = 4 jacobi/fused; n = 5..7 fused 1.21×–1.39×; n = 12, 20, 24 cta; n = 16 tie.) syev_cta is right for complex from n ≥ 9; Jacobi is 4×–6× off the pace at n ≥ 20, so the float rule would be wrong here. For complex<double>, fused leads cta by only 1.03×–1.08× at n = 4..16 (neutral), so that range is not split.
The complex<double> jobz=V rows of tuned/syev.sm_89.txt put vendor first for n = 25..32. 2026-08-07, median of 3, µs per matrix (cta / cta_fused / vendor): n = 20 3.0015 / 3.4597 / 3.1085 (cta); n = 24 4.0806 / 4.6301 / 3.8184 (tie); n = 26 5.4993 / 5.3048 / 4.4835 (vendor 1.23×); n = 32 8.0253 / 7.4929 / 5.6531 (vendor 1.42×). An FP64-rate artifact and machine-specific. complex<float> does not cross over (cta beats the vendor 1.23× at n = 32), so the rule is deliberately not "all complex".
syev_cta runs three launches (sytrd_cta, steqr_cta, ormqx_cta) plus two pack kernels and a copy on the eigenvector path, each stage round-tripping global memory. syev_cta_fused keeps one problem resident in one sub-group partition from load to store (one read of A, one write, against about seven round trips) with the same device code (sytrd_cta_device.hh, steqr_cta_device.hh), so a head-to-head measures fusion alone.
syev_jacobi_cta is partition-resident cyclic two-sided Jacobi: one SubGroupPartition<P> owns one problem, A and (optionally) Z stay in local memory, one launch per decomposition.
Results: the 2026-08-03 bake-off.
Historical; retired 2026-10-05. BATCHLAS_SYEV_CTA_MAX_N (0..32, default 32 = off) sent small eigenvector solves above the threshold to the vendor, motivated by n = 30, batch 8, float: cta 229.6 µs against cuSOLVER 103.7 (2.21×); nsys put 29.4% of LOBPCG GPU time (about 16% end to end) in that projected Rayleigh–Ritz solve. Through LOBPCG, eigenvectors, n = 256, µs/matrix, value 32 (off) / 16 / 8 / 0: batch 8 15563 / 14211 / 13746 / 13551; batch 64 1998.6 / 1890.4 / 1948.0 / 1945.9 (1.10×–1.15× and 1.03×–1.06×). It shipped off because moving the projected solve off CTA flips a near-tie in ILUKTests.SyevxInstrumentationAndPreconditioner (asserts lose_count == 0; at threshold 16 one case crosses to ratio 1.25 at an already-converged point, 4.2e-06). Open owner decision: tolerate a tie on an already-converged case, or pin that test's projected-solve provider.
Two files answer #include <batchlas/tuning_params.hh>: the committed include/batchlas/tuning_params.hh and one CMake generates into ${PROJECT_BINARY_DIR}/include/batchlas/. src/CMakeLists.txt lists the source include directory first, so the committed header wins on every build; the CMake defaults in cmake/BatchLASGeneratedHeaders.cmake and any profile-driven retune written to the binary directory are dead (2026-08-03: ormqr_block_size_for_n(512) = 16, sytrd_block_size_for_n(512) = 24; the generated header would say 128 and 64).
Warning: a retune takes effect only when its constants are ported by hand into the committed header (the 2026-08-07 retune, 924b3a59, did). Making the committed header the single source of truth is not done.
Every accessor is overridable through BATCHLAS_TUNE_*; overrides are read on every call and feed *_buffer_size as well as the solve, so they must not change between a buffer-size query and its call.
Hypothesis: shadowing pinned ORMQR at 16 where a generated header said 128, explaining why BatchLAS stopped beating cuSOLVER at n = 32..512. Refuted: (1) ORMQR has been 16 in every bucket since the header was created (7363746, 2026-03-02); what moved was SYTRD (91341c1, unmeasured). (2) A/B at n = 512, eigenvectors, batch 128, median of 5 interleaved rounds: ORMQR = 128 gave 452.32 µs/matrix (IQR 450.73–454.91) against 301.95 (IQR 300.79–303.75) unset, 1.50× slower. (3) 16 was the measured optimum wherever the knob mattered. No baseline exists to restore: git history holds no artifact of BatchLAS beating cuSOLVER at n = 32..512, so nothing may calibrate against that claim. The values came from a standalone ormqr_blocked microbenchmark whose n ≤ 64 search offered only block sizes 4–16 (16 was the ceiling), and the generator's _derive_param_buckets keeps the parameters of whichever case had the lower absolute time, not the best block size per case.
Superseded by the 2026-08-07 retune. First measurement of BATCHLAS_TUNE_ORMQR_BLOCK_SIZE and BATCHLAS_TUNE_SYTRD_BLOCK_SIZE through syev: blocked provider, float, knee batch capped at 512 by the workspace limit. ORMQR 16 was optimal at n = 32..128 (n = 128: 7.446 µs against 8.142 with ORMQR 32 and 22.805 with ORMQR 128); at n = 256 and 512 ORMQR 32 won by only 1.016× and 1.059× (neutral). SYTRD matched the committed 8/8/16/24 at n = 64/128/256/512 (n = 512: 24 gives 339.9, 8 gives 380.8). Double agreed at n ≤ 128 and preferred SYTRD 16 at n = 512 (1.05×–1.07×, neutral). ORMQR = 128 cost 2.1× at n = 64 and 3.1× at n = 128. Gap: n = 256 missing in both modes and precisions (syev_blocked exceeded the 3 GB per-invocation abort even at batch 476). Raw: benchmarks/results/syev_retune_phaseB_float.jsonl.
Commit 924b3a59 retuned on a trimmed space and A/B'd each group end to end at the consumers. CUDA/float, µs per matrix, before to after: syev n=64 b=4096 1.0798 to 0.9938 (1.087×); syev n=512 b=512 166.95 to 151.01 (1.106×); gesvd n=128 b=1024 9.2465 to 8.0684 (1.146×); gesvd n=256 b=512 1.032×; gesvd n=512 b=256 1.030×. Everything else was within 1%.
Tuning is float-only and the header is regenerated from a float-only bench, so per-type corrections live at the consumer where the next retune cannot overwrite them.
sytrd_block_size_default<T> (src/ops/syev/kernels/syev_blocked.cc) returns 32 for complex at 256 < n ≤ 512. The harness value 8 costs complex 1.16×–1.20×: a complex panel column moves twice the bytes and does four times the flops per element, so it needs a wider panel to amortise the O(n) barrier-and-reduction chain of latrd_lower_panel. complex<float>, eigenvectors, blocked, µs/matrix, median of 5, device 1, 2026-08-07:
| n / batch | nb=8 | nb=16 | nb=24 | nb=32 | nb=48 | nb=64 |
|---|---|---|---|---|---|---|
| 320 / 819 | 197.25 | 173.24 | 166.27 | 168.55 | 175.92 | 190.95 |
| 512 / 512 | 942.68 | 838.73 | 815.30 | 802.74 | 811.89 | 827.70 |
(n = 384 and 448 also have nb = 32 optimal.) nb = 32 is optimal or within 1.5% at every n in the bucket. Float at n = 512 spreads only 1.04× across nb 8–48 (inside the 1.10× neutral band), so float keeps the harness value. The override lives at the consumer, not in tuning_params.hh, because that header is regenerated from a float-only bench. Other complex buckets: n ≤ 128 prefers 8; 128 < n ≤ 256 prefers 16 (n = 192: 33.83 against 36.09 at 8); n > 512 is within 1.03× of the committed 48 at nb = 32.
The prediction that removing the split-WY penalty (f7f3c57, nb = kd) would move the eigenvector optimum to kd = 96–128 (MAGMA; Gates, Tomov and Dongarra 2018) is disproved. Total ms, eigenvectors, median of 3, 2026-08-04:
| n / batch | kd=16 | kd=32 | kd=64 | kd=128 | blocked |
|---|---|---|---|---|---|
| 256 / 1024 | 69.5 | 61.3 | 65.7 | 75.7 | 40.9 |
| 512 / 512 | 217.7 | 194.2 | 236.1 | - | 191.4 |
| 1024 / 128 | 445.3 | 369.9 | 433.2 | 715.3 | 475.5 |
kd = 32 is optimal at every n ≥ 256 (kd = 48 and 96 are worse throughout); wide kd wins only at n = 128 (kd = 128: 16.9 against 22.0 at kd = 32; blocked 14.6), where two-stage loses to blocked anyway. The nb hint helps narrow bands (1.02×–1.06× for kd ≤ 64) and hurts kd = 128 (0.93×) because LARFT work is \(O(m k\,n_b)\) (see sytrd: the dense-to-band block-width hint).
"Two-stage wins at n ≥ 1024" is shape-dependent. At n = 2048 / batch 32 (kd = 32: 1066 ms; kd = 128: 1943), legacy latrd gave blocked 1235.1 against two-stage 1066 (two-stage 1.16×); grid latrd (the default; 1.41× over legacy) gives blocked 875.8 (blocked 1.22×). At n = 1024 / batch 128 the batch saturates the 128 SMs and two-stage wins 1.29× either way.
Mode. Both modes reduce at a real band width and share choose_two_stage_kd (the Householder chase retains \(Q_2\); see sytrd: the Householder chase against the Givens chase). The literature puts the eigenvector optimum higher (96/128 on GPU); this table does not reproduce it.
latrd_grid_min_n = 768 was derived in eigenvalues-only mode and applies in both. legacy/grid ratio in eigenvector mode (> 1 means grid wins), batch 1 / 8 / 64: n = 512 1.012 / 1.018 / 0.956; n = 768 1.379 / 1.375 / 1.090; n = 1024 1.843 / 1.781 / -; n = 256 0.734 / 0.734 / 0.788 (n = 384: 0.881–0.894). The crossover is the same in both modes, so the constant is not split. The win shrinks with batch: once the batch alone saturates the SMs there is no starvation left to absorb.
The plan turned a profile of float and complex<float> syev at saturating batch into work packages with an end-to-end accept gate (a kernel win is not a solver win: a 2.16× kernel win became an 11% gesvd loss). Rig: build-cuda (sm_89, RelWithDebInfo), device 1, one process, idle GPU, --warmup=2 --min_iters=5; positional arguments n batch nb fuse jobz uplo. Rig validation at n = 512, batch 512, two-stage, eigenvectors: float 334.28 against the research profile's 332.5, cfloat 971.17 against 972.4; differences under about 1% are noise. Accuracy runs planned for WP1, WP3, WP8 are not recorded.
| WP | item | gate | measured | verdict |
|---|---|---|---|---|
| WP0 | harness unblockers | cfloat rows appear | rows appear | passed |
| WP1 | A2: stebz values-only | ≥ 1.20× float n=256 | 1.44× | passed |
| WP4 | A3: per-type tile/subs | ≥ 1.10× cfloat, float unchanged | 1.139×, float unchanged | passed |
| WP3 | A1: her2k trailing update | ≥ 1.05× cfloat | 1.043× | marginal miss, kept |
| WP2 | B3: cta-large-n | ≥ 1.15× vs blocked | 0.012× (85× slower) | rejected |
sb2st_hh_benchmark registered only float, so --type=cfloat silently produced zero rows (fixed; n = 512, batch 512, kd = 32: CHASE float 56.578, cfloat 267.70 ms, 4.73×; BACK 99.095 / 317.03, 3.20×). BATCHLAS_LATRD_GRID_GROUPS cannot force the grid latrd path: it is clamped to MAX_COMPUTE_UNITS / batch, so at batch ≥ 128 no value reaches the grid kernel. BATCHLAS_LATRD_GRID_FORCE_UNSAFE bypasses the cap but is deadlock-capable: the software grid barrier hangs, without failing, when the launch exceeds residency (it looks like slow JIT; run under a timeout). It was added and never exercised.
In values mode the blocked solver ran a full eigenvector divide-and-conquer and discarded Z (28.3% of the float solve at n = 256, batch 1024). stebz bisection (indices 0..n-1, ascending) now runs straight off the real tridiagonal, as two-stage already did (f28a3fcd). Values only, provider=blocked, batch 1024, µs per matrix (before to after): float n = 64 1.1614 to 0.62228 (1.87×), n = 256 20.774 to 14.476 (1.44×), n = 320 42.804 to 31.438 (1.36×); cfloat n = 64 1.5055 to 1.0465 (1.44×), n = 256 47.218 to 38.250 (1.23×), n = 320 101.45 to 87.666 (1.16×). The arithmetic bound was 1.39× against a gate of 1.20×. stedc runs in real arithmetic, so it is 28.3% of the float solve but about 12% of the cfloat solve at n = 256, hence float gains more. Values mode also stops allocating eigenvector-only workspace (n = 256, batch 1024, cfloat: z_span 268 MB, zc_span 537 MB, plus the stedc workspace).
provider=two_stage, eigenvectors, n = 512, batch 512, µs per matrix: cfloat shipped constants 971.17, forced tile=2 subs=4 852.77 (1.139× faster), per-type code 852.46; float shipped 334.28, forced tile=2 subs=4 384.30 (1.15× slower), per-type code 331.25 (unchanged). A global flip trades a 1.14× complex win for a 1.15× float loss; the per-type selection (8c24626f) keeps the complex win and leaves float at its optimum. The complex back-transform runs at SM throughput 35.3% and occupancy 49.8% (float 93.7% and 82.8%). Only n = 512 was measured, so the rule is n == 512 exactly. The instantiated set is tile ∈ {1,2,4,8} × subs ∈ {4,8,16}; any other pair silently falls through to the slower tiled kernel. WP4 and WP6 are not additive: measure any complex stage-2 estimate against 852, not 971.
For complex, sytrd_blocked's trailing update issued two full \(n_2 \times n_2\) GEMMs where float issues one triangle-only syr2k. her2k does half the arithmetic ( \(\alpha A B^H\) and \(\bar\alpha B A^H\) are conjugate transposes) via one batched GEMM into scratch plus a Hermitian fold (405b4120). Primitive, \(n_2\) = 480, batch 512, cfloat, k = 32: her2k 3.3010 ms against the GEMM pair 4.3866 (1.33×; 1.32× at k = 16, 24), nearly independent of k, so a larger nb will not grow the win. End to end, provider=blocked, eigenvectors, n = 256, batch 1024: cfloat 66.628 to 63.899 (1.043×); float unchanged. The gate was 1.05×; the miss was predicted and the change kept (real 1.04× at no runtime cost; it forced a single-definition route prediction in expansion_budget.hh). The gate should read 1.04×, not "met".
The unmerged cta-large-n branch lifted the CTA eigensolver above n = 32 (local-memory limits n = 128 float, 64 cfloat). Tests passed but it had no performance measurement. Forced provider=cta against blocked, float, eigenvectors, batch 256: n = 33 150.17 against 1.7715 µs (85× slower); n = 64 760.68 against 3.6021 (211× slower). At batch 2048 an n = 33..128 sweep finished no iteration in ten minutes while holding 24 GB. A CTA-resident solve saves about 15 launches but runs an unblocked level-2 algorithm in one work-group per matrix against a blocked level-3 pipeline. Reverted (14bad498; it survives on worktree-wf_3b4af334-426-9). The branch also rewrote the live n ≤ 32 CTA kernels and the sub-group partition (src/kernels/common/sg_partition/), never measured; if wanted it returns as its own change with n ≤ 32 numbers.
Note: benchmark through the forced provider. Under flat selection the cap is kSmallMaxN = 32 in the small families' can_run (src/ops/syev/syev.cc). A measured dead end in the agent guide.
Measured 2026-08-05. PR 55 (0bb92fb, the STEDC level-driver rework) cost syev 1.05×–3.25× at every size (the unchanged vendor column reproduced the grid, so the deltas are code). Float, eigenvectors, baseline / after PR 55 / fixed, µs per matrix: n = 64 1.64 / 2.25 / 1.70; n = 320 74.53 / 242.43 (3.25×) / 70.96; n = 640 727.48 / 1095.9 (1.51×) / 714.99; n = 1024 2441.93 / 2559.1 / 2394.7. Fixed is within 0.96×–1.05× of baseline at all eight measured sizes (the four not quoted: n = 128 1.01×, 256 1.03×, 512 1.02×, 768 1.02×).
Defect A: the merge variant was flipped the wrong way. 3072ea6 moved STEDC_MERGE_VARIANT_* from FusedCta (2) to Fused (1), claiming FusedCta was numerically wrong and 2–12% slower; neither reproduced. With variant 2 all 16 CUDA stedc_tests pass in both precisions, plus syev_tests 8/8, syev_blocked_tests 44/44 and syev_two_stage_tests 20/20; the NaN was a symptom of the deadlock fixed in 4bde59a. FusedCta is 11–35% faster on stedc on both drivers (n = 64: Levels+Fused 0.672, Levels+FusedCta 0.499; n = 512: 24.744 against 22.228).
Warning: FusedCtaConditionedHeavyDeflation asserts only finite-and-sorted, never accuracy. Strengthen it before trusting variant 2 on unmeasured hardware. Current merge constants are on the stedc page.
Defect B: the level planner chose a leaf past the STEQR cliff. plan_stedc_levels searched leaves in [threshold/2, threshold·2], tie-breaking toward the leaf nearest the threshold; at n = 320 and 640 it chose leaf 40 over 20. steqr takes the fast steqr_cta path only for n ≤ the sub-group width and falls back to steqr_wg above it. steqr alone, batch 10416, eigenvectors: n = 32 0.262 µs, n = 36 3.761, n = 40 4.870, n = 64 36.78. A 12.5% width increase costs 14×: a code-path cliff, not scaling. At n = 320 thresholds 16..24 chose leaf 20 (stedc 12.5–13.1) while threshold 32 chose leaf 40 (168.4; 649.7 at threshold 64). STEDC_RECURSION_THRESHOLD = 32 is the sub-group width, an invariant. The fix caps the leaf at the threshold (hi = threshold, stedc_levels_plan.hh); the level driver is then faster at these sizes (10.8 against 11.1 Recursive at n = 320). It escaped review because stedc_benchmark registered only n ∈ {64, 128, 256} (320 and 640 now are), and the plan shape is asserted directly in the host-only StedcLevelPlan.* tests, since a bad leaf still gives correct eigenvalues.
Open defect, measured 2026-08-03. syev_cta, n = 32, eigenvectors, peak device memory by batch: 512 445 MiB; 1024 2237 MiB (0.320 µs/matrix); 4096 7505 MiB (0.192); 16384 24083 MiB (whole card; 0.592, regresses). That is about 1.5–2.2 MiB of workspace per 32×32 matrix at batch 1024–16384 (1.8 MiB at 4096; roughly 380–560× the 4 KB of data), growing about linearly with batch above 512. Eigenvalues-only is unaffected (about 343 MiB at batch 16384). syev_blocked has the same problem: eigenvalues-only at batch 16384 (n = 32), 15258 (64), 3814 (128) and 953 (256) each exceeded 3 GB (predates WP1, not re-measured). Usable batch is capped near 4096 for n = 32 with vectors on a 24 GB card, below saturation; it caused a real OOM on a shared machine, binds every small-n eigenvector solve, is uninvestigated, and makes n = 2048 unsaturable.
Past the knee extra batch buys little (n = 8, float, syev_cta): batch 65536 gives 0.0075321 µs/matrix in 487 MiB; batch 1048576 gives 0.0071809 (4.7% better) in 24083 MiB (the run that caused the OOM). Saturation therefore means the knee, where a 4× batch step buys less than 5%. The bake-off driver enforces about 1.5 GB estimated footprint, a 3 GB hard abort per invocation and a 120 s timeout, and reports a cell that does not reach its knee within them.
Plan WP5 (B1), open. The "open 2.7×" in the agent guide. The latrd panel is 35–71% of every blocked solve. Counters on LatrdLowerPanel (ib = 32, j₀ = 0), as a multiple of ideal traffic for float n=256 b=1024 / float n=512 b=256 / cfloat n=512 b=256: DRAM 0.14× / 1.13× / 1.43×; L2 1.95× / 2.34× / 2.50×; L1TEX 11.9× / 12.2× / 15.8×; SM throughput 52.7% / 10.8% / 10.7%. ncu on LatrdLowerPanelKernelLegacy<float, 256, 0>, n = 512, batch 256 (2026-08-08): L1TEX 52,499 MB, SM 10.84%, occupancy 33.19%, 11.444 ms; grid (256,1,1), one work-group per matrix.
Mechanism. Thread r reads \(A(r,c)\) on its row walk and thread c reads the same element as \(A(c,r)\) on its column walk (2.3× L2); on the column walk consecutive lanes are lda apart (12–16× L1). The DRAM floor at n = 512 caps the gain at 2.7×. The kernel reaches only about 330 GB/s of roughly 1000 GB/s, so it is latency- or occupancy-bound.
Design constraint. The symv sits inside the per-reflector loop of one work-group per matrix, with v_local[n] and wcol_local[n] in local memory; reflectors are sequentially dependent, so a MAGMA-style multi-block symv is not a drop-in. Proposed: stage 32×32 off-diagonal tiles with a coalesced load, then compute \(y_R \mathrel{+}= A_{RC} v_C\) and \(y_C \mathrel{+}= A_{RC}^H v_R\) from that one copy (a cfloat tile is 8 KB; n = 512 cfloat needs 16 KB, still ≥ 2 blocks per SM). An earlier "one sub-group per column" attempt changed the access pattern without reducing reads. Predicted end-to-end gain at 2.7× on the panel: cfloat n=512 vectors 1.82×, float n=256 vectors 1.29×, float n=256 values 1.50× (1.56×, 1.22×, 1.36× at 2.0×). Gate: ≥ 1.6× on the panel kernel at n = 512, with ncu showing L1TEX falling from about 12× toward 2× and L2 from 2.3× toward 1×.
Plan WP6 (B2), open. At n = 512, batch 512 the float back-transform is saturated (SM 93.7%, occupancy 82.8%) and the complex one is not (35.3% / 49.8%); the chase runs 45.0% / 65.5% float and 29.8% / 45.1% complex. Halved occupancy is a register-pressure signature, so the 3.90× complex cost is not the price of complex arithmetic.
Prize: reaching the roughly 2.2× complex/float ratio the panel reached would cut the two stage-2 kernels from 1353 ms to 704 ms, about 1.49× on the complex two-stage solve (about 650 µs/matrix against 973.8), ahead of blocked (698) and the vendor (707) at n = 512. Any routing change waits for the measurement.
The values-mode paths of syev_blocked and syev_two_stage solve the tridiagonal with stebz bisection, which records no status: its three exits (tolerance met, budget exhausted, midpoint stopped advancing) look alike, so both tiers clear info to 0 for every item. stebz needs a convergence flag. The secular solver's budget-exhaustion arm also uses a literal the tuning parameter does not reach (see stedc: open debts).
| idea | why rejected | evidence |
|---|---|---|
| restore ORMQR = 128 | 1.50× slower at n = 512; 2.1× and 3.1× slower at n = 64, 128 | refutation |
| lift the CTA eigensolver above n = 32 (B3) | 85–211× slower than blocked | WP2 |
| global sb2st tile=2/subs=4 | costs float 1.15× | WP4 |
| wide two-stage kd (96–128) | kd = 32 still optimal at n ≥ 256 | kd sweep |
| freeze SB2ST_BACK_*; adopt stedc's own tpr/wgm | heuristic within 0.6%; stedc's winners cost syev 2.7% at n = 256 | retune |
| native Upper reductions | recover at most the mirror's 0.3–2.0% | Upper |
| the float small-n rule applied to complex | Jacobi 4–6× off the pace at n ≥ 20 | complex small n |
| a batch floor on the values-only rule | sent n = 1024 at batch 254 to the vendor, 2.75× | values routing |
| split latrd_grid_min_n per mode | same crossover in both modes | latrd gate |
| block Jacobi as a speed path (C1/WP8) | 2–11× slower; accuracy-only opt-in. At 1.65–2.1 TFLOP/s (about 3.5% of ~47 TFLOP/s SGEMM) a ~10× flop premium needs 17–21 TFLOP/s to break even | agent guide |
| C5 spectral divide and conquer; C6 real embedding of the Hermitian problem | cluster parallelism does not apply at saturating batch; 2× flops and memory | plan, not measured |
| data | where |
|---|---|
| small-n bake-off (phase A), both precisions (973 measurements) | benchmarks/results/syev_retune_phaseA.jsonl (Git LFS) |
| syev-context block-size sweep (phase B), float | benchmarks/results/syev_retune_phaseB_float.jsonl (Git LFS) |
| driver for both, with its resource caps | a one-off driver, removed |
| builds | 7911847 (bake-off), 12963a8 (grids), 0bb92fb (stedc regression), preset build/presets/cuda |
| commits | 80737301 plan, afc8ae35 baselines, c42e5ce9, 8c24626f WP4, 405b4120 WP3, f28a3fcd WP1, af8768b0 / be327f5b results (PR 68); 73eb5510 per-type routing (PR 65); 924b3a59 constant retune |
Phase B is complete except n = 256; only the float file is committed (double phase-B figures appear only on this page). Run git lfs pull to fetch the data (see the raw-data rules).
Reproduce one shape (positional n batch nb fuse jobz uplo; nb = 0 means the shipped default):
Kernel attribution that sees cuBLAS: nsys profile -t cuda -s none, then nsys stats --report cuda_gpu_kern_sum. BATCHLAS_KERNEL_TRACE does not see cuBLAS. Counters on the panel kernel:
1.18.0
