Status: current.
steqr, steqr_wg and steqr_cta solve batches of symmetric tridiagonal eigenproblems by implicit shifted QR/QL with deflation. This page gives the mathematics: deflation, scaling, shifts, the Givens bulge chase, split bookkeeping, the two CTA update schemes and eigenvector accumulation. It also covers the leaf solve of stedc and the fused small-n syev. CTA performance is in STEQR: the CTA tridiagonal solver; the stedc merge is in STEDC: the divide-and-conquer merge.
The solver takes a batch of tridiagonal matrices
\[T = \operatorname{tridiag}(e_1, \dots, e_{n-1}; d_1, \dots, d_n) \]
and computes
\[T = Q \Lambda Q^T, \]
where \(\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n)\) and, if requested, \(Q\) is the orthogonal eigenvector matrix.
Which path runs (steqr_dispatch, src/ops/syev/steqr/steqr.cc): steqr_cta when \(0 < n \le\) the device's largest sub-group size, otherwise steqr_wg.
Note: On ROCm the CTA path is disabled (should_use_cta). Its chunked sub-group operations gave wrong eigenvalues on AMD gfx1200.
Each matrix is stored as its diagonal \(d\) and off-diagonal \(e\) only:
\[T = \begin{bmatrix} d_1 & e_1 & & & 0 \\ e_1 & d_2 & e_2 & & \\ & e_2 & d_3 & \ddots & \\ & & \ddots & \ddots & e_{n-1} \\ 0 & & & e_{n-1} & d_n \end{bmatrix}. \]
No dense similarity transform is formed. Givens rotations chase a bulge through the band.
An off-diagonal entry deflates when
\[|e_i|^2 \le \varepsilon^2 |d_i| |d_{i+1}| + \mathrm{safmin}, \]
where \(\varepsilon\) is the machine precision and \(\mathrm{safmin}\) the safe minimum (LAPACK style). The code then sets \(e_i \leftarrow 0\), and \(T\) becomes block diagonal:
\[T = \begin{bmatrix} T_{11} & 0 \\ 0 & T_{22} \end{bmatrix}. \]
All paths use the relative criterion; there is no absolute deflation threshold.
Before sweeping an active block, the solver computes its norm
\[\|T\|_\infty^{(\mathrm{block})} \approx \max\left( |d_i|, |e_i| \right) \]
and rescales against overflow and underflow: by \(\alpha = \mathrm{ssfmax} / \|T\|\) if the norm is too large, and by \(\alpha = \mathrm{ssfmin} / \|T\|\) if it is too small but nonzero. After convergence the block is scaled back by \(\alpha^{-1}\).
Each implicit step uses a shift \(\mu\) from a \(2 \times 2\) principal block at the trailing (QR) or leading (QL) end.
For \(B = \begin{bmatrix} a & b \\ b & c \end{bmatrix}\),
\[\lambda_{\pm} = \frac{a+c}{2} \pm \sqrt{\left(\frac{a-c}{2}\right)^2 + b^2}, \qquad \mu = \begin{cases} \lambda_+ & \text{if } |\lambda_+ - c| < |\lambda_- - c|, \\ \lambda_- & \text{otherwise.} \end{cases} \]
The QL case uses the leading block and the eigenvalue closest to the top endpoint.
For a local \(2 \times 2\) block with leading entries \(d_1, d_2, e_1\), the CTA path can use
\[g = \frac{d_2 - d_1}{2 e_1}, \qquad r = \sqrt{g^2 + 1}, \qquad \mu = d_1 - \frac{e_1}{g + \operatorname{copysign}(r, g)}, \]
which avoids the cancellation of the quadratic formula.
\[G(c,s) = \begin{bmatrix} c & -s \\ s & c \end{bmatrix}, \qquad c^2 + s^2 = 1, \qquad G(c,s)^T \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}, \]
with \((c, s, r)\) from lartg(x, y). An implicit step applies a sequence of these as similarities \(T \leftarrow G_k^T T G_k\). Each similarity preserves symmetry, introduces a single bulge and chases it.
A similarity with an adjacent rotation on rows and columns \(i, i+1\) is
\[\widehat G_i = \operatorname{diag}(I_{i-1}, G(c,s), I_{n-i-1}), \qquad T^+ = \widehat G_i^T T \widehat G_i , \]
which changes only the band around indices \(i-1, \dots, i+2\). The local window before the update is
\[T_{\mathrm{loc}} = \begin{bmatrix} \alpha & \beta & 0 & 0 \\ \beta & d_i & e_i & 0 \\ 0 & e_i & d_{i+1} & \gamma \\ 0 & 0 & \gamma & \delta \end{bmatrix}, \qquad \alpha = d_{i-1},\ \beta = e_{i-1},\ \gamma = e_{i+1},\ \delta = d_{i+2}, \]
with entries that do not exist omitted. After the similarity,
\[T_{\mathrm{loc}}^+ = \begin{bmatrix} \alpha & c\beta - s b_{\mathrm{prev}} & -s\beta & 0 \\ c\beta - s b_{\mathrm{prev}} & d_i^+ & e_i^+ & -s\gamma \\ -s\beta & e_i^+ & d_{i+1}^+ & c\gamma \\ 0 & -s\gamma & c\gamma & \delta \end{bmatrix}, \]
with \(b_{\mathrm{prev}} = 0\) for the first rotation of a sweep. The entry \(-s\gamma\) is the new bulge. One step touches only \((e_{i-1}, d_i, e_i, d_{i+1}, e_{i+1})\).
The work-group path and the CTA EXP scheme use
\[\begin{aligned} d_i^+ &= c(c d_i - e_i s) - s(e_i c - s d_{i+1}), \\ d_{i+1}^+ &= c(c d_{i+1} + e_i s) + s(e_i c + s d_i), \\ e_i^+ &= c(c e_i + s d_i) - s(c d_{i+1} + s e_i), \\ e_{i-1}^+ &= c e_{i-1} - s b_{\mathrm{prev}}, \\ e_{i+1}^+ &= c e_{i+1}, \end{aligned} \]
and the outgoing bulge is \(b_{\mathrm{new}} = -s e_{i+1}\).
The first rotation is generated from the shifted corner \((d_{\mathrm{edge}} - \mu,\ e_{\mathrm{edge}})^T\), so the shift creates the bulge. Each later rotation is generated from \((e_{i-1}^{(\mathrm{current})},\ b_{\mathrm{prev}})^T\). The bulge moves one position and is not recreated.
For a block with endpoints \(d_\ell\) and \(d_r\), the solver converges from the endpoint of smaller magnitude: QR if \(|d_\ell| \le |d_r|\), QL otherwise. QR chases the bulge top to bottom and converges the bottom eigenvalue first. QL is the reverse. This is the LAPACK policy.
A deflated block of size 2, \(T_2 = \begin{bmatrix} a & b \\ b & c \end{bmatrix}\), is solved analytically. With eigenvectors requested, it also yields a rotation
\[Q_2 = \begin{bmatrix} c_2 & -s_2 \\ s_2 & c_2 \end{bmatrix}, \qquad Q_2^T T_2 Q_2 = \operatorname{diag}(\lambda_1, \lambda_2), \]
which is accumulated into the eigenvector matrix.
The solver iterates on maximal active blocks, never on the whole matrix. Let \(\mathcal Z = \{ i : e_i = 0 \}\). The active blocks are the maximal intervals \(I_k = [s_k, t_k]\) with \(e_j \neq 0\) for \(j = s_k, \dots, t_k-1\), so that
\[T = \operatorname{diag}(T_{I_1}, T_{I_2}, \dots, T_{I_p}), \]
and each \(T_{I_k}\) is processed independently.
On each outer pass, steqr_wg scans \(e\) and records each maximal run of nonzeros \(e_s, \dots, e_{t-1}\) as the block \([s, t]\) ((start_ix, end_ix), end_ix exclusive). An inclusive scan compacts the descriptors into one list of (batch index, start index, end index) triples, and all blocks across the batch run in parallel. Each outer pass has two stages:
A sweep can create new zeros, so the list is rebuilt on every pass.
steqr_cta builds no global list. Each sub-group partition holds one matrix in registers and walks its blocks with a running left boundary \(b_{\mathrm{next}}\). Each iteration sets \(b_{\mathrm{begin}} \leftarrow b_{\mathrm{next}}\), enforces \(e_{b_{\mathrm{begin}}-1} \leftarrow 0\) when that index exists, and takes the first zero to the right as the right boundary:
\[b_{\mathrm{end}} = \min\{ i \ge b_{\mathrm{begin}} : e_i = 0 \} \quad (b_{\mathrm{end}} = n \text{ if none}), \qquad b_{\mathrm{next}} \leftarrow b_{\mathrm{end}} + 1. \]
After each implicit step the solver re-tests deflation and locates the first (QL) or last (QR) new zero.
The nested loops (split, \(\ell\), sweep) are what the lockstep flat solver restructures.
QR uses the physical ordering. The first rotation annihilates the first subdiagonal of the shifted matrix,
\[G_1^T \begin{bmatrix} d_1 - \mu \\ e_1 \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}, \]
and each later rotation annihilates the bulge left by the previous one.
QL runs the same formulas on reversed local indices. Per active block, define
\[\tilde d_k = d_{n-k+1}, \qquad \tilde e_k = e_{n-k}, \]
and apply the QR-style chase to \((\tilde d_k)\) and \((\tilde e_k)\). This is a QL iteration in the original coordinates.
steqr_cta writes one implicit step in one of two equivalent forms, chosen by SteqrParams::cta_update_scheme. The default is EXP. (The enum comment that calls PG the "current default" is stale.)
PG (Parlett-Gray) propagates the scalars \((g, p, c, s)\) across the sweep instead of writing each local similarity. In QR form, step \(i\) is
\[\begin{aligned} f_i &= s_{i-1} e_i, \\ (c_i, s_i, r_i) &= \operatorname{lartg}(g_{i-1}, f_i), \\ g_i^{\star} &= d_i - p_{i-1}, \\ r_i^{\star} &= (d_{i+1} - g_i^{\star}) s_i + 2 c_i c_{i-1} e_i, \\ p_i &= s_i r_i^{\star}, \\ d_i &\leftarrow g_i^{\star} + p_i, \\ g_i &= c_i r_i^{\star} - c_{i-1} e_i . \end{aligned} \]
QL is the mirrored recurrence with the sweep reversed.
The CTA solver runs one tridiagonal per chunk of \(P\) lanes, several chunks per warp. In the nested loop structure, chunks realign only at loop exits, so the warp pays the slowest chunk's sweep count per eigenvalue. steqr_cta_solve_flat moves the sweep out of the nest. Each chunk first settles (opens a block, tests, advances past deflated eigenvalues, solves 2x2 blocks) until it has a sweep \([\ell, m]\) or is done. Then all chunks with a sweep run it in the same pass. The per-chunk operation sequence is unchanged, so results are bitwise equal to the nested solver. The flat solver ships only for some (type, \(P\)) pairs. The kSteqrCtaFlatPays gate, the sub-group-wide form and the measured losses are in Lockstep flat solver.
With back_transform = false the solver starts from \(Q \leftarrow I\). Otherwise it treats the input matrix as a basis and right-applies the rotations. Each local rotation on columns \(k, k+1\) does
\[Q_{:,\{k,k+1\}} \leftarrow Q_{:,\{k,k+1\}} G(c,s). \]
In QL the physical column order is reversed relative to the sweep, but each operation is still a right multiplication by an adjacent rotation. steqr_wg stores a sweep's rotations and applies them afterwards. steqr_cta usually applies them at once to a shared-memory cache of \(Q\).
After convergence \(\lambda_i \leftarrow d_i\). With sort = true the eigenpairs are sorted ascending, \(\lambda_1 \le \lambda_2 \le \dots \le \lambda_n\), and the same permutation is applied to the columns of \(Q\).
In short,
\[T_0 \xrightarrow{\text{deflation + implicit shifted QR/QL steps}} T_k = \Lambda, \qquad Q = G_1 G_2 \cdots G_p, \qquad Q^T T_0 Q = \Lambda . \]
1.18.0
