Status: current.
Batched SVD uses one-sided Jacobi (gesvdj_cta) up to n = 64 (complex<double> with vectors up to 32). Above that it bidiagonalises with gebrd, solves the bidiagonal problem with bdsdc on the batched stedc, and back-transforms with ormbr. This page gives the design rules, the kernel's local-memory and scaling rules, and the traps. Measurements are in the GESVD performance page.
The accuracy target is cuSOLVER gesvdjBatched, which gives relative accuracy (Demmel–Veselić):
\[ \frac{|\delta\sigma_i|}{\sigma_i} \le O(\varepsilon)\,\kappa(A_c), \]
where \(A_c\) is the column-equilibrated matrix. Paths that form the tridiagonal of \(B^T B\) take \(\sigma = \sqrt{\lambda}\) and square the condition number (see the normal-equations defect).
Hestenes one-sided Jacobi applies 2x2 rotations to column pairs of A until the columns are orthogonal. Then \(\sigma_i = \|A_i\|\), \(U_i = A_i / \sigma_i\), and V is the accumulated rotation. Each pair's Gram is recomputed from the current columns, so the test is relative:
\[ |a_{pq}| > \mathrm{tol}\cdot\sqrt{|a_{pp}|\,|a_{qq}|}, \qquad \mathrm{tol} = t\cdot n\cdot\varepsilon \]
An absolute test \(|a_{pq}| \le \mathrm{tol}\cdot\max|a_{kl}|\) loses this accuracy (LAWN 169, Remark 2.2).
Why Jacobi here:
The kernel keeps the relative stopping test and skips converged pairs from LAPACK xGESVJ. It does not use column equilibration (scaling) or de Rijk ordering (evidence).
References: Hestenes (1958); Demmel & Veselić, SIAM J. Matrix Anal. Appl. 13(4) (1992); Drmač & Veselić, LAWN 169/170; Golub & Van Loan, Alg. 8.5.1.
bdsdc is not a port of LAPACK dlasd*. For B upper bidiagonal with diagonal d and superdiagonal e,
\[ T = \begin{pmatrix} 0 & B^T \\ B & 0 \end{pmatrix} \]
has eigenvalues \(\pm\sigma_i\). Under the perfect-shuffle permutation it is a symmetric tridiagonal matrix of order 2n with zero diagonal and off-diagonal \((d_0, e_0, d_1, e_1, \dots, d_{n-1})\), which the batched stedc solves. Nothing is squared. A native bidiagonal divide and conquer would be about 2000 lines of new device code. Solving the doubled problem with stedc is about 1000x faster than the sequential bdsqr (table).
Extraction is exact. For \(\sigma \ne 0\) the \(+\sigma\) eigenvector is \((v; u)/\sqrt{2}\) interleaved, so even and odd rows give v and u. The only residue is a possible sign flip on u when the \(\pm\sigma\) vectors mix by more than 45°. It costs at most \(2\sigma\) in the residual, and only when \(\sigma \approx \varepsilon\|B\|\).
Null space. Repeated zero singular values, or a \(\pm\sigma\) pair too close to resolve, make eigenvectors collapse onto the same (v, u). bdsdc_repair_degenerate rebuilds each column with \(\sigma_i < 500\,\varepsilon\,\sigma_{\max}\) by Gram–Schmidt against the final columns. Any orthonormal completion satisfies \(B = U S V^T\) to within \(2\sigma\). The threshold has no factor of n. The rejected criteria are in the perf page.
bdsdc writes the vectors it computes, so \(U^T U = I\) is a property of the method. tests/bdsdc_tests.cc checks it on graded spectra, a hostile gebrd-like generator, rank deficiency of multiplicity 1-5, and the all-zero matrix. bdsqr stays behind BATCHLAS_GESVD_BIDIAG=bdsqr: it keeps relative accuracy for tiny \(\sigma\) at 300-400x the cost.
Each round has n/2 disjoint pairs, and each needs three length-m dot products before its rotation is known. Three mappings were compared at m = n = 32, P = 32, float, with vectors. These counts are projected, not measured. syev_jacobi_cta is 266.
| design | per round | vs syev |
|---|---|---|
| lane = row, reduce-scatter Gram (chosen) | ~181 | 0.68 |
| lane = column | ~251-264 | 0.94-1.00 |
| Gram-resident ( \(G = A^H A\) kept by congruence) | 334 + refresh | 1.26-1.39 |
Lane = row wins because both operands of a rotation sit in one lane, so no partner shuffles are needed. Converged pairs are free: the skip if (sk == 0) continue; is warp-uniform, which lane = column would not be. Lane r alone reads and writes row r of A and V, so the update has no cross-lane hazard.
The Gram-resident design is rejected. Congruence updates give normwise error \(\sim\varepsilon\sigma_{\max}^2\) rather than columnwise error, so the relative threshold degrades to an absolute test. It also needs 51,200 B per complex<double> problem, above the 49,152 B the design assumed.
\(P\) is the partition width (lanes, at most 32, the sub-group width). \(C\) is the tile capacity. For every rung up to 32, \(P = C\). The 64 rung is \((P, C) = (32, 64)\), with kRPL = C/P rows per lane. The rung is keyed on \(\max(m,n)\): 4, 8, 16, 32, then (32, 64).
Bytes per problem with V resident, plus the work-group pair table:
| T | C = 32 | C = 64 (per problem + table) |
|---|---|---|
| float | 8,768 (+992) | 33,920 + 4,032 = 37,952 |
| double | 17,472 (+992) | 67,712 + 4,032 = 71,744 |
| complex<float> | 17,344 (+992) | 67,712 + 4,032 = 71,744 |
| complex<double> | 34,624 (+992) | 134,784 + 4,032 = 138,816 |
The device reports 101,376 B of local memory. complex<double> with vectors therefore does not launch at C=64 and is capped at 32. gesvd_jacobi_max_dim (src/ops/gesvd/gesvd_native.hh) sets the cap, and both the kernel guard gesvdj_cta_max_dim and can_run read it. Other types drop to 1-2 work-groups per SM, against 10 at C=32. Occupancy, not the cap, is the binding limit.
The Gram phase computes kGramChunk = P/2 dot products per chunk and scatters them so that pair k lands in lanes 2k and 2k+1. Scattering V values over L lanes takes \(\log_2 V\) scatter steps plus \(\log_2(L/V)\) all-reduce steps. At P=32, V=16 that is 4 + 1, not 5 + 0:
Written as five halving steps, half is 0 at the fifth step, the loop never runs, and each dot product sums only 16 of the 32 rows. max(1, half) reads a stale g[1]. The exact-norm reduction (32 values over 32 lanes) is 5 + 0 and correct.
The cost is 16 shuffles for all sixteen dot products, against 80 for independent butterflies. The tree sum's error is \(\sim\log_2(m)\,\varepsilon\), against \(\mathrm{tol} = n\varepsilon\). A lane-sequential sum would give \(\sim m\varepsilon\), which causes threshold churn.
At C=64 a round runs as two chunks of 16 pairs. The chunks touch disjoint columns, and the live register set is 80 T rather than 160. The exact-norm reduction stays 32 lanes wide, over C/P passes.
Per-column equilibration is wrong as an algorithm step. Scaling columns by D computes the SVD of \(AD\), and \(D^{-1}W\) is not orthogonal, so the result does not factor A. \(\kappa(A_c)\) is an analysis quantity only.
The kernel applies one global power-of-two rescale centred on the geometric mean of the column norms:
\[ \beta = 2^{-e}, \qquad e = \operatorname{round}\!\Big(\tfrac14\big(\log_2 \max_c\|A_c\|^2 + \log_2 \min_{c,\ \|A_c\| > 0}\|A_c\|^2\big)\Big). \]
A power of two keeps the round trip exact. For float, geometric centring tolerates a column-norm ratio of \(1.69\times10^{38}\), against \(9.2\times10^{18}\) for centring on the maximum. An all-zero input keeps \(\beta = 1\). The column-norm epilogue uses a scaled sum of squares. Without it, an underflowing norm gives \(\sigma = 0\) and fabricates a U column.
Note: the \(\max\)/ \(\min\) reductions read Nrm_local only for the first 32 columns. On the C=64 rung, columns 32..63 do not affect \(\beta\). Correctness holds, but the overflow headroom covers only those 32 columns. Observed, not measured.
The zero threshold is relative to \(\sigma_{\max}\) only: \(\mathrm{tol}_0 = \mathtt{zero\_sigma\_multiplier}\cdot\varepsilon\cdot\sigma_{\max}\). The absolute form \(\varepsilon\max(1, \sigma_{\max})\) declares every \(\sigma\) zero when \(\sigma_{\max} = 10^{-10}\).
Completion. Columns CC..RR-1 of a full U, and any column below \(\mathrm{tol}_0\), are filled in-kernel by two passes of classical Gram–Schmidt against the accepted columns. The trial cursor runs across output columns and accepts above \(1/(2 \cdot RR)\). A per-column restart leaves the last columns unfilled, because a canonical vector's residual is about \(d/RR\) when d dimensions remain. The branch is gated on a warp-uniform predicate.
Thin (SvdVectors::Thin) sets left_cols = CC and skips completion. A numerically deficient column is still repaired.
m < n solves \(A^H\), transposed at load time. Solving \(A^T\) instead gives \(\overline{V'} S U'^T\), whose conjugations are invisible in real arithmetic and wrong for complex.
SvdVectors (include/batchlas/blas/enums.hh) follows LAPACK jobu/jobvt: None, All ('A') and Thin ('S'). For tall-skinny input All is unusable. A 10000 x 32 problem needs a 10000 x 10000 U, 400 MB per matrix in float.
For \(k = \min(m, n)\), a thin U equals the full U when \(m \le n\), and a thin \(V^H\) equals the full one when \(m \ge n\). Entry points canonicalise Thin to All when the shapes coincide (canonical_jobu / canonical_jobvh).
Warning: a *_buffer_size and its run path must canonicalise identically, or the workspace is sized for a different computation.
LAPACK 'O' is absent. If added, append an enumerator so existing ordinals stay stable.
1.18.0
