Add DBDSVDMR3: MRRR-based bidiagonal SVD (Golub-Kahan TGK) with DBDSDC safety fallback

[saisuryadv]

Update (naming): per review, the driver is now DBDSVDMR3 (was DBDSVR) and the computational routine DBDSVDMR3_WORK (was DBDSVDMR3); the battery below was re-run after the rename with identical results.

Add DBDSVDMR3: an MRRR-based bidiagonal SVD driver

Summary

DBDSVDMR3 computes the SVD of an n-by-n bidiagonal matrix by the MRRR algorithm applied to the Golub-Kahan tridiagonal, with the same argument list as DBDSVDX (plus LWORK/LIWORK) and a built-in safety net: every output is audited for nonfiniteness and, on audit or nonzero INFO, inputs are restored and DBDSDC recomputes, with the execution path reported in IWORK(2) (1 native / 2 fallback-on-INFO / 3 fallback-on-audit). Motivation: DBDSVDX is 7-228x slower than DBDSVDMR3 on a 16-variety battery and returns unorthogonalized vectors with INFO=0 on clustered spectra (ratios to 2.2e14); DBDSDC remains the deflation-friendly workhorse, and DBDSVDMR3 complements it with the best worst-case time of the four solvers tested.

Integration pattern (mirrors the DBDSVDX integration)

• SRC/: dbdsvdmr3.f (driver), dbdsvdmr3_work.f, dbdtgk.f, dlar1v_tgk.f, dlarrf_tgk.f, dlarrv_tgk.f, dlarrb_tgk.f (bundled bisection renamed from DLARRB - its extended signature collides with SRC/dlarrb.f; linking the stock routine segfaults).
• SRC/Makefile, SRC/CMakeLists.txt: DLASRC lists extended.
• TESTING/EIG/dchkbd.f: tests 35-49, an exact programmatic mirror of the DBDSVDX tests 20-34 (RANGE='A','I','V' x residual, U/VT orthogonality, ordering, values-vs-vectors agreement).
• TESTING/EIG/derrbd.f: 12 error-exit checks.

Validation on current master (svd.in)

• All tests for DBD routines passed the threshold ( 14820 tests run)
• All tests for DBD drivers passed the threshold ( 14820 tests run)
• All tests for DBD routines passed the threshold ( 14820 tests run)
• All tests for DBD drivers passed the threshold ( 14820 tests run)
• All tests for DBD routines passed the threshold ( 14820 tests run)
• All tests for DBD drivers passed the threshold ( 14820 tests run)
• All tests for DBD routines passed the threshold ( 14820 tests run)
• All tests for DBD drivers passed the threshold ( 14820 tests run)
• All tests for DBD routines passed the threshold ( 14820 tests run)
• All tests for DBD drivers passed the threshold ( 14820 tests run) - includes the 15 new DBDSVDMR3 tests across all sizes/types.
• DBD routines passed the tests of the error exits ( 55 tests done) - includes the 12 new checks.

Accuracy on a harder 16-variety battery (separate DCHKBSVR harness)

Sizes 2-100, DBDT04/DORT01 in n.ulp units, threshold 20; identical 240-case RANGE='A' basis; DBDSVDMR3/DBDSVDX additionally exercised on I/V modes:

Method	n	exact-0 share	share > 20	p50	p90	p99	max
DBDSVDMR3 (ours)	240	6%	2.1%	0.20	6.38	25.76	30.20
DBDSDC	240	45%	0.0%	0.01	0.76	1.41	1.53
DBDSVDX	240	32%	0.0%	0.03	0.91	3.41	7.80
MR3GK (Willems & Lang)	240	44%	0.0%	0.02	1.66	4.65	14.50

DBDSVDMR3: 1195/1200 residuals below threshold; the exceedances (max 30.2, i.e. <=1.6x threshold) are confined to Cholesky-lifted geometric spectra where the construction fuses ~80% of sigma into one 1e-13-tight cluster (child-representation selection in dlarrf_tgk; fix path = block factorizations, validated externally). DBDSVDX on the same battery fails variety 8 at up to 2.2e14 with INFO=0 in its I/V modes. Wilkinson-class zero-diagonal inputs route to fallback path 3 by design.

Performance (same battery; ratios are the reliable objects)

