From 941f78cb750a116c7d8dbacd0f2f72b3d3877e15 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 20:21:32 -0700 Subject: [PATCH] Let an exact zero always split the matrix in xBDSQR so that an Inf cannot hang it xBDSQR with singular vectors never returns for a bidiagonal matrix with an Inf in one of its first two diagonal entries, and neither do xBDSDC with COMPQ = 'I' and the drivers xGESVD and xGESDD built on them. The relative-accuracy threshold THRESH comes from a recurrence that evaluates Inf/Inf for such input and becomes NaN; the block search IF( ABSE.LE.THRESH ) then never recognises a split, while the convergence test zeroes E( M-1 ) and jumps back to the top of the loop without advancing the iteration counter. The two chase each other forever. Without vectors xBDSQR calls xLASQ1 instead and returns. Make an exact zero in E a split regardless of THRESH. For a finite THRESH, which is at least MAXITR*N*N*UNFL > 0, the extra clause is implied by the existing test, so nothing changes for finite input; for a NaN threshold every path back to the loop head now shrinks the block or advances ITER, and the routine returns within its iteration bound like it does for other non-finite input. xERRBD gets the case as a regression test: a 4 by 4 bidiagonal with an infinite D(1), passed to xBDSQR with both sets of singular vectors, which must return without calling XERBLA. The error-exit tests are where it belongs, since the routine returns no meaningful output for such a matrix and the point of the test is that it returns at all; on the parent, all four xeigtst binaries stop responding at that call. Over 480 Inf/NaN cases (the four xBDSQR, DBDSDC, DGESVD, DGESDD, ZGESVD; five positions; n = 2, 3, 5, 8) the parent hangs in 66, this branch returns in all of them. 160 finite bidiagonal matrices give bit-identical singular values and vectors on both. The full LAPACK test suite passes: 5441901 LAPACK tests, 0 numerical errors, 0 other errors, the same totals as the parent. Co-Authored-By: Claude Fable 5.1 --- SRC/cbdsqr.f | 6 +++++- SRC/dbdsqr.f | 6 +++++- SRC/sbdsqr.f | 6 +++++- SRC/zbdsqr.f | 6 +++++- TESTING/EIG/cerrbd.f | 35 +++++++++++++++++++++++++++++++++-- TESTING/EIG/derrbd.f | 33 ++++++++++++++++++++++++++++++++- TESTING/EIG/serrbd.f | 33 ++++++++++++++++++++++++++++++++- TESTING/EIG/zerrbd.f | 35 +++++++++++++++++++++++++++++++++-- 8 files changed, 150 insertions(+), 10 deletions(-) diff --git a/SRC/cbdsqr.f b/SRC/cbdsqr.f index fd6d4731b..df0ee5076 100644 --- a/SRC/cbdsqr.f +++ b/SRC/cbdsqr.f @@ -455,7 +455,11 @@ SUBROUTINE CBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, ABSE = ABS( E( LL ) ) IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) $ D( LL ) = ZERO - IF( ABSE.LE.THRESH ) +* An exact zero always marks a split, also when THRESH is NaN +* because the input contains an Inf; otherwise the convergence +* tests below zero E( M-1 ) and return here without the split +* ever being found. + IF( ABSE.LE.THRESH .OR. ABSE.EQ.ZERO ) $ GO TO 80 SMAX = MAX( SMAX, ABSS, ABSE ) 70 CONTINUE diff --git a/SRC/dbdsqr.f b/SRC/dbdsqr.f index 99a12aa79..e9dd2ae2b 100644 --- a/SRC/dbdsqr.f +++ b/SRC/dbdsqr.f @@ -462,7 +462,11 @@ SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, ABSE = ABS( E( LL ) ) IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) $ D( LL ) = ZERO - IF( ABSE.LE.THRESH ) +* An exact zero always marks a split, also when THRESH is NaN +* because the input contains an Inf; otherwise the convergence +* tests below zero E( M-1 ) and return here without the split +* ever being found. + IF( ABSE.LE.THRESH .OR. ABSE.EQ.ZERO ) $ GO TO 80 SMAX = MAX( SMAX, ABSS, ABSE ) 70 CONTINUE diff --git a/SRC/sbdsqr.f b/SRC/sbdsqr.f index 866b13295..9c06037a1 100644 --- a/SRC/sbdsqr.f +++ b/SRC/sbdsqr.f @@ -462,7 +462,11 @@ SUBROUTINE SBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, ABSE = ABS( E( LL ) ) IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) $ D( LL ) = ZERO - IF( ABSE.LE.THRESH ) +* An exact zero always marks a split, also when THRESH is NaN +* because the input contains an Inf; otherwise the convergence +* tests below zero E( M-1 ) and return here without the split +* ever being found. + IF( ABSE.LE.THRESH .OR. ABSE.EQ.ZERO ) $ GO TO 80 SMAX = MAX( SMAX, ABSS, ABSE ) 70 CONTINUE diff --git a/SRC/zbdsqr.f b/SRC/zbdsqr.f index ec67efeb0..74eb8776a 100644 --- a/SRC/zbdsqr.f +++ b/SRC/zbdsqr.f @@ -453,7 +453,11 @@ SUBROUTINE ZBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, ABSE = ABS( E( LL ) ) IF( TOL.LT.ZERO .AND. ABSS.LE.THRESH ) $ D( LL ) = ZERO - IF( ABSE.LE.THRESH ) +* An exact zero always marks a split, also when THRESH is NaN +* because the input contains an Inf; otherwise the convergence +* tests below zero E( M-1 ) and return here without the split +* ever being found. + IF( ABSE.LE.THRESH .OR. ABSE.EQ.ZERO ) $ GO TO 80 SMAX = MAX( SMAX, ABSS, ABSE ) 70 CONTINUE diff --git a/TESTING/EIG/cerrbd.f b/TESTING/EIG/cerrbd.f index 71f49e80a..ed4fa15db 100644 --- a/TESTING/EIG/cerrbd.f +++ b/TESTING/EIG/cerrbd.f @@ -68,12 +68,13 @@ SUBROUTINE CERRBD( PATH, NUNIT ) * .. Parameters .. INTEGER NMAX, LW PARAMETER ( NMAX = 4, LW = NMAX ) - REAL ONE - PARAMETER ( ONE = 1.0E+0 ) + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. CHARACTER*2 C2 INTEGER I, INFO, J, NT + REAL RZERO * .. * .. Local Arrays .. REAL D( NMAX ), E( NMAX ), RW( 4*NMAX ) @@ -279,6 +280,34 @@ SUBROUTINE CERRBD( PATH, NUNIT ) $ INFO ) CALL CHKXER( 'CBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 +* +* CBDSQR with singular vectors must return when D contains an +* infinity, which makes its convergence threshold a NaN, +* instead of iterating forever. +* + RZERO = ZERO + DO 40 J = 1, NMAX + D( J ) = REAL( J ) + E( J ) = ONE / REAL( J+1 ) + DO 30 I = 1, NMAX + U( I, J ) = ZERO + V( I, J ) = ZERO + 30 CONTINUE + U( J, J ) = ONE + V( J, J ) = ONE + 40 CONTINUE + D( 1 ) = ONE / RZERO + E( NMAX ) = ZERO + SRNAMT = 'CBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL CBDSQR( 'U', NMAX, NMAX, NMAX, 0, D, E, V, NMAX, U, + $ NMAX, A, NMAX, RW, INFO ) + IF( LERR ) THEN + WRITE( NOUT, FMT = 9997 )'CBDSQR' + OK = .FALSE. + END IF + NT = NT + 1 END IF * * Print a summary line. @@ -293,6 +322,8 @@ SUBROUTINE CERRBD( PATH, NUNIT ) $ I3, ' tests done)' ) 9998 FORMAT( ' *** ', A3, ' routines failed the tests of the error ', $ 'exits ***' ) + 9997 FORMAT( ' *** ', A6, ' called XERBLA for a matrix with an ', + $ 'infinity ***' ) * RETURN * diff --git a/TESTING/EIG/derrbd.f b/TESTING/EIG/derrbd.f index 8f2b6603b..af8a21b23 100644 --- a/TESTING/EIG/derrbd.f +++ b/TESTING/EIG/derrbd.f @@ -67,13 +67,14 @@ SUBROUTINE DERRBD( PATH, NUNIT ) * * .. Parameters .. INTEGER NMAX, LW - PARAMETER ( NMAX = 4, LW = NMAX ) + PARAMETER ( NMAX = 4, LW = 4*NMAX ) DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. CHARACTER*2 C2 INTEGER I, INFO, J, NS, NT + DOUBLE PRECISION RZERO * .. * .. Local Arrays .. INTEGER IQ( NMAX, NMAX ), IW( NMAX ) @@ -278,6 +279,34 @@ SUBROUTINE DERRBD( PATH, NUNIT ) CALL CHKXER( 'DBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 * +* DBDSQR with singular vectors must return when D contains an +* infinity, which makes its convergence threshold a NaN, +* instead of iterating forever. +* + RZERO = ZERO + DO 40 J = 1, NMAX + D( J ) = DBLE( J ) + E( J ) = ONE / DBLE( J+1 ) + DO 30 I = 1, NMAX + U( I, J ) = ZERO + V( I, J ) = ZERO + 30 CONTINUE + U( J, J ) = ONE + V( J, J ) = ONE + 40 CONTINUE + D( 1 ) = ONE / RZERO + E( NMAX ) = ZERO + SRNAMT = 'DBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL DBDSQR( 'U', NMAX, NMAX, NMAX, 0, D, E, V, NMAX, U, + $ NMAX, A, NMAX, W, INFO ) + IF( LERR ) THEN + WRITE( NOUT, FMT = 9997 )'DBDSQR' + OK = .FALSE. + END IF + NT = NT + 1 +* * DBDSDC * SRNAMT = 'DBDSDC' @@ -369,6 +398,8 @@ SUBROUTINE DERRBD( PATH, NUNIT ) $ ' (', I3, ' tests done)' ) 9998 FORMAT( ' *** ', A3, ' routines failed the tests of the error ', $ 'exits ***' ) + 9997 FORMAT( ' *** ', A6, ' called XERBLA for a matrix with an ', + $ 'infinity ***' ) * RETURN * diff --git a/TESTING/EIG/serrbd.f b/TESTING/EIG/serrbd.f index 683b8952e..0d0202e83 100644 --- a/TESTING/EIG/serrbd.f +++ b/TESTING/EIG/serrbd.f @@ -67,13 +67,14 @@ SUBROUTINE SERRBD( PATH, NUNIT ) * * .. Parameters .. INTEGER NMAX, LW - PARAMETER ( NMAX = 4, LW = NMAX ) + PARAMETER ( NMAX = 4, LW = 4*NMAX ) REAL ZERO, ONE PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) * .. * .. Local Scalars .. CHARACTER*2 C2 INTEGER I, INFO, J, NS, NT + REAL RZERO * .. * .. Local Arrays .. INTEGER IQ( NMAX, NMAX ), IW( NMAX ) @@ -278,6 +279,34 @@ SUBROUTINE SERRBD( PATH, NUNIT ) CALL CHKXER( 'SBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 * +* SBDSQR with singular vectors must return when D contains an +* infinity, which makes its convergence threshold a NaN, +* instead of iterating forever. +* + RZERO = ZERO + DO 40 J = 1, NMAX + D( J ) = REAL( J ) + E( J ) = ONE / REAL( J+1 ) + DO 30 I = 1, NMAX + U( I, J ) = ZERO + V( I, J ) = ZERO + 30 CONTINUE + U( J, J ) = ONE + V( J, J ) = ONE + 40 CONTINUE + D( 1 ) = ONE / RZERO + E( NMAX ) = ZERO + SRNAMT = 'SBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL SBDSQR( 'U', NMAX, NMAX, NMAX, 0, D, E, V, NMAX, U, + $ NMAX, A, NMAX, W, INFO ) + IF( LERR ) THEN + WRITE( NOUT, FMT = 9997 )'SBDSQR' + OK = .FALSE. + END IF + NT = NT + 1 +* * SBDSDC * SRNAMT = 'SBDSDC' @@ -369,6 +398,8 @@ SUBROUTINE SERRBD( PATH, NUNIT ) $ ' (', I3, ' tests done)' ) 9998 FORMAT( ' *** ', A3, ' routines failed the tests of the error ', $ 'exits ***' ) + 9997 FORMAT( ' *** ', A6, ' called XERBLA for a matrix with an ', + $ 'infinity ***' ) * RETURN * diff --git a/TESTING/EIG/zerrbd.f b/TESTING/EIG/zerrbd.f index 08e7821f5..d18654d12 100644 --- a/TESTING/EIG/zerrbd.f +++ b/TESTING/EIG/zerrbd.f @@ -68,12 +68,13 @@ SUBROUTINE ZERRBD( PATH, NUNIT ) * .. Parameters .. INTEGER NMAX, LW PARAMETER ( NMAX = 4, LW = NMAX ) - DOUBLE PRECISION ONE - PARAMETER ( ONE = 1.0D+0 ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. CHARACTER*2 C2 INTEGER I, INFO, J, NT + DOUBLE PRECISION RZERO * .. * .. Local Arrays .. DOUBLE PRECISION D( NMAX ), E( NMAX ), RW( 4*NMAX ) @@ -279,6 +280,34 @@ SUBROUTINE ZERRBD( PATH, NUNIT ) $ INFO ) CALL CHKXER( 'ZBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 +* +* ZBDSQR with singular vectors must return when D contains an +* infinity, which makes its convergence threshold a NaN, +* instead of iterating forever. +* + RZERO = ZERO + DO 40 J = 1, NMAX + D( J ) = DBLE( J ) + E( J ) = ONE / DBLE( J+1 ) + DO 30 I = 1, NMAX + U( I, J ) = ZERO + V( I, J ) = ZERO + 30 CONTINUE + U( J, J ) = ONE + V( J, J ) = ONE + 40 CONTINUE + D( 1 ) = ONE / RZERO + E( NMAX ) = ZERO + SRNAMT = 'ZBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL ZBDSQR( 'U', NMAX, NMAX, NMAX, 0, D, E, V, NMAX, U, + $ NMAX, A, NMAX, RW, INFO ) + IF( LERR ) THEN + WRITE( NOUT, FMT = 9997 )'ZBDSQR' + OK = .FALSE. + END IF + NT = NT + 1 END IF * * Print a summary line. @@ -293,6 +322,8 @@ SUBROUTINE ZERRBD( PATH, NUNIT ) $ I3, ' tests done)' ) 9998 FORMAT( ' *** ', A3, ' routines failed the tests of the error ', $ 'exits ***' ) + 9997 FORMAT( ' *** ', A6, ' called XERBLA for a matrix with an ', + $ 'infinity ***' ) * RETURN *