From deebd05fc47f980d9f80155dc74d1f181ad11a68 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 20:22:57 -0700 Subject: [PATCH] Return INFO = 1 from xLASQ1 instead of stopping in XERBLA when SIGMX is a NaN xLASQ1 scales the bidiagonal matrix by its largest entry SIGMX before running dqds. When D or E contains a NaN that survives the MAX reductions computing SIGMX (with gfortran, a NaN in D(N)), SIGMX is a NaN and xLASCL rejects it by stopping the process in XERBLA. xBDSQR without singular vectors and xGESVD with JOBU = JOBVT = 'N' take this path; a NaN elsewhere in D returns INFO = 0 with NaN output. Test SIGMX with xISNAN and return INFO = 1, the mechanism #1382 uses for xLALSD and xBDSDC. xBDSQR treats every nonzero INFO from xLASQ1 as a request to finish with the QR algorithm, so it then reports the non-convergence on the NaN data through its own INFO like it does for other non-finite input. Finite input never takes the new branch. xERRBD gets the case as a regression test: a 4 by 4 bidiagonal with a NaN in D(N), passed to xBDSQR without 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. All four xERRBD files carry it, because the complex xBDSQR takes the same path through the real xLASQ1. Over 120 NaN/Inf cases of DBDSQR without vectors and DGESVD the parent stops in XERBLA 3 times and this branch never; every other case returns the same INFO 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/dlasq1.f | 14 ++++++++++++-- SRC/slasq1.f | 14 ++++++++++++-- TESTING/EIG/cerrbd.f | 32 +++++++++++++++++++++++++++++--- TESTING/EIG/derrbd.f | 30 ++++++++++++++++++++++++++++-- TESTING/EIG/serrbd.f | 30 ++++++++++++++++++++++++++++-- TESTING/EIG/zerrbd.f | 32 +++++++++++++++++++++++++++++--- 6 files changed, 138 insertions(+), 14 deletions(-) diff --git a/SRC/dlasq1.f b/SRC/dlasq1.f index 0b4903dd1..62bf85df8 100644 --- a/SRC/dlasq1.f +++ b/SRC/dlasq1.f @@ -81,7 +81,8 @@ *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: the algorithm failed -*> = 1, a split was marked by a positive value in E +*> = 1, a split was marked by a positive value in E, or +*> the input contains a NaN *> = 2, current block of Z not diagonalized after 100*N *> iterations (in inner while loop) On exit D and E *> represent a matrix with the same singular values @@ -132,7 +133,8 @@ SUBROUTINE DLASQ1( N, D, E, WORK, INFO ) * .. * .. External Functions .. DOUBLE PRECISION DLAMCH - EXTERNAL DLAMCH + LOGICAL DISNAN + EXTERNAL DISNAN, DLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT @@ -176,6 +178,14 @@ SUBROUTINE DLASQ1( N, D, E, WORK, INFO ) SIGMX = MAX( SIGMX, D( I ) ) 20 CONTINUE * +* A NaN in D or E can make SIGMX a NaN, which DLASCL would reject +* by stopping in XERBLA; report it through INFO instead. +* + IF( DISNAN( SIGMX ) ) THEN + INFO = 1 + RETURN + END IF +* * Copy D and E into WORK (in the Z format) and scale (squaring the * input data makes scaling by a power of the radix pointless). * diff --git a/SRC/slasq1.f b/SRC/slasq1.f index 4d8380874..42ae2c163 100644 --- a/SRC/slasq1.f +++ b/SRC/slasq1.f @@ -81,7 +81,8 @@ *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: the algorithm failed -*> = 1, a split was marked by a positive value in E +*> = 1, a split was marked by a positive value in E, or +*> the input contains a NaN *> = 2, current block of Z not diagonalized after 100*N *> iterations (in inner while loop) On exit D and E *> represent a matrix with the same singular values @@ -132,7 +133,8 @@ SUBROUTINE SLASQ1( N, D, E, WORK, INFO ) * .. * .. External Functions .. REAL SLAMCH - EXTERNAL SLAMCH + LOGICAL SISNAN + EXTERNAL SISNAN, SLAMCH * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, SQRT @@ -176,6 +178,14 @@ SUBROUTINE SLASQ1( N, D, E, WORK, INFO ) SIGMX = MAX( SIGMX, D( I ) ) 20 CONTINUE * +* A NaN in D or E can make SIGMX a NaN, which SLASCL would reject +* by stopping in XERBLA; report it through INFO instead. +* + IF( SISNAN( SIGMX ) ) THEN + INFO = 1 + RETURN + END IF +* * Copy D and E into WORK (in the Z format) and scale (squaring the * input data makes scaling by a power of the radix pointless). * diff --git a/TESTING/EIG/cerrbd.f b/TESTING/EIG/cerrbd.f index 71f49e80a..c4d21e0bb 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 RNAN, RONE * .. * .. Local Arrays .. REAL D( NMAX ), E( NMAX ), RW( 4*NMAX ) @@ -98,7 +99,7 @@ SUBROUTINE CERRBD( PATH, NUNIT ) COMMON / SRNAMC / SRNAMT * .. * .. Intrinsic Functions .. - INTRINSIC REAL + INTRINSIC REAL, SQRT * .. * .. Executable Statements .. * @@ -279,6 +280,29 @@ SUBROUTINE CERRBD( PATH, NUNIT ) $ INFO ) CALL CHKXER( 'CBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 +* +* CBDSQR without singular vectors must return when D contains a +* NaN, which the dqds path would otherwise pass to SLASCL as a +* scaling factor, instead of stopping in XERBLA. +* + RONE = ONE + RNAN = SQRT( -RONE ) + DO 30 J = 1, NMAX + D( J ) = REAL( J ) + E( J ) = ONE / REAL( J+1 ) + 30 CONTINUE + D( NMAX ) = RNAN + E( NMAX ) = ZERO + SRNAMT = 'CBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL CBDSQR( 'U', NMAX, 0, 0, 0, D, E, V, 1, U, 1, A, 1, 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 +317,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 a NaN', + $ ' ***' ) * RETURN * diff --git a/TESTING/EIG/derrbd.f b/TESTING/EIG/derrbd.f index 8f2b6603b..c6ce95a6f 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 RNAN, RONE * .. * .. Local Arrays .. INTEGER IQ( NMAX, NMAX ), IW( NMAX ) @@ -100,7 +101,7 @@ SUBROUTINE DERRBD( PATH, NUNIT ) COMMON / SRNAMC / SRNAMT * .. * .. Intrinsic Functions .. - INTRINSIC DBLE + INTRINSIC DBLE, SQRT * .. * .. Executable Statements .. * @@ -278,6 +279,29 @@ SUBROUTINE DERRBD( PATH, NUNIT ) CALL CHKXER( 'DBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 * +* DBDSQR without singular vectors must return when D contains a +* NaN, which the dqds path would otherwise pass to DLASCL as a +* scaling factor, instead of stopping in XERBLA. +* + RONE = ONE + RNAN = SQRT( -RONE ) + DO 30 J = 1, NMAX + D( J ) = DBLE( J ) + E( J ) = ONE / DBLE( J+1 ) + 30 CONTINUE + D( NMAX ) = RNAN + E( NMAX ) = ZERO + SRNAMT = 'DBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL DBDSQR( 'U', NMAX, 0, 0, 0, D, E, V, 1, U, 1, A, 1, W, + $ INFO ) + IF( LERR ) THEN + WRITE( NOUT, FMT = 9997 )'DBDSQR' + OK = .FALSE. + END IF + NT = NT + 1 +* * DBDSDC * SRNAMT = 'DBDSDC' @@ -369,6 +393,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 a NaN', + $ ' ***' ) * RETURN * diff --git a/TESTING/EIG/serrbd.f b/TESTING/EIG/serrbd.f index 683b8952e..861fbe033 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 RNAN, RONE * .. * .. Local Arrays .. INTEGER IQ( NMAX, NMAX ), IW( NMAX ) @@ -100,7 +101,7 @@ SUBROUTINE SERRBD( PATH, NUNIT ) COMMON / SRNAMC / SRNAMT * .. * .. Intrinsic Functions .. - INTRINSIC REAL + INTRINSIC REAL, SQRT * .. * .. Executable Statements .. * @@ -278,6 +279,29 @@ SUBROUTINE SERRBD( PATH, NUNIT ) CALL CHKXER( 'SBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 * +* SBDSQR without singular vectors must return when D contains a +* NaN, which the dqds path would otherwise pass to SLASCL as a +* scaling factor, instead of stopping in XERBLA. +* + RONE = ONE + RNAN = SQRT( -RONE ) + DO 30 J = 1, NMAX + D( J ) = REAL( J ) + E( J ) = ONE / REAL( J+1 ) + 30 CONTINUE + D( NMAX ) = RNAN + E( NMAX ) = ZERO + SRNAMT = 'SBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL SBDSQR( 'U', NMAX, 0, 0, 0, D, E, V, 1, U, 1, A, 1, W, + $ INFO ) + IF( LERR ) THEN + WRITE( NOUT, FMT = 9997 )'SBDSQR' + OK = .FALSE. + END IF + NT = NT + 1 +* * SBDSDC * SRNAMT = 'SBDSDC' @@ -369,6 +393,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 a NaN', + $ ' ***' ) * RETURN * diff --git a/TESTING/EIG/zerrbd.f b/TESTING/EIG/zerrbd.f index 08e7821f5..8bb9234ff 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 RNAN, RONE * .. * .. Local Arrays .. DOUBLE PRECISION D( NMAX ), E( NMAX ), RW( 4*NMAX ) @@ -98,7 +99,7 @@ SUBROUTINE ZERRBD( PATH, NUNIT ) COMMON / SRNAMC / SRNAMT * .. * .. Intrinsic Functions .. - INTRINSIC DBLE + INTRINSIC DBLE, SQRT * .. * .. Executable Statements .. * @@ -279,6 +280,29 @@ SUBROUTINE ZERRBD( PATH, NUNIT ) $ INFO ) CALL CHKXER( 'ZBDSQR', INFOT, NOUT, LERR, OK ) NT = NT + 8 +* +* ZBDSQR without singular vectors must return when D contains a +* NaN, which the dqds path would otherwise pass to DLASCL as a +* scaling factor, instead of stopping in XERBLA. +* + RONE = ONE + RNAN = SQRT( -RONE ) + DO 30 J = 1, NMAX + D( J ) = DBLE( J ) + E( J ) = ONE / DBLE( J+1 ) + 30 CONTINUE + D( NMAX ) = RNAN + E( NMAX ) = ZERO + SRNAMT = 'ZBDSQR' + INFOT = 0 + LERR = .FALSE. + CALL ZBDSQR( 'U', NMAX, 0, 0, 0, D, E, V, 1, U, 1, A, 1, 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 +317,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 a NaN', + $ ' ***' ) * RETURN *