N=2000, ms per call (unpinned VM; median 39% absolute drift measured, 11% on ratios; a pinned EPYC 9454 confirms all rankings and reaches 65x over DBDSDC on Clement at n=3000):

Ty	Variety	DBDSVDMR3 (ours)	DBDSDC	DBDSVDX	MR3GK (Willems & Lang)
1	zero	39.0	3.1	12.3	5.1
2	identity	50.7	6.4	16.5	5.2
3	arith bidiag	296.5	3,451.0†	27,200.0†	8,691.4†
4	geom bidiag	255.6	206.3	23,634.0†	2,968.5†
5	log-dist	88.8	41.2	25,042.0†	287.6
6	arith tri-chol	700.7	1,665.4†	8,585.4†	8,102.7†
7	geom tri-chol	1,124.1†	177.8	19,765.0†	3,276.3†
8	clust tri-chol	85.3	35.0	18.4	17.0
9	arith sqrt-ovfl	278.3	3,715.1†	27,154.0†	8,899.9†
10	arith sqrt-unfl	352.0	3,466.9†	28,675.0†	8,537.3†
11	arith diag	75.5	6.1	15.4	5.0
12	geom diag	72.2	5.9	16.0	5.0
13	clust diag	59.4	6.3	15.0	5.0
14	diag sqrt-ovfl	74.0	7.5	16.3	5.0
15	diag sqrt-unfl	70.7	5.8	19.6	5.0
16	SPD geom-chol	830.6	292.2	20,374.0†	4,169.3†

Pooled t/n^2 percentiles over 80 (variety, size) cases:

Method	p25	p50	p75	p90	max
DBDSVDMR3 (ours)	8.6	18.7	63.8	176.5	281.0
DBDSDC	1.6	4.9	83.9	288.3	928.8
DBDSVDX	5.7	390.4	1,389.9	5,139.6	7,168.8
MR3GK (Willems & Lang)	1.3	38.8	1,237.0	1,955.9	2,241.7

DBDSVDMR3 wins exactly the deflation-free spectra (9.8-13.4x over DBDSDC at n=2000), concedes the deflation-friendly ones (2.2-6.3x), has the smallest maximum of any solver, and the best suite total (2.9x over DBDSDC at n=2000, growing with n).

Known limitations / review questions

Double precision only (SBDSVDMR3 to follow); RANGE='I'/'V' currently post-filters a full computation (index pushdown into dlarrv_tgk is the planned follow-up); the two clustered-spectrum exceedance classes above; LAPACKE bindings not included; the fallback introduces a dependency on DBDSDC.

[mgates3]

Interesting approach to an SVD with MRRR.

I think names bdsvr and bdsvdmr3 are confusing. If users are looking for an SVD with MRRR, they will likely find the low-level bdsvdmr3 and could miss that the high-level bdsvr is even an SVD. The archaic 6-char limit produced some inscrutable abbreviations; let's not perpetuate the trend. (Note confusion with gesvj and gejsv, which have rather unfortunate names; gesvdj or gesvd_jacobi would be better.)

I'm not sure the best way to distinguish them, but I think both should start bdsvd. I would probably use something like bdsvdmr3 for the high-level driver (currently bdsvr) and bdsvdmr3_work for the low-level computational routine (currently bdsvdmr3). There's not really an LAPACK naming convention for high-level vs. low-level routines. Or bdsvdmrrr or bdsvd_mrrr; the 3 saves only 1 char and is less recognizable.

If the "classic" SVD with MRRR -- which has had numerical stability issues, hence never made it to LAPACK -- were ever added to LAPACK, would it replace these names, or how would it be distinguished? Something to consider.

[saisuryadv]

Thanks @mgates3 — agreed on all points, and I've adopted exactly the pair you suggested. The routines are now:

• DBDSVR → DBDSVDMR3 (high-level driver)
• DBDSVDMR3 → DBDSVDMR3_WORK (computational routine)

so both start bdsvd and the driver is discoverable as an SVD. The internal MRRR kernels keep their dlarr*_tgk names, which describe the Golub–Kahan adaptation rather than user-facing functionality. The testers, error-exit strings (derrbd now checks DBDSVDMR3), and both build lists were updated, and the full svd.in battery was re-run after the rename: all DBD ratios pass the threshold across all five data sets, error exits 55/55. PR title and description updated to match.

