Status: current · written 2026-08-01, corrected in place where measurement disagreed.
syevx returns \(k\) eigenpairs of a symmetric or Hermitian matrix, by extreme end or by an index or value range. It dispatches to one of four algorithms. This page gives the algorithms, the routing rule, the cost model, the verdict on each candidate, and the rules the dispatcher enforces. Measured crossovers are on SYEVX. Range semantics are on syevx range selection: design record.
| Algorithm | Input | Method | Chosen by Auto when |
|---|---|---|---|
| Direct | dense; all types, ranges and jobz | syev, then select | \(n \le 64\), or eigenvalues only |
| DirectSubset | dense | band reduction, stebz, stein, back-transforms on \(k\) columns | eigenvectors at large \(n\) and batch |
| LOBPCG | CSR (dense on explicit request) | batched LOBPCG, optional ILU(k) or Jacobi preconditioner | sparse input, extremal ranges only |
| Filtered | CSR or dense | scaled Chebyshev filter with Rayleigh-Ritz | never; explicit request only |
On dense input with eigenvectors, Auto picks DirectSubset when \(n \ge 1024\) and \(n \cdot \mathrm{batch} \ge 128 \cdot 1024\), and Direct otherwise. The gates are kSyevxSubsetMinN and kSyevxSubsetMinWork in src/ops/syevx/syevx.cc. The rule has no \(k/n\) term: the DirectSubset advantage is flat in \(k\) from 0.8 % to 25 % of the spectrum.
Routing gave most of the practical gain (3.0x to 3.7x geometric mean). DirectSubset adds up to 2.4x, at large \(n\) and batch only.
Leading-order flops for real symmetric \(n \times n\) with \(k\) wanted eigenpairs. These are estimates; measurement did not confirm their performance predictions.
| Path | Flops |
|---|---|
| Full, one-stage, with vectors | \(\approx 4.7n^3\) (sytrd \(\tfrac43 n^3\), stedc \(\tfrac43 n^3\), ormtr \(2n^3\)) |
| Full, two-stage, with vectors | \(\approx 6.7n^3\); the back-transform is applied twice |
| DirectSubset | \(\approx \tfrac43 n^3 + 4n^2k\) |
| Chebyshev subspace iteration, dense, degree \(d\) | \(4dn^2k\) per outer pass |
| LOBPCG, dense, \(it\) iterations | \(2\,it\,n^2k + O(it\,nk^2)\) |
For \(n \le 64\) a CTA-resident solver is bound by occupancy, shared-memory traffic and sweep count, not flops. A subset kernel would land within noise of syev_cta, so syevx calls syev and selects for \(n \le\) kSyevxSmallN (64). The path exists for API compatibility, not speed.
Not built:
This is where DirectSubset applies: \(n\) from about 128 to 4096, batch 8 to 1000.
| Stage | Component |
|---|---|
| dense to band | sytrd_sy2sb |
| band to tridiagonal (keeps Q2) | sytrd_sb2st_hh |
| subset eigenvalues | stebz |
| subset eigenvectors | stein |
| back-transform through Q2 | unmqr_hb2st, on \(k\) columns |
| back-transform through Q1 | ormqr_blocked, on \(k\) columns |
No vendor routine is a batched subset eigensolver. That gap is the product case for DirectSubset.
Bisection counts eigenvalues below a shift from Sturm sign changes. Each eigenvalue is independent, so it maps to one lane with no communication. Inverse iteration's weakness, loss of orthogonality in tight clusters, costs \(O(nk^2)\) to repair. With \(k \ll n\) that does not bind. Built as stebz and stein.
\(O(n^2)\) with no reorthogonalization, but hard to implement correctly. For medium \(n\) the saving is negligible against the \(O(n^3)\) tridiagonalization. LAPACK's dsyevr falls back to dstebz and dstein on some subset paths.
All work is GEMM (dense) or SpMM (sparse). It needs a spectral bound estimate, reasonable separation, and a degree uniform across the batch. It is the strong path on sparse input, not dense.
Not built. The polynomial filter serves as the expansion operator, so no preconditioner is needed. Filtered covers much of the same ground. Unpreconditioned LOBPCG stagnates, and that is the default configuration unless the caller supplies ILU(k) or Jacobi.
Kept for sparse \(A\) with a good preconditioner, where no \(O(n^3)\) alternative exists. The hardening follows Duersch, Shao, Yang and Gu (SISC 2018): Hetmaniuk-Lehoucq basis selection, svqb/svqbDrop with condition-based dropping, and a backward-stable criterion. The old criterion \(\|r\| / (\|x\|\,|\lambda|)\) was not backward stable. Open: staircase locking, warm start, and a Chebyshev preconditioner.
Slicing adds orthogonality bookkeeping that the batch dimension already makes unnecessary. FEAST and shift-and-invert need a solve with \(A - zI\) per shift per matrix, strictly more expensive than the direct route for batched dense input. The syevr and generalized sygvx non-goals are on syevx range selection: design record.
validate_syevx_preconditioner_params and validate_syevx_range_params run in syevx and syevx_buffer_size before syevx_select_algorithm. They check the problem, not the algorithm. Otherwise an illegal combination on dense input would reach Direct or DirectSubset, which ignore preconditioners.
value_range_reportable is false only for solve overloads without an m output. A short m span is rejected on the host, because the device write it would cause is unchecked in release builds.
algorithm_from_env resolves BATCHLAS_SYEVX_ALGORITHM against SyevxParams::method. An environment default degrades where an explicit request throws, so that a whole suite can run on one algorithm for diagnosis.
An algorithm can be substituted, because only performance changes. A preconditioner cannot: an ILU(k) factor the caller built has no substitute, and silently ignoring it is a correctness surprise. Only LOBPCG uses the preconditioner.
Auto never picks a Jacobi form, since neither is a free win on the measured matrices. JacobiShifted is offered because it is safe in general and is the only Jacobi form legal with find_largest.
Direct and DirectSubset need dense \(A\), so on CSR they degrade to LOBPCG. Filtered works on CSR through spmm, so an explicit Filtered on CSR is honoured. Auto still lands on LOBPCG.
The CSR branch must stay below the non-extremal range check. Otherwise sparse input with an interior range would return extremal eigenpairs instead of throwing. LOBPCG and Filtered only see extremal ranges, so syevx fills m with the resolved block size before the solve, on the same in-order queue.
syevx_direct and syevx_direct_subset are public, so they re-check Index and Value ranges. An out-of-range block would index past the \(n\)-entry eigenvalue array. The check throws std::invalid_argument; syevx_resolve_range clamps as a second defence. DirectSubset validates before the reduction, since stebz would otherwise reject the range after \(O(n^3)\) work.
neigs above \(n\) is accepted and clamped. neigs = 0 is rejected by DirectSubset, because stein needs \(k \ge 1\).
info is 0 for every item. That is a caveat, not a check: stebz gives no exit status, and stein runs a fixed number of inverse iterations with no convergence test. A real flag needs a residual check, which is not implemented. Writing 0 is deliberate, so that an untouched span cannot be read as "converged".
One outer iteration:
\(p_m\) maps the unwanted interval into \([-1, 1]\), where \(|T_m| \le 1\). The wanted end lies outside, where \(T_m\) grows like \(\cosh(m\,\mathrm{acosh}\,x)\). The cut is the first unwanted Ritz value. The filter needs only matvecs, so no preconditioner or factorization.
The recurrence is the scaled Zhou/Saad form:
\[Y_1 = \frac{\sigma_1}{e}(A - cI)X, \qquad Y_{j+1} = \frac{2\sigma'}{e}(A - cI)Y_j - \sigma\sigma' Y_{j-1}. \]
Unscaled \(T_m\) overflows: \(T_{25}(3) \approx 10^{38}\), which is float infinity. The \(\sigma\) factors keep intermediates near unit magnitude. The block is orthonormalized afterwards.
Guard directions (extra_directions, default max(2, k/4)) keep the last wanted pair from converging slowly. The start block is a hash of the index, so it is reproducible across runs and backends.
Normalise at the extreme Ritz value, not the Gershgorin bound. Gershgorin overestimates a random symmetric matrix's spectral radius ( \(O(n)\) against \(O(\sqrt n)\)). Normalising there underflows the block at degree 40, and orthogonalizing a zero block gives NaN.
Cap the degree by precision. The filter amplifies the most-wanted direction over the least-wanted one by
\[\frac{\cosh(d\,\mathrm{acosh}\,y_{far})}{\cosh(d\,\mathrm{acosh}\,y_{edge})}. \]
Past what the working precision holds, the block is rank-deficient and the Cholesky-based ortho returns NaN. The cap is
\[d_{cap} = \frac{-\tfrac14 \log \epsilon}{\mathrm{acosh}\,y_{far} - \mathrm{acosh}\,y_{edge}}, \]
a ratio budget of \(\epsilon^{-1/4}\) (about 300 in float). A budget of \(\epsilon^{-1/2}\) still failed at degree 40. When \(y_{edge} \le 1\) the edge term is taken as 0, the worst case, so the cap still applies. Each item's cap is reduced by min across the batch, so every GEMM shape stays uniform.
Behind BATCHLAS_SYEVX_FILTER_DEGREE_AUTO=1, off by default. The degree that finishes in one more outer iteration is
\[d_{need} = \frac{\log(\text{worst residual} / \text{target})}{r}, \qquad r = \mathrm{acosh}\,y_{edge}, \]
in the units of the convergence test. Each item's \(d_{need}\) is clamped to [kAutoDegreeMin, kAutoDegreeMax] = [10, 40], reduced by max across the batch, then trimmed by the batch-minimum cap. The lower bound is the old constant, so enabling the derivation never makes a step less selective.
Precedence: BATCHLAS_SYEVX_FILTER_DEGREE > params.filter_degree > derived > kDefaultFilterDegree. An explicit value turns the derivation off. BATCHLAS_SYEVX_FILTER_DEGREE_AUTO=0 restores the constant exactly.
Trap: BATCHLAS_SYEVX_FILTER_DEGREE_AUTO is parsed with atoi, not env_truthy, so true and on read as false. This is kept, so existing spellings keep their meaning.
fp64 CPU failures on the primary machine are usually the broken OpenBLAS Cooperlake dgemm. Check a known-good BLAS first.
| Component | State |
|---|---|
| Direct, DirectSubset, stebz, stein (Index and Value ranges) | done |
| Filtered | done; opt-in only |
| LOBPCG hardening | mostly done; open: staircase locking, warm start, Chebyshev preconditioner |
| Fused filtered CTA solver for small \(n\) | not built; gated on a gap estimate |
Known risks: inverse iteration on variable-size clusters (orthogonality is a hard test gate); a bad spectral bound in Filtered gives no acceleration or overflow (bounds are conservative, the degree is clamped, and divergence is detected; two NaN modes were found and fixed); and per-feature choices across the batch (max-over-batch or bucketing, not single-matrix adaptivity).
Bisection, inverse iteration, MRRR:
Polynomial filtering:
LOBPCG:
Dense eigensolver performance on GPUs:
1.18.0
