Status: current · RTX 4090 (sm_89), CUDA backend, float unless a row says otherwise
Measurements behind the partial symmetric eigensolver syevx: the dense crossovers that set Auto routing, the sparse Filtered-vs-LOBPCG sweep, the Chebyshev-degree A/B, and the LOBPCG hardening experiments. The design is in syevx: design of the partial symmetric eigensolver; range selection is in syevx range selection: design record.
Raw data:
LOBPCG iteration counts in the hardening sections are CPU (NETLIB) builds with BATCHLAS_SYEVX_CHECK_EVERY=1, unless a row says otherwise.
The dense Auto rule, from syevx_select_algorithm (src/ops/syevx/syevx.cc:302-320):
| condition | route |
|---|---|
| n <= kSyevxSmallN (64) | Direct |
| eigenvalues only (jobz != EigenVectors) | Direct |
| eigenvectors, real T, n >= kSyevxSubsetMinN (1024) and n * batch >= kSyevxSubsetMinWork (128*1024) | DirectSubset |
| everything else dense | Direct |
| sparse (MatrixFormat::CSR) | LOBPCG, or Filtered only when explicitly requested |
k does not enter the rule. Auto never picks Filtered, even where it wins on sparse input (Filtered vs LOBPCG on CSR).
Before commit fc4fa7a9, syev's Auto order put the in-house blocked solver ahead of the vendor, so every dense Direct solve ran it. The fix, syev_prefer_vendor (syev.hh:181 when written; no longer in include/ or src/, routing is now flat selection, see syev), sent syev to the vendor on CUDA. Over n = 64..1024, batch = 1..512, the vendor won 24 of 25 shapes (1.27x to 15.4x). Blocked leads by up to 1.37x on 320 <= n <= 640 with batch >= 128; that carve-out is kept as measured, since the boundary is not monotone in n.
End to end on syevx Auto (neigs = 8, n = 64..2048, batch = 1..128): geometric mean 3.72x (eigenvalues only) and 3.00x (with eigenvectors); best 33.7x, worst 0.99x.
BM_SYEVX_CrossoverVectors (benchmarks/syevx_benchmark.cc), RTX 4090, float, eigenvectors. Cells are Direct time / DirectSubset time, so above 1 means DirectSubset wins.
| shape | b=1 | b=4 | b=16 | b=64 | b=256 |
|---|---|---|---|---|---|
| n=1024, k=8 | 0.09 | 0.34 | 0.36 | 1.00 | 2.40 |
| n=2048, k=8 | 0.06 | 0.28 | 0.43 | 1.12 | 1.98 |
Spectrum position does not enter either path's cost (syevx range selection: design record). BM_SYEVX_RangePosition that would check this has not been run (open debts).
Superseded. Recorded 2026-08-02 (commit c25f5050). Its Direct column is the blocked syev, not cuSOLVER, so ratios against Direct from this grid are invalid. CSV columns are n, batch, neigs, algorithm with algorithm 1..4 = Direct, DirectSubset, Filtered, LOBPCG. The DirectSubset and Filtered columns remain valid:
BM_SYEVX_SparseShape and BM_SYEVX_SparseGap (benchmarks/syevx_benchmark.cc), 2026-08-06, commit 826bc5cb. CUDA, float, eigenvectors, 16 nonzeros per row, diagonal_boost 10, k about 1 % of n, tolerance 1e-6, iterations = 200, Chol2 ortho, guard block max(1, neigs/4). Ratio is LOBPCG time / Filtered time, so above 1 means Filtered wins.
| n | batch | largest end | smallest end |
|---|---|---|---|
| 1024 | 1 | 8.8 | 2.7 |
| 1024 | 64 | 2.2 | 1.5 |
| 2048 | 1 | 10.9 | 4.6 |
| 2048 | 64 | 1.4 | 0.73 |
| 4096 | 64 | 0.88 | 0.51 |
Full 9-shape grid: benchmarks/results/syevx_sparse_shape_rtx4090.csv.
Auto stays on LOBPCG for sparse input: the batch-64 crossover is real, and one operator family (random sparse Hermitian with an imposed dominant diagonal) is too thin a basis for a routing rule. The next step is a mesh/stencil family.
Note: Explicit params.method = Filtered on CSR was silently ignored before this sweep, so runs labelled Filtered on CSR earlier actually ran LOBPCG. Guarded by SyevxRangeRoutingTest.SparseHonoursExplicitFilteredRequest.
BM_SYEVX_SparseGap: n = 2048, batch 8, k = 20, diagonal_boost 1 to 30 (raw: benchmarks/results/syevx_sparse_gap_rtx4090.csv).
syevx_filtered uses kDefaultFilterDegree = 10 (src/ops/syevx/kernels/syevx_filtered.cc:56) unless the degree is pinned (BATCHLAS_SYEVX_FILTER_DEGREE or params.filter_degree). The per-matrix derived degree behind BATCHLAS_SYEVX_FILTER_DEGREE_AUTO=1 must not ship.
RTX 4090, float, n = 1024, neigs = 8, derived / fixed-10 (above 1 means derived wins). EigenVectors agrees with NoEigenVectors to within 0.02x.
| batch | 1 | 2 | 4 | 8 | 16 |
|---|---|---|---|---|---|
| ratio | 2.18x | 2.25x | 0.70x | 0.50x | 0.48x |
The loss tracks batch size: the source takes the maximum derived need across the batch, clamps to [10, 40], then trims by the minimum precision cap, so every matrix runs at the most demanding item's degree. Unverified; a per-iteration degree trace at batch 8 would settle it. The constant 10 has never been tuned on this hardware.
Every tier except Tier 5 exists (plan: syevx: design of the partial symmetric eigensolver).
syevx_direct_subset (src/ops/syevx/kernels/syevx_direct_subset.cc): sytrd_sy2sb → band-to-tridiagonal → stebz → stein → back-transforms on the k selected columns. Real types, dense input only.
syevx_filtered (src/ops/syevx/kernels/syevx_filtered.cc): spectral bounds → Chebyshev filter → ortho → Rayleigh-Ritz, looping until the wanted residuals converge. Dense and CSR, real and complex, no preconditioner, not limited to the smallest eigenpairs. Reachable only explicitly (params.method or BATCHLAS_SYEVX_ALGORITHM=filtered).
Two numerical rules, each caught by a test:
Tests assert eigenvalues against a reference syev, because a self-consistent pair does not prove it is one of the wanted pairs.
Known cost: the convergence test reads a device flag on the host once per outer iteration, the same defect as LOBPCG's host synchronization.
Derivations and degree precedence: the Filtered design notes.
syevx_filtered starts from Gershgorin bounds (lo/hi are the min/max over rows of \(a_{ii} \mp \sum_{j \ne i}|a_{ij}|\)). The bounds must be conservative; a loose interval costs filter sharpness, never correctness.
The shipped kernel gives each work-item one row (kBoundsWG = 64 rows per work-group), so a work-group reads a contiguous column run; a second kernel reduces the per-(matrix, row-group) partials. The legacy one-work-item-per-matrix kernel took 80.2 % of the solve at n = 1024 (RTX 4090).
Tier 4. Iteration counts are CPU (NETLIB) unless stated.
The loop used to read unified memory on the host every iteration (about 30 forced drains per solve). It now drains only at the convergence check (every check_every, default 4, BATCHLAS_SYEVX_CHECK_EVERY) or for instrumentation that cannot be staged on the device. Overshooting by up to check_every - 1 iterations costs less than the drains. Verified at {1, 4, 16}.
The residual kernel writes a staging buffer carved from the pool, and one host pass after the loop scatters into the caller's spans.
One pass forms R and accumulates the per-column norms, replacing 2 * neigs sequential group reductions (128 barriers at neigs = 64).
The criterion is \(\|r\| / (\|Ax\| + |\lambda|\,\|x\|)\) (Duersch et al.). The former \(\|r\| / (\|x\|\,|\lambda|)\) never converged for \(\lambda \to 0\).
Column masking (variant a) is implemented and off by default: BATCHLAS_SYEVX_SOFT_LOCK=1 enables it (the parser for this knob is inverted; see include/batchlas/settings.hh). The batch-wide staircase (variant b) is not implemented.
The residual kernel flags col_converged[bid][j] when residual_j <= 0.1 * tol (BATCHLAS_SYEVX_LOCK_FACTOR). ortho mixes columns (R becomes \(R\,U^{-1}\), U upper triangular), so masking on one side only fails: masking after ortho deletes content the unconverged columns were orthogonalised against (NaN by iteration 26 at n = 128, k = 8, tol = 1e-4, where the unlocked run converged at 14), and masking before ortho makes \(R^T R\) singular, so potrf returns NaN.
The implementation masks both sides. Before ortho, converged columns get a deterministic pseudo-random vector scaled to the largest surviving residual; after ortho they are zeroed. A zeroed column left in StAS adds a spurious eigenvalue 0, which is fatal for find_largest = false on positive-definite A, so each locked diagonal gets a sentinel \(\pm(\|StAS\|_F + 1)\) that sorts it to the unwanted end.
The locking threshold must be stricter than the tolerance. A column locked at exactly tol drifts back above it and oscillates; that made the final residual worse in 3 of 8 ILU(k) sweep cases (tests/iluk_tests.cc).
Measured, no win: CPU iteration counts moved by at most +1 across six cases (n = 256, b = 4, k = 16, tol 1e-5: 19 to 20). Masking saves no flops, since masked columns still occupy the fixed shapes. The regime where it should help is unreachable on CPU, because Backend::NETLIB forces Householder. A GPU A/B under Chol2 is still owed before the default could flip.
Variant (b) would need the active blocks to be non-contiguous, which breaks the single S({0,n},{0,Nvecs}) slice every gemm uses. Nvecs is computed in three places that must agree, so collapse those into one value first.
extra_directions == 0 means max(2, neigs/4), the convention filter_degree also uses; an explicit width wins. Iterations to converge (find_largest, tol 1e-5, batch 4):
| n / neigs | no guard | guard | change |
|---|---|---|---|
| 64 / 4 | 32 | 20 | -38 % |
| 256 / 4 | 78 | 39 | -50 % |
| 256 / 8 | 50 | 33 | -34 % |
Across eight shapes (n = 64 to 256, neigs = 4 to 16) the guard cut iterations by 34 % to 50 %, at a cost linear in the extra width. Comparisons of LOBPCG with other paths use a guard block too.
iluk_apply reads R's strided slice directly, so the R_contiguous repack is gone. R_preconditioned stays separate: the forward solve writes out while still reading rhs.
The two per-iteration batch-wide launches that reversed columns for find_largest are gone. X stays in syev's ascending order; no replacement was available (MatrixView clamps ld, so there is no negative stride for vendor gemm; syev returns ascending only; gemm has no permutation). Column j pairs with lambda[eig_offset + j], where eig_offset = num_eigvals - block_vectors when find_largest. The largest-first order comes from the X_best snapshot. One cold-path launch remains, for params.iterations == 0.
SyevxLobpcgVectorsTest checks per-pair residuals against the original A. Tests that request only NoEigenVectors cannot catch a column-to-eigenvalue mismatch.
| n, k | p=0 | p=1 | p=4 | p=8 |
|---|---|---|---|---|
| 64, 4 | 22.67 | 20.67 | 18.67 | 16.00 |
| 64, 16 | 8.33 | 7.00 | 4.67 | 3.00 |
| 256, 8 | 36.00 | 35.33 | 31.00 | 27.33 |
p = 4 gives 12 % to 44 % fewer iterations (n = 64, k = 16: 23.4 to 11.3 ms). p = 8 is better still, but the default stays at 4: an over-compressed block is rank-deficient, and the Cholesky ortho then returns NaN.
Negative result: the same start with \(\sigma I - A\) for find_largest = false was flat (23.00 to 21.67 at n = 64, k = 4) and is not implemented.
SyevxParams::preconditioner_type (SyevxPreconditioner; Auto keeps the old behaviour). diag(A) is extracted once, and the apply is one in-place kernel over R. DISABLED_IterationSweep in tests/syevx_tests.cc reproduces the fp32 counts.
The projected problem is a batched 3k x 3k dense eigensolve per iteration. syev's Auto gave every n <= 32 to CTA, so it never reached the vendor. At n = 30, batch 8, float, eigenvectors: CTA 229.6 µs per call against cuSOLVER's 103.7 µs (2.21x), and nsys attributes 29.4 % of LOBPCG GPU time to it.
The retired (2026-10-05) syev_prefer_vendor_over_cta knob (BATCHLAS_SYEV_CTA_MAX_N, default 32 via syev_cta_max_n_default_for<T>(); none of these exist in include/ or src/ now) would have sent eigenvector-mode CUDA solves with n above the threshold to the vendor. It shipped off: it measured 1.10x to 1.15x through LOBPCG at batch 8 but flipped a marginal case in ILUKTests.SyevxInstrumentationAndPreconditioner. The open owner decision is on the syev page. BATCHLAS_SYEVX_PROJECTED_VENDOR forces the vendor for the projected solve; it changes workspace size, so it must not change between a buffer-size query and the call.
Since 2026-10-09 (deep sm_120 tables). The tuned syev tables rank families by time only, and the deep sm_120 run ranks cuSOLVER (vendor) or jacobi first for many small eigenvector-mode cells. LOBPCG cannot afford that trade. On threadripper02 (RTX PRO 6000, sm_120), float, n = 16, cuSOLVER's eigenvalues are about 1.9e-5 off, against 2.7e-6 for the CTA solver. With the projected solves on vendor or jacobi, ILUKTests.SyevxInstrumentationAndPreconditioner (iluk_k3, case d0.06_b0.5) stalls at a residual of 5.2e-6 against its 4.2e-6 target. It passes on CTA.
projected_syev_pin (src/ops/syevx/kernels/syevx_lobpcg.cc) therefore pins every projected solve (XtAX, StAS and its restart size) to cta_fused when the CTA solvers can take the shape (n <= 32, a GPU with 32-wide sub-groups). The buffer-size query uses the same pin, so the workspace and the run branch identically. cta has the same can_run as cta_fused, so it is never a fallback rung. Above n = 32 the table decides. An explicit route (BATCHLAS_SYEV_ROUTE or an outer select::ScopedPin) still wins, which keeps the five Pinned/SyevxProjectedSyevRouteTest.* cases exercising each family. SyevxProjectedSyevPin.CtaFusedAtOrBelow32EvenWhenTheTableRanksTheVendorFirst guards the pin through the selection trace, under a synthetic table that ranks the vendor first.
1.18.0