On your forward-compatibility question: the driver is engine-agnostic — its interface (the DBDSVDX argument list plus LWORK/LIWORK) makes no commitment to how the decomposition is computed. If a "classic" structure-preserving MRRR SVD were added later, the cleanest path would be replacing or augmenting _WORK behind the same DBDSVDMR3 name, since the accuracy characteristics would only improve; if the community preferred both engines exposed, a distinguishing suffix on the worker level keeps the single driver name intact. Happy to adjust further if the team prefers bdsvdmrrr spelled out.

[rmlarsen]

Thanks for this contribution. One thing up front: I am not a LAPACK maintainer, so please read what follows as review and advice from an independent contributor with a general interest in SVD algorithms, not as a statement of what the project requires. An MRRR-based bidiagonal SVD with dqds singular values is an interesting direction, and the wall-clock advantage on deflation-free spectra is real: at n = 2000 I measure 0.39 s vs 5.6 s for DBDSDC on the all-ones bidiagonal and 0.38 s vs 3.0 s on a Lipshitz-type spectrum, and your observation that DBDSVDX loses orthogonality on tight clusters reproduces too.

I built the branch (6d0ac24, gfortran 13.3, Release), ran the full suite (215/215; DBD 14820 tests and 55 error exits pass), and then exercised the driver on a battery of generated bidiagonals (random, graded, glued blocks, tight clusters, zero diagonals, extreme scaling, NaN/Inf) with valgrind, AddressSanitizer and -fcheck=all. The probe program, scripts and raw logs are in this gist: https://gist.github.com/rmlarsen/a50e1563c7a4b0e861c0cd4d1942a998. Below is what came out, roughly in order of importance, with file:line pointers.

1. RANGE='I' selects different singular values than DBDSVDX. dbdsvdmr3.f:489-491 uses IL/IU as indices into the ascending spectrum, so IL = IU = 1 returns the smallest singular value; DBDSVDX returns the largest (it targets -s(IL)..-s(IU) on the TGK). For diag(5,4,3,2,1) the two routines return 1.0 and 5.0 respectively. Since the driver advertises DBDSVDX's argument list, I think it should map ILO = N-IU+1, IHI = N-IL+1. dchkbd cannot catch this because tests 25-29 / 40-44 only compare 'V' against 'N' for the same IL, IU.

2. Heap overflow in the DSTEIN fallback of DLARRV_TGK. With LIWORK = 20N, the driver leaves DLARRV_TGK 12N integers, but the DSTEIN call at dlarrv_tgk.f:306-313 indexes up to IINDWK + 2*IN + NEWSIZ = 14N + NEWSIZ (IN is the TGK order 2N). ASan reports dstein.f:225 writing 0 bytes past the IWORK block; a NaN anywhere in D reaches this path and the process aborts in malloc (probe nan 20 5, nan 200 100). Two side notes on that path: at the root (REP='T') WORK(INDLD)/WORK(INDLLD) are identically zero (lines 205-211), so DSTEIN there would factor the zero matrix; and the driver then reports DBDSDC's INFO (19 in my run) although the documentation only says "every backend failed".

3. Out-of-bounds read at dlarrv_tgk.f:421. IF( KTOT2.GT.0 .AND. IWORK( KTOT2 ).EQ.-1 ) evaluates IWORK(0) when KTOT2 = 0 since Fortran does not short-circuit. -fcheck=all flags it for ordinary inputs (probe ones 3, zerod 100 3, tinye 100 1e-16). Benign in effect, but it needs a nested IF.

4. A cluster can be dropped silently (latent). At dlarrv_tgk.f:281-285, when DLARRF_TGK returns INFO = 1 and MINRGP >= MGSTOL, NOMGS = .FALSE. is set inside the branch that NOMGS guards, so the ELSE IF( NEWSIZ.GT.1 ) path at line 321 is never re-entered (with MGSSIZ = 1 it is unreachable anyway). The cluster is neither processed nor counted, the loop spins to NDEPTH > 4N, and the routine returns INFO = 0 with zero columns; the driver's audit only looks for NaN/Inf. LAPACK 3.0's DLARRV handled this via MGSCLS = .TRUE.. I could not trigger it, but a nonzero INFO there would be safer.

