Skip to content

Propagate a NaN through xLAE2 and xLAEV2 instead of returning finite eigenvalues #1405

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
rmlarsen wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Propagate a NaN through xLAE2 and xLAEV2 instead of returning finite eigenvalues #1405

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6f96b00
Propagate a NaN through xLAE2 and xLAEV2 instead of returning finite …
rmlarsen Sep 10, 2026
349f09f
TESTING: Cover LAEV2 NaN propagation with eigenvectors
rmlarsen Sep 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion SRC/dlae2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down
8 changes: 7 additions & 1 deletion SRC/dlaev2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down
8 changes: 7 additions & 1 deletion SRC/slae2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down
8 changes: 7 additions & 1 deletion SRC/slaev2.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 )
Expand Down
62 changes: 61 additions & 1 deletion TESTING/EIG/cchkst.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -651,16 +651,21 @@ 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 ),
$ KTYPE( MAXTYP )
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,
Expand Down Expand Up @@ -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, ')' )
*
Expand Down
62 changes: 61 additions & 1 deletion TESTING/EIG/dchkst.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -634,16 +634,21 @@ 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 ),
$ KTYPE( MAXTYP )
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,
Expand Down Expand Up @@ -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, ')' )
*
Expand Down
62 changes: 61 additions & 1 deletion TESTING/EIG/schkst.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -634,16 +634,21 @@ 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 ),
$ KTYPE( MAXTYP )
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,
Expand Down Expand Up @@ -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, ')' )
*
Expand Down
Loading
Loading