From 6f96b00e7d1a4ca0ba48d1b5aab9e27a5485db6f Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Thu, 10 Sep 2026 01:11:18 -0700 Subject: [PATCH 1/2] Propagate a NaN through xLAE2 and xLAEV2 instead of returning finite eigenvalues xLAE2 and xLAEV2 compute the eigenvalues of the symmetric 2 by 2 matrix [a b; b c] from sm = a + c, adf = |a - c| and ab = |2b|: IF( ADF.GT.AB ) THEN RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) ELSE RT = AB*SQRT( TWO ) END IF The last branch is meant for adf = ab, but a NaN diagonal entry makes adf a NaN, both comparisons false, and rt = |2b| sqrt(2), finite. The sign test on sm that follows is false for a NaN as well, so the routines returned rt1 = -rt2 = |b| sqrt(2) with no trace of the NaN. Every 2 by 2 tridiagonal eigenproblem goes through this code: xSTEQR, xSTERF, xSTEV, xSTEVD, xSTEDC and xSTEMR returned INFO = 0 and finite eigenvalues for a 2 by 2 matrix with a NaN on the diagonal, while a NaN off the diagonal was propagated. Take the last branch only for adf = ab and let a NaN in either operand propagate, rt = adf + ab, so that the eigenvalues, and in xLAEV2 the eigenvector, come out as NaN. Finite arguments are not affected: the branch they take is unchanged and so is the value it computes. xCHKST gets the case as a regression test: it calls xSTEQR, xSTERF and xSTEMR on the 2 by 2 matrix [1 1; 1 2] with a NaN in place of either diagonal entry and reports a finite eigenvalue as a failure. On the parent commit all six calls return INFO = 0 and finite eigenvalues in every precision; here every eigenvalue is a NaN. Over the 1 by 1 and 2 by 2 NaN cases of a sweep of every tridiagonal solver the parent returned an all-finite result in 62 of 111 cases and this branch in 16, all of them calls that return no eigenvalue at all (RANGE = 'V' or 'I' with M = 0) or bisection paths, which do the same for larger n; every finite case is bit-identical. The full LAPACK test suite passes: 0 numerical errors, 0 other errors, 120 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/dlae2.f | 8 ++++++- SRC/dlaev2.f | 8 ++++++- SRC/slae2.f | 8 ++++++- SRC/slaev2.f | 8 ++++++- TESTING/EIG/cchkst.f | 57 +++++++++++++++++++++++++++++++++++++++++++- TESTING/EIG/dchkst.f | 57 +++++++++++++++++++++++++++++++++++++++++++- TESTING/EIG/schkst.f | 57 +++++++++++++++++++++++++++++++++++++++++++- TESTING/EIG/zchkst.f | 57 +++++++++++++++++++++++++++++++++++++++++++- 8 files changed, 252 insertions(+), 8 deletions(-) diff --git a/SRC/dlae2.f b/SRC/dlae2.f index 388b82786..11d8a2dbf 100644 --- a/SRC/dlae2.f +++ b/SRC/dlae2.f @@ -145,11 +145,17 @@ SUBROUTINE DLAE2( A, B, C, RT1, RT2 ) RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) - ELSE + ELSE IF( ADF.EQ.AB ) THEN * * Includes case AB=ADF=0 * RT = AB*SQRT( TWO ) + ELSE +* +* ADF or AB is a NaN; propagate it instead of returning +* finite eigenvalues for a matrix that contains a NaN. +* + RT = ADF + AB END IF IF( SM.LT.ZERO ) THEN RT1 = HALF*( SM-RT ) diff --git a/SRC/dlaev2.f b/SRC/dlaev2.f index b808e708d..1a73fc858 100644 --- a/SRC/dlaev2.f +++ b/SRC/dlaev2.f @@ -165,11 +165,17 @@ SUBROUTINE DLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) - ELSE + ELSE IF( ADF.EQ.AB ) THEN * * Includes case AB=ADF=0 * RT = AB*SQRT( TWO ) + ELSE +* +* ADF or AB is a NaN; propagate it instead of returning +* finite eigenvalues for a matrix that contains a NaN. +* + RT = ADF + AB END IF IF( SM.LT.ZERO ) THEN RT1 = HALF*( SM-RT ) diff --git a/SRC/slae2.f b/SRC/slae2.f index 2e3b9ac67..6cb115bcd 100644 --- a/SRC/slae2.f +++ b/SRC/slae2.f @@ -145,11 +145,17 @@ SUBROUTINE SLAE2( A, B, C, RT1, RT2 ) RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) - ELSE + ELSE IF( ADF.EQ.AB ) THEN * * Includes case AB=ADF=0 * RT = AB*SQRT( TWO ) + ELSE +* +* ADF or AB is a NaN; propagate it instead of returning +* finite eigenvalues for a matrix that contains a NaN. +* + RT = ADF + AB END IF IF( SM.LT.ZERO ) THEN RT1 = HALF*( SM-RT ) diff --git a/SRC/slaev2.f b/SRC/slaev2.f index ed4030198..64ad6d37f 100644 --- a/SRC/slaev2.f +++ b/SRC/slaev2.f @@ -165,11 +165,17 @@ SUBROUTINE SLAEV2( A, B, C, RT1, RT2, CS1, SN1 ) RT = ADF*SQRT( ONE+( AB / ADF )**2 ) ELSE IF( ADF.LT.AB ) THEN RT = AB*SQRT( ONE+( ADF / AB )**2 ) - ELSE + ELSE IF( ADF.EQ.AB ) THEN * * Includes case AB=ADF=0 * RT = AB*SQRT( TWO ) + ELSE +* +* ADF or AB is a NaN; propagate it instead of returning +* finite eigenvalues for a matrix that contains a NaN. +* + RT = ADF + AB END IF IF( SM.LT.ZERO ) THEN RT1 = HALF*( SM-RT ) diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 790c7f5f1..62cd49fc3 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*6 RNAME + INTEGER JNAN, JROUT + REAL RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -658,9 +661,10 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, REAL DUMMA( 1 ) * .. * .. External Functions .. + LOGICAL SISNAN INTEGER ILAENV REAL SLAMCH, SLARND, SSXT1 - EXTERNAL ILAENV, SLAMCH, SLARND, SSXT1 + EXTERNAL SISNAN, ILAENV, SLAMCH, SLARND, SSXT1 * .. * .. External Subroutines .. EXTERNAL CCOPY, CHET21, CHETRD, CHPT21, CHPTRD, CLACPY, @@ -1953,11 +1957,62 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* A 2 by 2 matrix with a NaN on its diagonal must not get finite +* eigenvalues: SLAE2 and SLAEV2 fell through the comparisons of +* |a-c| with |2b| and of a+c with zero, which are all false for a +* NaN, and returned +/- |b| sqrt(2). +* + IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. + $ ILAENV( 10, 'CSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'CSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 ) THEN + N = 2 + RONE = ONE + RNAN = SQRT( -RONE ) + DO 380 JNAN = 1, 2 + DO 370 JROUT = 1, 3 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + IF( JROUT.EQ.1 ) THEN + RNAME = 'CSTEQR' + CALL CSTEQR( 'N', N, SD, SE, Z, LDU, RWORK, IINFO ) + ELSE IF( JROUT.EQ.2 ) THEN + RNAME = 'SSTERF' + CALL SSTERF( N, SD, SE, IINFO ) + ELSE + RNAME = 'CSTEMR' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL CSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK, LRWORK, + $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.EQ.0 .AND. M.EQ.N ) + $ CALL SCOPY( N, WR, 1, SD, 1 ) + END IF + IF( IINFO.EQ.0 .AND. .NOT.( SISNAN( SD( 1 ) ) .AND. + $ SISNAN( SD( 2 ) ) ) ) THEN + WRITE( NOUNIT, FMT = 9983 )RNAME, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 370 CONTINUE + 380 CONTINUE + END IF +* * Summary * CALL SLASUM( 'CST', NOUNIT, NERRS, NTESTT ) RETURN * + 9983 FORMAT( ' CCHKST: ', A6, ' returned a finite eigenvalue for a', + $ ' 2 by 2 matrix with a NaN at D(', 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..e32d6952f 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*6 RNAME + INTEGER JNAN, JROUT + DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -641,9 +644,10 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, DOUBLE PRECISION DUMMA( 1 ) * .. * .. External Functions .. + LOGICAL DISNAN INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLARND, DSXT1 - EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1 + EXTERNAL DISNAN, ILAENV, DLAMCH, DLARND, DSXT1 * .. * .. External Subroutines .. EXTERNAL DCOPY, DLACPY, DLASET, DLASUM, DLATMR, DLATMS, @@ -1931,11 +1935,62 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* A 2 by 2 matrix with a NaN on its diagonal must not get finite +* eigenvalues: DLAE2 and DLAEV2 fell through the comparisons of +* |a-c| with |2b| and of a+c with zero, which are all false for a +* NaN, and returned +/- |b| sqrt(2). +* + IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. + $ ILAENV( 10, 'DSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'DSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 ) THEN + N = 2 + RONE = ONE + RNAN = SQRT( -RONE ) + DO 380 JNAN = 1, 2 + DO 370 JROUT = 1, 3 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + IF( JROUT.EQ.1 ) THEN + RNAME = 'DSTEQR' + CALL DSTEQR( 'N', N, SD, SE, Z, LDU, WORK, IINFO ) + ELSE IF( JROUT.EQ.2 ) THEN + RNAME = 'DSTERF' + CALL DSTERF( N, SD, SE, IINFO ) + ELSE + RNAME = 'DSTEMR' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL DSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ WORK, LWORK, + $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.EQ.0 .AND. M.EQ.N ) + $ CALL DCOPY( N, WR, 1, SD, 1 ) + END IF + IF( IINFO.EQ.0 .AND. .NOT.( DISNAN( SD( 1 ) ) .AND. + $ DISNAN( SD( 2 ) ) ) ) THEN + WRITE( NOUNIT, FMT = 9983 )RNAME, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 370 CONTINUE + 380 CONTINUE + END IF +* * Summary * CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT ) RETURN * + 9983 FORMAT( ' DCHKST: ', A6, ' returned a finite eigenvalue for a', + $ ' 2 by 2 matrix with a NaN at D(', 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..48ef35197 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*6 RNAME + INTEGER JNAN, JROUT + REAL RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -641,9 +644,10 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, REAL DUMMA( 1 ) * .. * .. External Functions .. + LOGICAL SISNAN INTEGER ILAENV REAL SLAMCH, SLARND, SSXT1 - EXTERNAL ILAENV, SLAMCH, SLARND, SSXT1 + EXTERNAL SISNAN, ILAENV, SLAMCH, SLARND, SSXT1 * .. * .. External Subroutines .. EXTERNAL SCOPY, SLACPY, SLASET, SLASUM, SLATMR, SLATMS, @@ -1931,11 +1935,62 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* A 2 by 2 matrix with a NaN on its diagonal must not get finite +* eigenvalues: SLAE2 and SLAEV2 fell through the comparisons of +* |a-c| with |2b| and of a+c with zero, which are all false for a +* NaN, and returned +/- |b| sqrt(2). +* + IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. + $ ILAENV( 10, 'SSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'SSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 ) THEN + N = 2 + RONE = ONE + RNAN = SQRT( -RONE ) + DO 380 JNAN = 1, 2 + DO 370 JROUT = 1, 3 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + IF( JROUT.EQ.1 ) THEN + RNAME = 'SSTEQR' + CALL SSTEQR( 'N', N, SD, SE, Z, LDU, WORK, IINFO ) + ELSE IF( JROUT.EQ.2 ) THEN + RNAME = 'SSTERF' + CALL SSTERF( N, SD, SE, IINFO ) + ELSE + RNAME = 'SSTEMR' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL SSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ WORK, LWORK, + $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.EQ.0 .AND. M.EQ.N ) + $ CALL SCOPY( N, WR, 1, SD, 1 ) + END IF + IF( IINFO.EQ.0 .AND. .NOT.( SISNAN( SD( 1 ) ) .AND. + $ SISNAN( SD( 2 ) ) ) ) THEN + WRITE( NOUNIT, FMT = 9983 )RNAME, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 370 CONTINUE + 380 CONTINUE + END IF +* * Summary * CALL SLASUM( 'SST', NOUNIT, NERRS, NTESTT ) RETURN * + 9983 FORMAT( ' SCHKST: ', A6, ' returned a finite eigenvalue for a', + $ ' 2 by 2 matrix with a NaN at D(', 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..fc81a3dbd 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*6 RNAME + INTEGER JNAN, JROUT + DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -658,9 +661,10 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, DOUBLE PRECISION DUMMA( 1 ) * .. * .. External Functions .. + LOGICAL DISNAN INTEGER ILAENV DOUBLE PRECISION DLAMCH, DLARND, DSXT1 - EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1 + EXTERNAL DISNAN, ILAENV, DLAMCH, DLARND, DSXT1 * .. * .. External Subroutines .. EXTERNAL DCOPY, DLASUM, DSTEBZ, DSTECH, DSTERF, XERBLA, @@ -1952,11 +1956,62 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 300 CONTINUE 310 CONTINUE * +* +* A 2 by 2 matrix with a NaN on its diagonal must not get finite +* eigenvalues: DLAE2 and DLAEV2 fell through the comparisons of +* |a-c| with |2b| and of a+c with zero, which are all false for a +* NaN, and returned +/- |b| sqrt(2). +* + IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. + $ ILAENV( 10, 'ZSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'ZSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 ) THEN + N = 2 + RONE = ONE + RNAN = SQRT( -RONE ) + DO 380 JNAN = 1, 2 + DO 370 JROUT = 1, 3 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + IF( JROUT.EQ.1 ) THEN + RNAME = 'ZSTEQR' + CALL ZSTEQR( 'N', N, SD, SE, Z, LDU, RWORK, IINFO ) + ELSE IF( JROUT.EQ.2 ) THEN + RNAME = 'DSTERF' + CALL DSTERF( N, SD, SE, IINFO ) + ELSE + RNAME = 'ZSTEMR' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL ZSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, + $ RWORK, LRWORK, + $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.EQ.0 .AND. M.EQ.N ) + $ CALL DCOPY( N, WR, 1, SD, 1 ) + END IF + IF( IINFO.EQ.0 .AND. .NOT.( DISNAN( SD( 1 ) ) .AND. + $ DISNAN( SD( 2 ) ) ) ) THEN + WRITE( NOUNIT, FMT = 9983 )RNAME, JNAN + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 370 CONTINUE + 380 CONTINUE + END IF +* * Summary * CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT ) RETURN * + 9983 FORMAT( ' ZCHKST: ', A6, ' returned a finite eigenvalue for a', + $ ' 2 by 2 matrix with a NaN at D(', I1, ')' ) 9999 FORMAT( ' ZCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * From 349f09fc532d07dcbf997e3e6afa7eaa1f3a28a6 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: Cover LAEV2 NaN propagation with eigenvectors Add STEQR(I) and STEMR(V,A) cases for both diagonal NaN positions in all four precisions. Include the job option in diagnostics and check the real workspace capacity for the complex drivers. Validation: 8 sep/se2 driver runs and 8 AddressSanitizer cases passed. All four precisions detect reverting only SLAEV2/DLAEV2, and all four also retain the expected failures against the parent kernels. --- TESTING/EIG/cchkst.f | 25 +++++++++++++++---------- TESTING/EIG/dchkst.f | 23 ++++++++++++++--------- TESTING/EIG/schkst.f | 23 ++++++++++++++--------- TESTING/EIG/zchkst.f | 25 +++++++++++++++---------- 4 files changed, 58 insertions(+), 38 deletions(-) diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 62cd49fc3..01329947f 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -651,7 +651,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. - CHARACTER*6 RNAME + CHARACTER OPT + CHARACTER*9 RNAME INTEGER JNAN, JROUT REAL RNAN, RONE * .. Local Arrays .. @@ -1961,35 +1962,39 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * A 2 by 2 matrix with a NaN on its diagonal must not get finite * eigenvalues: SLAE2 and SLAEV2 fell through the comparisons of * |a-c| with |2b| and of a+c with zero, which are all false for a -* NaN, and returned +/- |b| sqrt(2). +* NaN, and returned +/- |b| sqrt(2). Eigenvalue-only calls test +* xLAE2; eigenvector-producing calls also test xLAEV2. * - IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. + IF( NMAX.GE.2 .AND. LRWORK.GE.36 .AND. LIWORK.GE.24 .AND. $ ILAENV( 10, 'CSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. $ ILAENV( 11, 'CSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 ) THEN N = 2 RONE = ONE RNAN = SQRT( -RONE ) DO 380 JNAN = 1, 2 - DO 370 JROUT = 1, 3 + DO 370 JROUT = 1, 5 SD( 1 ) = ONE SD( 2 ) = ONE + ONE SE( 1 ) = ONE SE( 2 ) = ZERO SD( JNAN ) = RNAN - IF( JROUT.EQ.1 ) THEN - RNAME = 'CSTEQR' - CALL CSTEQR( 'N', N, SD, SE, Z, LDU, RWORK, IINFO ) + OPT = 'N' + IF( JROUT.EQ.4 ) OPT = 'I' + IF( JROUT.EQ.5 ) OPT = 'V' + IF( JROUT.EQ.1 .OR. JROUT.EQ.4 ) THEN + RNAME = 'CSTEQR('//OPT//')' + CALL CSTEQR( OPT, N, SD, SE, Z, LDU, RWORK, IINFO ) ELSE IF( JROUT.EQ.2 ) THEN RNAME = 'SSTERF' CALL SSTERF( N, SD, SE, IINFO ) ELSE - RNAME = 'CSTEMR' + RNAME = 'CSTEMR('//OPT//')' VL = ZERO VU = ZERO IL = 0 IU = 0 TRYRAC = .TRUE. - CALL CSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + CALL CSTEMR( OPT, 'A', N, SD, SE, VL, VU, IL, IU, M, $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, $ RWORK, LRWORK, $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) @@ -2011,7 +2016,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, CALL SLASUM( 'CST', NOUNIT, NERRS, NTESTT ) RETURN * - 9983 FORMAT( ' CCHKST: ', A6, ' returned a finite eigenvalue for a', + 9983 FORMAT( ' CCHKST: ', A9, ' returned a finite eigenvalue for a', $ ' 2 by 2 matrix with a NaN at D(', 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 e32d6952f..a8b3433b7 100644 --- a/TESTING/EIG/dchkst.f +++ b/TESTING/EIG/dchkst.f @@ -634,7 +634,8 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. - CHARACTER*6 RNAME + CHARACTER OPT + CHARACTER*9 RNAME INTEGER JNAN, JROUT DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. @@ -1939,7 +1940,8 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * A 2 by 2 matrix with a NaN on its diagonal must not get finite * eigenvalues: DLAE2 and DLAEV2 fell through the comparisons of * |a-c| with |2b| and of a+c with zero, which are all false for a -* NaN, and returned +/- |b| sqrt(2). +* NaN, and returned +/- |b| sqrt(2). Eigenvalue-only calls test +* xLAE2; eigenvector-producing calls also test xLAEV2. * IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. $ ILAENV( 10, 'DSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. @@ -1948,26 +1950,29 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, RONE = ONE RNAN = SQRT( -RONE ) DO 380 JNAN = 1, 2 - DO 370 JROUT = 1, 3 + DO 370 JROUT = 1, 5 SD( 1 ) = ONE SD( 2 ) = ONE + ONE SE( 1 ) = ONE SE( 2 ) = ZERO SD( JNAN ) = RNAN - IF( JROUT.EQ.1 ) THEN - RNAME = 'DSTEQR' - CALL DSTEQR( 'N', N, SD, SE, Z, LDU, WORK, IINFO ) + OPT = 'N' + IF( JROUT.EQ.4 ) OPT = 'I' + IF( JROUT.EQ.5 ) OPT = 'V' + IF( JROUT.EQ.1 .OR. JROUT.EQ.4 ) THEN + RNAME = 'DSTEQR('//OPT//')' + CALL DSTEQR( OPT, N, SD, SE, Z, LDU, WORK, IINFO ) ELSE IF( JROUT.EQ.2 ) THEN RNAME = 'DSTERF' CALL DSTERF( N, SD, SE, IINFO ) ELSE - RNAME = 'DSTEMR' + RNAME = 'DSTEMR('//OPT//')' VL = ZERO VU = ZERO IL = 0 IU = 0 TRYRAC = .TRUE. - CALL DSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + CALL DSTEMR( OPT, 'A', N, SD, SE, VL, VU, IL, IU, M, $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, $ WORK, LWORK, $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) @@ -1989,7 +1994,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT ) RETURN * - 9983 FORMAT( ' DCHKST: ', A6, ' returned a finite eigenvalue for a', + 9983 FORMAT( ' DCHKST: ', A9, ' returned a finite eigenvalue for a', $ ' 2 by 2 matrix with a NaN at D(', 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 48ef35197..bdeb6b8eb 100644 --- a/TESTING/EIG/schkst.f +++ b/TESTING/EIG/schkst.f @@ -634,7 +634,8 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. - CHARACTER*6 RNAME + CHARACTER OPT + CHARACTER*9 RNAME INTEGER JNAN, JROUT REAL RNAN, RONE * .. Local Arrays .. @@ -1939,7 +1940,8 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * A 2 by 2 matrix with a NaN on its diagonal must not get finite * eigenvalues: SLAE2 and SLAEV2 fell through the comparisons of * |a-c| with |2b| and of a+c with zero, which are all false for a -* NaN, and returned +/- |b| sqrt(2). +* NaN, and returned +/- |b| sqrt(2). Eigenvalue-only calls test +* xLAE2; eigenvector-producing calls also test xLAEV2. * IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. $ ILAENV( 10, 'SSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. @@ -1948,26 +1950,29 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, RONE = ONE RNAN = SQRT( -RONE ) DO 380 JNAN = 1, 2 - DO 370 JROUT = 1, 3 + DO 370 JROUT = 1, 5 SD( 1 ) = ONE SD( 2 ) = ONE + ONE SE( 1 ) = ONE SE( 2 ) = ZERO SD( JNAN ) = RNAN - IF( JROUT.EQ.1 ) THEN - RNAME = 'SSTEQR' - CALL SSTEQR( 'N', N, SD, SE, Z, LDU, WORK, IINFO ) + OPT = 'N' + IF( JROUT.EQ.4 ) OPT = 'I' + IF( JROUT.EQ.5 ) OPT = 'V' + IF( JROUT.EQ.1 .OR. JROUT.EQ.4 ) THEN + RNAME = 'SSTEQR('//OPT//')' + CALL SSTEQR( OPT, N, SD, SE, Z, LDU, WORK, IINFO ) ELSE IF( JROUT.EQ.2 ) THEN RNAME = 'SSTERF' CALL SSTERF( N, SD, SE, IINFO ) ELSE - RNAME = 'SSTEMR' + RNAME = 'SSTEMR('//OPT//')' VL = ZERO VU = ZERO IL = 0 IU = 0 TRYRAC = .TRUE. - CALL SSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + CALL SSTEMR( OPT, 'A', N, SD, SE, VL, VU, IL, IU, M, $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, $ WORK, LWORK, $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) @@ -1989,7 +1994,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, CALL SLASUM( 'SST', NOUNIT, NERRS, NTESTT ) RETURN * - 9983 FORMAT( ' SCHKST: ', A6, ' returned a finite eigenvalue for a', + 9983 FORMAT( ' SCHKST: ', A9, ' returned a finite eigenvalue for a', $ ' 2 by 2 matrix with a NaN at D(', 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 fc81a3dbd..6b7674fef 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -651,7 +651,8 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. - CHARACTER*6 RNAME + CHARACTER OPT + CHARACTER*9 RNAME INTEGER JNAN, JROUT DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. @@ -1960,35 +1961,39 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * A 2 by 2 matrix with a NaN on its diagonal must not get finite * eigenvalues: DLAE2 and DLAEV2 fell through the comparisons of * |a-c| with |2b| and of a+c with zero, which are all false for a -* NaN, and returned +/- |b| sqrt(2). +* NaN, and returned +/- |b| sqrt(2). Eigenvalue-only calls test +* xLAE2; eigenvector-producing calls also test xLAEV2. * - IF( NMAX.GE.2 .AND. LWORK.GE.36 .AND. LIWORK.GE.24 .AND. + IF( NMAX.GE.2 .AND. LRWORK.GE.36 .AND. LIWORK.GE.24 .AND. $ ILAENV( 10, 'ZSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 .AND. $ ILAENV( 11, 'ZSTEQR', 'N', 1, 0, 0, 0 ).EQ.1 ) THEN N = 2 RONE = ONE RNAN = SQRT( -RONE ) DO 380 JNAN = 1, 2 - DO 370 JROUT = 1, 3 + DO 370 JROUT = 1, 5 SD( 1 ) = ONE SD( 2 ) = ONE + ONE SE( 1 ) = ONE SE( 2 ) = ZERO SD( JNAN ) = RNAN - IF( JROUT.EQ.1 ) THEN - RNAME = 'ZSTEQR' - CALL ZSTEQR( 'N', N, SD, SE, Z, LDU, RWORK, IINFO ) + OPT = 'N' + IF( JROUT.EQ.4 ) OPT = 'I' + IF( JROUT.EQ.5 ) OPT = 'V' + IF( JROUT.EQ.1 .OR. JROUT.EQ.4 ) THEN + RNAME = 'ZSTEQR('//OPT//')' + CALL ZSTEQR( OPT, N, SD, SE, Z, LDU, RWORK, IINFO ) ELSE IF( JROUT.EQ.2 ) THEN RNAME = 'DSTERF' CALL DSTERF( N, SD, SE, IINFO ) ELSE - RNAME = 'ZSTEMR' + RNAME = 'ZSTEMR('//OPT//')' VL = ZERO VU = ZERO IL = 0 IU = 0 TRYRAC = .TRUE. - CALL ZSTEMR( 'N', 'A', N, SD, SE, VL, VU, IL, IU, M, + CALL ZSTEMR( OPT, 'A', N, SD, SE, VL, VU, IL, IU, M, $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, $ RWORK, LRWORK, $ IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) @@ -2010,7 +2015,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT ) RETURN * - 9983 FORMAT( ' ZCHKST: ', A6, ' returned a finite eigenvalue for a', + 9983 FORMAT( ' ZCHKST: ', A9, ' returned a finite eigenvalue for a', $ ' 2 by 2 matrix with a NaN at D(', I1, ')' ) 9999 FORMAT( ' ZCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )