From b9ed114c5a94ef85630ce6c3f1a90e49dda9b5be Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Thu, 10 Sep 2026 01:11:18 -0700 Subject: [PATCH 1/2] Keep a NaN inside its block in xSTEDC and report it through INFO xSTEDC splits the tridiagonal matrix into independent blocks where an off-diagonal entry is negligible against its diagonal neighbours: TINY = EPS*SQRT( ABS( D( FINISH ) ) )*SQRT( ABS( D( FINISH+1 ) ) ) IF( ABS( E( FINISH ) ).GT.TINY ) THEN extend the block Both comparisons are false for a NaN, so a NaN off-diagonal entry split the matrix and was dropped, and a NaN diagonal entry made TINY a NaN, ended the block before it and left it as a 1 by 1 block. The routine then returned INFO = 0 with a finite spectrum for a matrix with a NaN off the diagonal, or with the NaN reported as an eigenvalue and the coupling to its neighbours ignored. The other tridiagonal solvers, xSTEQR, xSTERF, xSTEBZ and xSTEMR, keep a NaN in its block. Extend the block unless the off-diagonal entry is known to be small, so that a NaN stays with its neighbours. A block that reaches the divide and conquer recursion is scaled by its max-norm with xLASCL, which stops in XERBLA for a NaN, so the norm is tested first and a NaN block is reported as a failure on that block, INFO = START*(N+1) + FINISH, the encoding the QR fallback already uses. A block small enough for xSTEQR gets the NaN through the QR iteration, which returns INFO > 0 or, for a 2 by 2 block, NaN eigenvalues. The same code sits in cSTEDC and zSTEDC; their COMPZ = 'I' path calls the real routine. xCHKST gets the case as a regression test for the divide and conquer path, which the sizes in sep.in (N <= 20, below SMLSIZ = 25) never reach: after the size and type loops it calls xSTEDC with COMPZ = 'I' and 'V' on an N = LDU tridiagonal matrix, once with a NaN in the middle of E and once with a NaN at the end of D, and reports INFO = 0 as a failure. On the parent commit all four calls return INFO = 0 in every precision; here they return INFO > 0. Over a NaN and Inf sweep of xSTEDC and xSTEVD (six positions, n = 1 to 64, every COMPZ and JOBZ, real and complex) the parent returned INFO = 0 for a NaN matrix with n >= 3 in 90 cases (DSTEDC 36, ZSTEDC 36, DSTEVD 18) and this branch in none of the 280 such cases; every finite case is bit-identical. The full LAPACK test suite passes: 0 numerical errors, 0 other errors, 80 tests more than the parent from the new calls. The regression test fails on the parent and passes with the fix, and the reproducer prints the same before and after output, with gfortran 13 (x86-64 Release and Debug with -fcheck=all, and under QEMU on aarch64, ppc64le, s390x and riscv64), flang-19 and Intel ifx 2025.3. Co-Authored-By: Claude Fable 5.1 --- SRC/cstedc.f | 13 ++++++++--- SRC/dstedc.f | 13 ++++++++--- SRC/sstedc.f | 13 ++++++++--- SRC/zstedc.f | 13 ++++++++--- TESTING/EIG/cchkst.f | 54 ++++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/dchkst.f | 53 +++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/schkst.f | 53 +++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/zchkst.f | 54 ++++++++++++++++++++++++++++++++++++++++++++ 8 files changed, 254 insertions(+), 12 deletions(-) diff --git a/SRC/cstedc.f b/SRC/cstedc.f index 2a912cb96..bd76a9173 100644 --- a/SRC/cstedc.f +++ b/SRC/cstedc.f @@ -230,10 +230,10 @@ SUBROUTINE CSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, REAL EPS, ORGNRM, P, TINY * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL SISNAN, LSAME INTEGER ILAENV REAL SLAMCH, SLANST, SROUNDUP_LWORK - EXTERNAL ILAENV, LSAME, SLAMCH, SLANST, + EXTERNAL ILAENV, LSAME, SISNAN, SLAMCH, SLANST, $ SROUNDUP_LWORK * .. * .. External Subroutines .. @@ -395,7 +395,10 @@ SUBROUTINE CSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, IF( FINISH.LT.N ) THEN TINY = EPS*SQRT( ABS( D( FINISH ) ) )* $ SQRT( ABS( D( FINISH+1 ) ) ) - IF( ABS( E( FINISH ) ).GT.TINY ) THEN +* A NaN in D or E must not split the matrix: keep it in the +* block so that it is reported below instead of being +* dropped or isolated. + IF( .NOT.( ABS( E( FINISH ) ).LE.TINY ) ) THEN FINISH = FINISH + 1 GO TO 40 END IF @@ -409,6 +412,10 @@ SUBROUTINE CSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, * Scale. * ORGNRM = SLANST( 'M', M, D( START ), E( START ) ) + IF( SISNAN( ORGNRM ) ) THEN + INFO = START*( N+1 ) + FINISH + GO TO 70 + END IF CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), $ M, $ INFO ) diff --git a/SRC/dstedc.f b/SRC/dstedc.f index 02c3719fb..14ad9d0a4 100644 --- a/SRC/dstedc.f +++ b/SRC/dstedc.f @@ -205,10 +205,10 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, DOUBLE PRECISION EPS, ORGNRM, P, TINY * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL DISNAN, LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLANST - EXTERNAL LSAME, ILAENV, DLAMCH, DLANST + EXTERNAL DISNAN, LSAME, ILAENV, DLAMCH, DLANST * .. * .. External Subroutines .. EXTERNAL DGEMM, DLACPY, DLAED0, DLASCL, DLASET, @@ -359,7 +359,10 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, IF( FINISH.LT.N ) THEN TINY = EPS*SQRT( ABS( D( FINISH ) ) )* $ SQRT( ABS( D( FINISH+1 ) ) ) - IF( ABS( E( FINISH ) ).GT.TINY ) THEN +* A NaN in D or E must not split the matrix: keep it in the +* block so that it is reported below instead of being +* dropped or isolated. + IF( .NOT.( ABS( E( FINISH ) ).LE.TINY ) ) THEN FINISH = FINISH + 1 GO TO 20 END IF @@ -377,6 +380,10 @@ SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, * Scale. * ORGNRM = DLANST( 'M', M, D( START ), E( START ) ) + IF( DISNAN( ORGNRM ) ) THEN + INFO = START*( N+1 ) + FINISH + GO TO 50 + END IF CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), $ M, $ INFO ) diff --git a/SRC/sstedc.f b/SRC/sstedc.f index 140a51a88..bd18a3bab 100644 --- a/SRC/sstedc.f +++ b/SRC/sstedc.f @@ -205,10 +205,10 @@ SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, REAL EPS, ORGNRM, P, TINY * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL SISNAN, LSAME INTEGER ILAENV REAL SLAMCH, SLANST, SROUNDUP_LWORK - EXTERNAL ILAENV, LSAME, SLAMCH, SLANST, + EXTERNAL ILAENV, LSAME, SISNAN, SLAMCH, SLANST, $ SROUNDUP_LWORK * .. * .. External Subroutines .. @@ -360,7 +360,10 @@ SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, IF( FINISH.LT.N ) THEN TINY = EPS*SQRT( ABS( D( FINISH ) ) )* $ SQRT( ABS( D( FINISH+1 ) ) ) - IF( ABS( E( FINISH ) ).GT.TINY ) THEN +* A NaN in D or E must not split the matrix: keep it in the +* block so that it is reported below instead of being +* dropped or isolated. + IF( .NOT.( ABS( E( FINISH ) ).LE.TINY ) ) THEN FINISH = FINISH + 1 GO TO 20 END IF @@ -378,6 +381,10 @@ SUBROUTINE SSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK, * Scale. * ORGNRM = SLANST( 'M', M, D( START ), E( START ) ) + IF( SISNAN( ORGNRM ) ) THEN + INFO = START*( N+1 ) + FINISH + GO TO 50 + END IF CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), $ M, $ INFO ) diff --git a/SRC/zstedc.f b/SRC/zstedc.f index 4a5d9fa69..e1f710a27 100644 --- a/SRC/zstedc.f +++ b/SRC/zstedc.f @@ -230,10 +230,10 @@ SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, DOUBLE PRECISION EPS, ORGNRM, P, TINY * .. * .. External Functions .. - LOGICAL LSAME + LOGICAL DISNAN, LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLANST - EXTERNAL LSAME, ILAENV, DLAMCH, DLANST + EXTERNAL DISNAN, LSAME, ILAENV, DLAMCH, DLANST * .. * .. External Subroutines .. EXTERNAL DLASCL, DLASET, DSTEDC, DSTEQR, DSTERF, @@ -394,7 +394,10 @@ SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, IF( FINISH.LT.N ) THEN TINY = EPS*SQRT( ABS( D( FINISH ) ) )* $ SQRT( ABS( D( FINISH+1 ) ) ) - IF( ABS( E( FINISH ) ).GT.TINY ) THEN +* A NaN in D or E must not split the matrix: keep it in the +* block so that it is reported below instead of being +* dropped or isolated. + IF( .NOT.( ABS( E( FINISH ) ).LE.TINY ) ) THEN FINISH = FINISH + 1 GO TO 40 END IF @@ -408,6 +411,10 @@ SUBROUTINE ZSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, RWORK, * Scale. * ORGNRM = DLANST( 'M', M, D( START ), E( START ) ) + IF( DISNAN( ORGNRM ) ) THEN + INFO = START*( N+1 ) + FINISH + GO TO 70 + END IF CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), $ M, $ INFO ) diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 790c7f5f1..64b5f7e5b 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -651,6 +651,9 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER COMPZ + INTEGER JCOMPZ, JNAN, NSMLSZ + REAL RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -1953,11 +1956,62 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* CSTEDC must report a NaN in the matrix through INFO. A NaN +* off-diagonal entry used to split the matrix and was dropped, and +* a NaN diagonal entry was isolated as a 1 by 1 block, both with +* INFO = 0, whenever the block left over was large enough for the +* divide and conquer recursion, which the sizes in the input file +* do not reach; N = LDU is the capacity of the arrays. +* + NSMLSZ = ILAENV( 9, 'CSTEDC', ' ', 0, 0, 0, 0 ) + IF( LDU.GT.NSMLSZ .AND. + $ ILAENV( 10, 'CSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'CSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN + N = LDU + RONE = ONE + RNAN = SQRT( -RONE ) + CALL CSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, RWORK, -1, + $ IWORK, -1, IINFO ) + IF( IINFO.EQ.0 .AND. INT( REAL( WORK( 1 ) ) ).LE.LWORK .AND. + $ INT( RWORK( 1 ) ).LE.LRWORK .AND. IWORK( 1 ).LE.LIWORK ) + $ THEN + DO 340 JNAN = 1, 2 + DO 330 JCOMPZ = 1, 2 + IF( JCOMPZ.EQ.1 ) THEN + COMPZ = 'I' + ELSE + COMPZ = 'V' + END IF + DO 320 J = 1, N + SD( J ) = REAL( J ) + SE( J ) = ONE / REAL( J+1 ) + 320 CONTINUE + IF( JNAN.EQ.1 ) THEN + SE( N / 2 ) = RNAN + ELSE + SD( N ) = RNAN + END IF + CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) + CALL CSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, + $ RWORK, LRWORK, IWORK, LIWORK, IINFO ) + IF( IINFO.EQ.0 ) THEN + WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + 340 CONTINUE + END IF + END IF +* * Summary * CALL SLASUM( 'CST', NOUNIT, NERRS, NTESTT ) RETURN * + 9985 FORMAT( ' CCHKST: CSTEDC( ', A1, ' ) returned INFO=0 for a', + $ ' matrix with a NaN, case ', I1 ) 9999 FORMAT( ' CCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * diff --git a/TESTING/EIG/dchkst.f b/TESTING/EIG/dchkst.f index 012aa95d4..d0d51d455 100644 --- a/TESTING/EIG/dchkst.f +++ b/TESTING/EIG/dchkst.f @@ -634,6 +634,9 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER COMPZ + INTEGER JCOMPZ, JNAN, NSMLSZ + DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -1931,11 +1934,61 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* DSTEDC must report a NaN in the matrix through INFO. A NaN +* off-diagonal entry used to split the matrix and was dropped, and +* a NaN diagonal entry was isolated as a 1 by 1 block, both with +* INFO = 0, whenever the block left over was large enough for the +* divide and conquer recursion, which the sizes in the input file +* do not reach; N = LDU is the capacity of the arrays. +* + NSMLSZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 ) + IF( LDU.GT.NSMLSZ .AND. + $ ILAENV( 10, 'DSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'DSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN + N = LDU + RONE = ONE + RNAN = SQRT( -RONE ) + CALL DSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, IWORK, -1, + $ IINFO ) + IF( IINFO.EQ.0 .AND. INT( WORK( 1 ) ).LE.LWORK .AND. + $ IWORK( 1 ).LE.LIWORK ) THEN + DO 340 JNAN = 1, 2 + DO 330 JCOMPZ = 1, 2 + IF( JCOMPZ.EQ.1 ) THEN + COMPZ = 'I' + ELSE + COMPZ = 'V' + END IF + DO 320 J = 1, N + SD( J ) = DBLE( J ) + SE( J ) = ONE / DBLE( J+1 ) + 320 CONTINUE + IF( JNAN.EQ.1 ) THEN + SE( N / 2 ) = RNAN + ELSE + SD( N ) = RNAN + END IF + CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) + CALL DSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, + $ IWORK, LIWORK, IINFO ) + IF( IINFO.EQ.0 ) THEN + WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + 340 CONTINUE + END IF + END IF +* * Summary * CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT ) RETURN * + 9985 FORMAT( ' DCHKST: DSTEDC( ', A1, ' ) returned INFO=0 for a', + $ ' matrix with a NaN, case ', I1 ) 9999 FORMAT( ' DCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * diff --git a/TESTING/EIG/schkst.f b/TESTING/EIG/schkst.f index c5c7f57c3..f37215611 100644 --- a/TESTING/EIG/schkst.f +++ b/TESTING/EIG/schkst.f @@ -634,6 +634,9 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER COMPZ + INTEGER JCOMPZ, JNAN, NSMLSZ + REAL RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -1931,11 +1934,61 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* SSTEDC must report a NaN in the matrix through INFO. A NaN +* off-diagonal entry used to split the matrix and was dropped, and +* a NaN diagonal entry was isolated as a 1 by 1 block, both with +* INFO = 0, whenever the block left over was large enough for the +* divide and conquer recursion, which the sizes in the input file +* do not reach; N = LDU is the capacity of the arrays. +* + NSMLSZ = ILAENV( 9, 'SSTEDC', ' ', 0, 0, 0, 0 ) + IF( LDU.GT.NSMLSZ .AND. + $ ILAENV( 10, 'SSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'SSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN + N = LDU + RONE = ONE + RNAN = SQRT( -RONE ) + CALL SSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, IWORK, -1, + $ IINFO ) + IF( IINFO.EQ.0 .AND. INT( WORK( 1 ) ).LE.LWORK .AND. + $ IWORK( 1 ).LE.LIWORK ) THEN + DO 340 JNAN = 1, 2 + DO 330 JCOMPZ = 1, 2 + IF( JCOMPZ.EQ.1 ) THEN + COMPZ = 'I' + ELSE + COMPZ = 'V' + END IF + DO 320 J = 1, N + SD( J ) = REAL( J ) + SE( J ) = ONE / REAL( J+1 ) + 320 CONTINUE + IF( JNAN.EQ.1 ) THEN + SE( N / 2 ) = RNAN + ELSE + SD( N ) = RNAN + END IF + CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) + CALL SSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, + $ IWORK, LIWORK, IINFO ) + IF( IINFO.EQ.0 ) THEN + WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + 340 CONTINUE + END IF + END IF +* * Summary * CALL SLASUM( 'SST', NOUNIT, NERRS, NTESTT ) RETURN * + 9985 FORMAT( ' SCHKST: SSTEDC( ', A1, ' ) returned INFO=0 for a', + $ ' matrix with a NaN, case ', I1 ) 9999 FORMAT( ' SCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * diff --git a/TESTING/EIG/zchkst.f b/TESTING/EIG/zchkst.f index 4335e15f0..aae6d1e52 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -651,6 +651,9 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER COMPZ + INTEGER JCOMPZ, JNAN, NSMLSZ + DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -1952,11 +1955,62 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* ZSTEDC must report a NaN in the matrix through INFO. A NaN +* off-diagonal entry used to split the matrix and was dropped, and +* a NaN diagonal entry was isolated as a 1 by 1 block, both with +* INFO = 0, whenever the block left over was large enough for the +* divide and conquer recursion, which the sizes in the input file +* do not reach; N = LDU is the capacity of the arrays. +* + NSMLSZ = ILAENV( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 ) + IF( LDU.GT.NSMLSZ .AND. + $ ILAENV( 10, 'ZSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'ZSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN + N = LDU + RONE = ONE + RNAN = SQRT( -RONE ) + CALL ZSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, RWORK, -1, + $ IWORK, -1, IINFO ) + IF( IINFO.EQ.0 .AND. INT( DBLE( WORK( 1 ) ) ).LE.LWORK .AND. + $ INT( RWORK( 1 ) ).LE.LRWORK .AND. IWORK( 1 ).LE.LIWORK ) + $ THEN + DO 340 JNAN = 1, 2 + DO 330 JCOMPZ = 1, 2 + IF( JCOMPZ.EQ.1 ) THEN + COMPZ = 'I' + ELSE + COMPZ = 'V' + END IF + DO 320 J = 1, N + SD( J ) = DBLE( J ) + SE( J ) = ONE / DBLE( J+1 ) + 320 CONTINUE + IF( JNAN.EQ.1 ) THEN + SE( N / 2 ) = RNAN + ELSE + SD( N ) = RNAN + END IF + CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) + CALL ZSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, + $ RWORK, LRWORK, IWORK, LIWORK, IINFO ) + IF( IINFO.EQ.0 ) THEN + WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + 340 CONTINUE + END IF + END IF +* * Summary * CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT ) RETURN * + 9985 FORMAT( ' ZCHKST: ZSTEDC( ', A1, ' ) returned INFO=0 for a', + $ ' matrix with a NaN, case ', I1 ) 9999 FORMAT( ' ZCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * From 8cca9b516d2d2ca8a0a698befe0f56e350ffed59 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Thu, 10 Sep 2026 17:51:27 -0700 Subject: [PATCH 2/2] TESTING: Bound the STEDC NaN regression storage Use local arrays sized above SMLSIZ instead of treating LDU as the capacity of the caller's vectors and matrix columns. Preserve both NaN positions and both eigenvector options. Validation: 8 sep/se2 driver runs, 32 AddressSanitizer capacity cases, and 4 runs against the parent kernels that retain the expected failures. --- TESTING/EIG/cchkst.f | 29 +++++++++++++++++------------ TESTING/EIG/dchkst.f | 28 ++++++++++++++++------------ TESTING/EIG/schkst.f | 28 ++++++++++++++++------------ TESTING/EIG/zchkst.f | 29 +++++++++++++++++------------ 4 files changed, 66 insertions(+), 48 deletions(-) diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 64b5f7e5b..2f8c2415d 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -659,6 +659,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ KMAGN( MAXTYP ), KMODE( MAXTYP ), $ KTYPE( MAXTYP ) REAL DUMMA( 1 ) + REAL, ALLOCATABLE :: DREG( : ), EREG( : ) + COMPLEX, ALLOCATABLE :: ZREG( :, : ) * .. * .. External Functions .. INTEGER ILAENV @@ -1962,16 +1964,17 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * a NaN diagonal entry was isolated as a 1 by 1 block, both with * INFO = 0, whenever the block left over was large enough for the * divide and conquer recursion, which the sizes in the input file -* do not reach; N = LDU is the capacity of the arrays. +* do not reach. Use local arrays because LDU is a row stride, +* not the capacity of the caller's eigenvalue arrays or Z columns. * NSMLSZ = ILAENV( 9, 'CSTEDC', ' ', 0, 0, 0, 0 ) - IF( LDU.GT.NSMLSZ .AND. - $ ILAENV( 10, 'CSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + IF( ILAENV( 10, 'CSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. $ ILAENV( 11, 'CSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN - N = LDU + N = MAX( 2, NSMLSZ+1 ) + ALLOCATE( DREG( N ), EREG( N ), ZREG( N, N ) ) RONE = ONE RNAN = SQRT( -RONE ) - CALL CSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, RWORK, -1, + CALL CSTEDC( 'V', N, DREG, EREG, ZREG, N, WORK, -1, RWORK, -1, $ IWORK, -1, IINFO ) IF( IINFO.EQ.0 .AND. INT( REAL( WORK( 1 ) ) ).LE.LWORK .AND. $ INT( RWORK( 1 ) ).LE.LRWORK .AND. IWORK( 1 ).LE.LIWORK ) @@ -1984,17 +1987,18 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, COMPZ = 'V' END IF DO 320 J = 1, N - SD( J ) = REAL( J ) - SE( J ) = ONE / REAL( J+1 ) + DREG( J ) = REAL( J ) + EREG( J ) = ONE / REAL( J+1 ) 320 CONTINUE IF( JNAN.EQ.1 ) THEN - SE( N / 2 ) = RNAN + EREG( N / 2 ) = RNAN ELSE - SD( N ) = RNAN + DREG( N ) = RNAN END IF - CALL CLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) - CALL CSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, - $ RWORK, LRWORK, IWORK, LIWORK, IINFO ) + CALL CLASET( 'Full', N, N, CZERO, CONE, ZREG, N ) + CALL CSTEDC( COMPZ, N, DREG, EREG, ZREG, N, + $ WORK, LWORK, RWORK, LRWORK, IWORK, + $ LIWORK, IINFO ) IF( IINFO.EQ.0 ) THEN WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN NERRS = NERRS + 1 @@ -2003,6 +2007,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 330 CONTINUE 340 CONTINUE END IF + DEALLOCATE( DREG, EREG, ZREG ) END IF * * Summary diff --git a/TESTING/EIG/dchkst.f b/TESTING/EIG/dchkst.f index d0d51d455..d22355112 100644 --- a/TESTING/EIG/dchkst.f +++ b/TESTING/EIG/dchkst.f @@ -642,6 +642,8 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ KMAGN( MAXTYP ), KMODE( MAXTYP ), $ KTYPE( MAXTYP ) DOUBLE PRECISION DUMMA( 1 ) + DOUBLE PRECISION, ALLOCATABLE :: DREG( : ), EREG( : ) + DOUBLE PRECISION, ALLOCATABLE :: ZREG( :, : ) * .. * .. External Functions .. INTEGER ILAENV @@ -1940,16 +1942,17 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * a NaN diagonal entry was isolated as a 1 by 1 block, both with * INFO = 0, whenever the block left over was large enough for the * divide and conquer recursion, which the sizes in the input file -* do not reach; N = LDU is the capacity of the arrays. +* do not reach. Use local arrays because LDU is a row stride, +* not the capacity of the caller's eigenvalue arrays or Z columns. * NSMLSZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 ) - IF( LDU.GT.NSMLSZ .AND. - $ ILAENV( 10, 'DSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + IF( ILAENV( 10, 'DSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. $ ILAENV( 11, 'DSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN - N = LDU + N = MAX( 2, NSMLSZ+1 ) + ALLOCATE( DREG( N ), EREG( N ), ZREG( N, N ) ) RONE = ONE RNAN = SQRT( -RONE ) - CALL DSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, IWORK, -1, + CALL DSTEDC( 'V', N, DREG, EREG, ZREG, N, WORK, -1, IWORK, -1, $ IINFO ) IF( IINFO.EQ.0 .AND. INT( WORK( 1 ) ).LE.LWORK .AND. $ IWORK( 1 ).LE.LIWORK ) THEN @@ -1961,17 +1964,17 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, COMPZ = 'V' END IF DO 320 J = 1, N - SD( J ) = DBLE( J ) - SE( J ) = ONE / DBLE( J+1 ) + DREG( J ) = DBLE( J ) + EREG( J ) = ONE / DBLE( J+1 ) 320 CONTINUE IF( JNAN.EQ.1 ) THEN - SE( N / 2 ) = RNAN + EREG( N / 2 ) = RNAN ELSE - SD( N ) = RNAN + DREG( N ) = RNAN END IF - CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) - CALL DSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, - $ IWORK, LIWORK, IINFO ) + CALL DLASET( 'Full', N, N, ZERO, ONE, ZREG, N ) + CALL DSTEDC( COMPZ, N, DREG, EREG, ZREG, N, + $ WORK, LWORK, IWORK, LIWORK, IINFO ) IF( IINFO.EQ.0 ) THEN WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN NERRS = NERRS + 1 @@ -1980,6 +1983,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 330 CONTINUE 340 CONTINUE END IF + DEALLOCATE( DREG, EREG, ZREG ) END IF * * Summary diff --git a/TESTING/EIG/schkst.f b/TESTING/EIG/schkst.f index f37215611..7f7ced9f8 100644 --- a/TESTING/EIG/schkst.f +++ b/TESTING/EIG/schkst.f @@ -642,6 +642,8 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ KMAGN( MAXTYP ), KMODE( MAXTYP ), $ KTYPE( MAXTYP ) REAL DUMMA( 1 ) + REAL, ALLOCATABLE :: DREG( : ), EREG( : ) + REAL, ALLOCATABLE :: ZREG( :, : ) * .. * .. External Functions .. INTEGER ILAENV @@ -1940,16 +1942,17 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * a NaN diagonal entry was isolated as a 1 by 1 block, both with * INFO = 0, whenever the block left over was large enough for the * divide and conquer recursion, which the sizes in the input file -* do not reach; N = LDU is the capacity of the arrays. +* do not reach. Use local arrays because LDU is a row stride, +* not the capacity of the caller's eigenvalue arrays or Z columns. * NSMLSZ = ILAENV( 9, 'SSTEDC', ' ', 0, 0, 0, 0 ) - IF( LDU.GT.NSMLSZ .AND. - $ ILAENV( 10, 'SSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + IF( ILAENV( 10, 'SSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. $ ILAENV( 11, 'SSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN - N = LDU + N = MAX( 2, NSMLSZ+1 ) + ALLOCATE( DREG( N ), EREG( N ), ZREG( N, N ) ) RONE = ONE RNAN = SQRT( -RONE ) - CALL SSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, IWORK, -1, + CALL SSTEDC( 'V', N, DREG, EREG, ZREG, N, WORK, -1, IWORK, -1, $ IINFO ) IF( IINFO.EQ.0 .AND. INT( WORK( 1 ) ).LE.LWORK .AND. $ IWORK( 1 ).LE.LIWORK ) THEN @@ -1961,17 +1964,17 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, COMPZ = 'V' END IF DO 320 J = 1, N - SD( J ) = REAL( J ) - SE( J ) = ONE / REAL( J+1 ) + DREG( J ) = REAL( J ) + EREG( J ) = ONE / REAL( J+1 ) 320 CONTINUE IF( JNAN.EQ.1 ) THEN - SE( N / 2 ) = RNAN + EREG( N / 2 ) = RNAN ELSE - SD( N ) = RNAN + DREG( N ) = RNAN END IF - CALL SLASET( 'Full', N, N, ZERO, ONE, Z, LDU ) - CALL SSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, - $ IWORK, LIWORK, IINFO ) + CALL SLASET( 'Full', N, N, ZERO, ONE, ZREG, N ) + CALL SSTEDC( COMPZ, N, DREG, EREG, ZREG, N, + $ WORK, LWORK, IWORK, LIWORK, IINFO ) IF( IINFO.EQ.0 ) THEN WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN NERRS = NERRS + 1 @@ -1980,6 +1983,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 330 CONTINUE 340 CONTINUE END IF + DEALLOCATE( DREG, EREG, ZREG ) END IF * * Summary diff --git a/TESTING/EIG/zchkst.f b/TESTING/EIG/zchkst.f index aae6d1e52..8b9448455 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -659,6 +659,8 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ KMAGN( MAXTYP ), KMODE( MAXTYP ), $ KTYPE( MAXTYP ) DOUBLE PRECISION DUMMA( 1 ) + DOUBLE PRECISION, ALLOCATABLE :: DREG( : ), EREG( : ) + COMPLEX*16, ALLOCATABLE :: ZREG( :, : ) * .. * .. External Functions .. INTEGER ILAENV @@ -1961,16 +1963,17 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * a NaN diagonal entry was isolated as a 1 by 1 block, both with * INFO = 0, whenever the block left over was large enough for the * divide and conquer recursion, which the sizes in the input file -* do not reach; N = LDU is the capacity of the arrays. +* do not reach. Use local arrays because LDU is a row stride, +* not the capacity of the caller's eigenvalue arrays or Z columns. * NSMLSZ = ILAENV( 9, 'ZSTEDC', ' ', 0, 0, 0, 0 ) - IF( LDU.GT.NSMLSZ .AND. - $ ILAENV( 10, 'ZSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. + IF( ILAENV( 10, 'ZSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 .AND. $ ILAENV( 11, 'ZSTEDC', 'V', 1, 0, 0, 0 ).EQ.1 ) THEN - N = LDU + N = MAX( 2, NSMLSZ+1 ) + ALLOCATE( DREG( N ), EREG( N ), ZREG( N, N ) ) RONE = ONE RNAN = SQRT( -RONE ) - CALL ZSTEDC( 'V', N, SD, SE, Z, LDU, WORK, -1, RWORK, -1, + CALL ZSTEDC( 'V', N, DREG, EREG, ZREG, N, WORK, -1, RWORK, -1, $ IWORK, -1, IINFO ) IF( IINFO.EQ.0 .AND. INT( DBLE( WORK( 1 ) ) ).LE.LWORK .AND. $ INT( RWORK( 1 ) ).LE.LRWORK .AND. IWORK( 1 ).LE.LIWORK ) @@ -1983,17 +1986,18 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, COMPZ = 'V' END IF DO 320 J = 1, N - SD( J ) = DBLE( J ) - SE( J ) = ONE / DBLE( J+1 ) + DREG( J ) = DBLE( J ) + EREG( J ) = ONE / DBLE( J+1 ) 320 CONTINUE IF( JNAN.EQ.1 ) THEN - SE( N / 2 ) = RNAN + EREG( N / 2 ) = RNAN ELSE - SD( N ) = RNAN + DREG( N ) = RNAN END IF - CALL ZLASET( 'Full', N, N, CZERO, CONE, Z, LDU ) - CALL ZSTEDC( COMPZ, N, SD, SE, Z, LDU, WORK, LWORK, - $ RWORK, LRWORK, IWORK, LIWORK, IINFO ) + CALL ZLASET( 'Full', N, N, CZERO, CONE, ZREG, N ) + CALL ZSTEDC( COMPZ, N, DREG, EREG, ZREG, N, + $ WORK, LWORK, RWORK, LRWORK, IWORK, + $ LIWORK, IINFO ) IF( IINFO.EQ.0 ) THEN WRITE( NOUNIT, FMT = 9985 )COMPZ, JNAN NERRS = NERRS + 1 @@ -2002,6 +2006,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 330 CONTINUE 340 CONTINUE END IF + DEALLOCATE( DREG, EREG, ZREG ) END IF * * Summary