5. Orthogonality on clustered and glued spectra is 10-100x worse than DBDSDC and exceeds the dchkbd threshold (50) on at least one input. orthU = |I - U^T U| / (n ulp):

probe args	DBDSVDMR3	DBDSDC	DBDSVDX
glue 200 20 1e-14 L (10 glued 20-blocks, lower)	54.5	0.17	0.17
glue 100 10 1e-15	36.4	0.18	0.20
ones 1000 (d = e = 1)	34.1	0.36	1.4
lipshitz 500 1e-12	19.4	0.27	-
cluster 1000 1e-7 (d = 1, e = 1e-7)	13.7	0.33	0.18
rand 100	3.0	0.28	0.38

Residuals are good throughout and column norms are 1 to ~1e-15, so this is orthogonality within clusters (the child-representation issue you mention), not the TGK norm-coupling effect. Given that the PR's own battery reports p99 = 25.8 and max = 30.2 against a threshold of 20, it would be good to understand whether the block-factorization children you mention are feasible before this lands.

6. The split threshold is absolute. dbdsvdmr3_work.f:236-244 zeroes every TGK off-diagonal (d_i as well as e_i) below eps * max|B|, so small singular values that dqds would deliver to full relative accuracy come back as exactly 0. For d = (1, 1e-17), e = (0.5): DBDSVDMR3 gives S(2) = 0, DBDSQR and DBDSVDX give 8.944e-18. graded 100 30 (values only) has relative error 1.0 against DBDSQR for those entries. DBDSDC behaves the same way, so this may be an acceptable design choice, but it should be documented, or a relative criterion as in DBDSVDX used, since relative accuracy is the usual reason to reach for an MRRR/dqds solver.

7. NaN/Inf. Inf in D(1) or D(5) hangs (in DBDSQR inside the DBDSDC fallback; that is master's bug, see #1390), NaN in D(N) aborts through XERBLA in DLASQ1 (see #1387), and Inf in D(N) returns INFO = 0 with S(1) = Inf and NaN vectors via path 3, because the finiteness audit is only applied to the MR3 result, not to the fallback. Applying the audit to both would make the "safety net" story consistent.

8. dchkbd.f mirror. The per-type init DO 30 J = 1, 34 was not extended to 49, so RESULT(35:49) are undefined for the skipped types; the new skip block zeroes RESULT( 20:34 ) instead of 35:49; test 48 assigns RESULT( 43 ) (so it never fires and can overwrite test 43); DLAHD2 has no text for 35-49; and the test needs 7N^2 + 42N doubles, which exceeds dchkee's LWORK for N > 109 (svd.in stops at 40, so it is not visible today).

9. Smaller items, mostly documentation and integration:

• dbdsvdmr3.f: \brief DBDSVR, the Purpose text and \ingroup bdsvr are stale after the rename; line 508 cites /tmp/lapack-ref/SRC/dbdsvdx.f line 747; line 221 says "offline PR".
• The other six files have no Doxygen blocks, and comments refer to "the advisor's stegr_ID" and "thesis eq. 3.1.12 / Algorithm 3.3.1" without a citation; dlarrb_tgk.f still carries a "version TBA, November 11, 2003" header. The kernels descend from the LAPACK 3.0 (1999) DLARRV/DLAR1V/DLARRB/DLARRF rather than the current ones; worth stating explicitly, since the 3.1 rewrite fixed a number of robustness issues.
• dlarrv_tgk.f:69 is 75 columns wide, so the fixed-form EXTERNAL statement is silently truncated.
• Extra undocumented global routines bundled in the files (DBSORT, DLAR1V2_TGK, DLARRBTGK_STURM4, DLARRBTGK_STURM_SLOW); the hand-interleaved two- and four-lane recurrences are unusual for reference code and I would keep them out of a first version.
• dbdsvdmr3_work.f sets ISEED and PERTK but never uses them; IWORK(2) as an output channel is unusual for LAPACK.
• dbdsvdmr3.f:459-462 passes WORK(IWRK) as both Q and WORK (and IWORK as both IQ and IWORK) to DBDSDC, and for JOBZ='N' the engine gets U, VT and WORK all at WORK(IWRK): harmless today but nonconforming aliasing.
• Single precision and LAPACKE are usually expected alongside a new double routine.

Happy to re-run the battery on an updated branch; the gist has the exact commands.