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..01329947f 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -651,6 +651,10 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER OPT + CHARACTER*9 RNAME + INTEGER JNAN, JROUT + REAL RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -658,9 +662,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 +1958,66 @@ 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). Eigenvalue-only calls test +* xLAE2; eigenvector-producing calls also test xLAEV2. +* + 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, 5 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + 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('//OPT//')' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + 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 ) + 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: ', 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 012aa95d4..a8b3433b7 100644 --- a/TESTING/EIG/dchkst.f +++ b/TESTING/EIG/dchkst.f @@ -634,6 +634,10 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER OPT + CHARACTER*9 RNAME + INTEGER JNAN, JROUT + DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -641,9 +645,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 +1936,66 @@ 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). 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. + $ 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, 5 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + 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('//OPT//')' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + 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 ) + 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: ', 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 c5c7f57c3..bdeb6b8eb 100644 --- a/TESTING/EIG/schkst.f +++ b/TESTING/EIG/schkst.f @@ -634,6 +634,10 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER OPT + CHARACTER*9 RNAME + INTEGER JNAN, JROUT + REAL RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -641,9 +645,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 +1936,66 @@ 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). 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. + $ 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, 5 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + 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('//OPT//')' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + 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 ) + 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: ', 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 4335e15f0..6b7674fef 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -651,6 +651,10 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU * .. + CHARACTER OPT + CHARACTER*9 RNAME + INTEGER JNAN, JROUT + DOUBLE PRECISION RNAN, RONE * .. Local Arrays .. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ), $ KMAGN( MAXTYP ), KMODE( MAXTYP ), @@ -658,9 +662,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 +1957,66 @@ 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). Eigenvalue-only calls test +* xLAE2; eigenvector-producing calls also test xLAEV2. +* + 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, 5 + SD( 1 ) = ONE + SD( 2 ) = ONE + ONE + SE( 1 ) = ONE + SE( 2 ) = ZERO + SD( JNAN ) = RNAN + 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('//OPT//')' + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + 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 ) + 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: ', 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, ')' ) *