Skip to content

Drop the absolute floor from the splitting test in xSTEBZ #1385

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 1 commit into
base: master
Choose a base branch
from
Open

Drop the absolute floor from the splitting test in xSTEBZ #1385

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
2 changes: 1 addition & 1 deletion SRC/dlaebz.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@
*> Note that the intervals are in all cases half-open intervals,
*> i.e., of the form (a,b] , which includes b but not a .
*>
*> To avoid underflow, the matrix should be scaled so that its largest
*> To avoid overflow, the matrix should be scaled so that its largest
*> element is no greater than overflow**(1/2) * underflow**(1/4)
*> in absolute value. To assure the most accurate computation
*> of small eigenvalues, the matrix should be scaled to be
Expand Down
11 changes: 10 additions & 1 deletion SRC/dstebz.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -419,7 +419,16 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D,
*
DO 10 J = 2, N
TMP1 = E( J-1 )**2
IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
* Split where e(j)**2 < ulp**2 |d(j) d(j+1)| (Kahan's relative
* criterion, as in xLARRA) or e(j) = 0. No absolute floor: the
* drivers scale small matrices to sqrt(SAFEMN/ULP), where a floor
* of SAFEMN on e(j)**2 would discard off-diagonals as large as
* sqrt(ULP) relative to the matrix. The threshold is formed as
* a product of scaled factors so that it cannot overflow before
* e(j)**2 does.
*
TMP2 = ( ABS( D( J ) )*ULP )*( ABS( D( J-1 ) )*ULP )
IF( TMP1.LT.TMP2 .OR. TMP1.EQ.ZERO ) THEN
ISPLIT( NSPLIT ) = J - 1
NSPLIT = NSPLIT + 1
WORK( J-1 ) = ZERO
Expand Down
2 changes: 1 addition & 1 deletion SRC/slaebz.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@
*> Note that the intervals are in all cases half-open intervals,
*> i.e., of the form (a,b] , which includes b but not a .
*>
*> To avoid underflow, the matrix should be scaled so that its largest
*> To avoid overflow, the matrix should be scaled so that its largest
*> element is no greater than overflow**(1/2) * underflow**(1/4)
*> in absolute value. To assure the most accurate computation
*> of small eigenvalues, the matrix should be scaled to be
Expand Down
11 changes: 10 additions & 1 deletion SRC/sstebz.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -418,7 +418,16 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D,
*
DO 10 J = 2, N
TMP1 = E( J-1 )**2
IF( ABS( D( J )*D( J-1 ) )*ULP**2+SAFEMN.GT.TMP1 ) THEN
* Split where e(j)**2 < ulp**2 |d(j) d(j+1)| (Kahan's relative
* criterion, as in xLARRA) or e(j) = 0. No absolute floor: the
* drivers scale small matrices to sqrt(SAFEMN/ULP), where a floor
* of SAFEMN on e(j)**2 would discard off-diagonals as large as
* sqrt(ULP) relative to the matrix. The threshold is formed as
* a product of scaled factors so that it cannot overflow before
* e(j)**2 does.
*
TMP2 = ( ABS( D( J ) )*ULP )*( ABS( D( J-1 ) )*ULP )
IF( TMP1.LT.TMP2 .OR. TMP1.EQ.ZERO ) THEN
ISPLIT( NSPLIT ) = J - 1
NSPLIT = NSPLIT + 1
WORK( J-1 ) = ZERO
Expand Down
2 changes: 1 addition & 1 deletion TESTING/EIG/cchkee.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1861,7 +1861,7 @@ PROGRAM CCHKEE
* NBMIN = minimum block size
* NX = crossover point
*
MAXTYP = 21
MAXTYP = 22
NTYPES = MIN( MAXTYP, NTYPES )
CALL ALAREQ( C3, NTYPES, DOTYPE, MAXTYP, NIN, NOUT )
CALL XLAENV( 1, 1 )
Expand Down
42 changes: 37 additions & 5 deletions TESTING/EIG/cchkst.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -633,8 +633,10 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ CONE = ( 1.0E+0, 0.0E+0 ) )
REAL HALF
PARAMETER ( HALF = ONE / TWO )
REAL TENTH
PARAMETER ( TENTH = 0.1E0 )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 21 )
PARAMETER ( MAXTYP = 22 )
LOGICAL CRANGE
PARAMETER ( CRANGE = .FALSE. )
LOGICAL CREL
Expand Down Expand Up @@ -673,11 +675,11 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
* ..
* .. Data statements ..
DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
$ 8, 8, 9, 9, 9, 9, 9, 10 /
$ 8, 8, 9, 9, 9, 9, 9, 10, 11 /
DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
$ 2, 3, 1, 1, 1, 2, 3, 1 /
$ 2, 3, 1, 1, 1, 2, 3, 1, 4 /
DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
$ 0, 0, 4, 3, 1, 4, 4, 3 /
$ 0, 0, 4, 3, 1, 4, 4, 3, 0 /
* ..
* .. Executable Statements ..
*
Expand Down Expand Up @@ -806,7 +808,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
* Compute norm
*
GO TO ( 40, 50, 60 )KMAGN( JTYPE )
GO TO ( 40, 50, 60, 65 )KMAGN( JTYPE )
*
40 CONTINUE
ANORM = ONE
Expand All @@ -819,6 +821,15 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
60 CONTINUE
ANORM = RTUNFL*REAL( N )*ULPINV
GO TO 70
*
65 CONTINUE
*
* The smallest norm the eigenvalue drivers scale a matrix up
* to, below which the splitting test of SSTEBZ has to be
* relative to see an off-diagonal entry.
*
ANORM = RTUNFL / SQRT( ULP )
GO TO 70
*
70 CONTINUE
*
Expand Down Expand Up @@ -900,6 +911,20 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
A( I, I-1 ) = CONJG( A( I-1, I ) )
END IF
90 CONTINUE
*
ELSE IF( ITYPE.EQ.11 ) THEN
*
* Tridiagonal with equal diagonal entries and an
* off-diagonal entry whose square underflows, but which is
* far above the relative splitting threshold.
*
DO 95 JC = 1, N
A( JC, JC ) = ANORM
95 CONTINUE
DO 96 JC = 1, N - 1
A( JC, JC+1 ) = ANORM*( TENTH*SQRT( ULP ) )
A( JC+1, JC ) = A( JC, JC+1 )
96 CONTINUE
*
ELSE
*
Expand Down Expand Up @@ -1432,6 +1457,13 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
*
* The eigenvalues of type 22 are closer together than the
* floor the interval above is built with, so its 'V' range
* holds the whole spectrum and cannot match the 'I' range.
*
IF( JTYPE.EQ.22 )
$ RESULT( 19 ) = ZERO
*
* Call CSTEIN to compute eigenvectors corresponding to
* eigenvalues in WA1. (First call SSTEBZ again, to make sure
* it returns these eigenvalues in the correct order.)
Expand Down
2 changes: 1 addition & 1 deletion TESTING/EIG/dchkee.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1866,7 +1866,7 @@ PROGRAM DCHKEE
* NBMIN = minimum block size
* NX = crossover point
*
MAXTYP = 21
MAXTYP = 22
NTYPES = MIN( MAXTYP, NTYPES )
CALL ALAREQ( C3, NTYPES, DOTYPE, MAXTYP, NIN, NOUT )
CALL XLAENV( 1, 1 )
Expand Down
42 changes: 37 additions & 5 deletions TESTING/EIG/dchkst.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -617,8 +617,10 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ EIGHT = 8.0D0, TEN = 10.0D0, HUN = 100.0D0 )
DOUBLE PRECISION HALF
PARAMETER ( HALF = ONE / TWO )
DOUBLE PRECISION TENTH
PARAMETER ( TENTH = 0.1D0 )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 21 )
PARAMETER ( MAXTYP = 22 )
LOGICAL SRANGE
PARAMETER ( SRANGE = .FALSE. )
LOGICAL SREL
Expand Down Expand Up @@ -656,11 +658,11 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
* ..
* .. Data statements ..
DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
$ 8, 8, 9, 9, 9, 9, 9, 10 /
$ 8, 8, 9, 9, 9, 9, 9, 10, 11 /
DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
$ 2, 3, 1, 1, 1, 2, 3, 1 /
$ 2, 3, 1, 1, 1, 2, 3, 1, 4 /
DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
$ 0, 0, 4, 3, 1, 4, 4, 3 /
$ 0, 0, 4, 3, 1, 4, 4, 3, 0 /
* ..
* .. Executable Statements ..
*
Expand Down Expand Up @@ -787,7 +789,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
* Compute norm
*
GO TO ( 40, 50, 60 )KMAGN( JTYPE )
GO TO ( 40, 50, 60, 65 )KMAGN( JTYPE )
*
40 CONTINUE
ANORM = ONE
Expand All @@ -800,6 +802,15 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
60 CONTINUE
ANORM = RTUNFL*N*ULPINV
GO TO 70
*
65 CONTINUE
*
* The smallest norm the eigenvalue drivers scale a matrix up
* to, below which the splitting test of DSTEBZ has to be
* relative to see an off-diagonal entry.
*
ANORM = RTUNFL / SQRT( ULP )
GO TO 70
*
70 CONTINUE
*
Expand Down Expand Up @@ -885,6 +896,20 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
A( I, I-1 ) = A( I-1, I )
END IF
90 CONTINUE
*
ELSE IF( ITYPE.EQ.11 ) THEN
*
* Tridiagonal with equal diagonal entries and an
* off-diagonal entry whose square underflows, but which is
* far above the relative splitting threshold.
*
DO 95 JC = 1, N
A( JC, JC ) = ANORM
95 CONTINUE
DO 96 JC = 1, N - 1
A( JC, JC+1 ) = ANORM*( TENTH*SQRT( ULP ) )
A( JC+1, JC ) = A( JC, JC+1 )
96 CONTINUE
*
ELSE
*
Expand Down Expand Up @@ -1416,6 +1441,13 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
*
* The eigenvalues of type 22 are closer together than the
* floor the interval above is built with, so its 'V' range
* holds the whole spectrum and cannot match the 'I' range.
*
IF( JTYPE.EQ.22 )
$ RESULT( 19 ) = ZERO
*
* Call DSTEIN to compute eigenvectors corresponding to
* eigenvalues in WA1. (First call DSTEBZ again, to make sure
* it returns these eigenvalues in the correct order.)
Expand Down
2 changes: 1 addition & 1 deletion TESTING/EIG/schkee.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1867,7 +1867,7 @@ PROGRAM SCHKEE
* NBMIN = minimum block size
* NX = crossover point
*
MAXTYP = 21
MAXTYP = 22
NTYPES = MIN( MAXTYP, NTYPES )
CALL ALAREQ( C3, NTYPES, DOTYPE, MAXTYP, NIN, NOUT )
CALL XLAENV( 1, 1 )
Expand Down
42 changes: 37 additions & 5 deletions TESTING/EIG/schkst.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -617,8 +617,10 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
$ EIGHT = 8.0E0, TEN = 10.0E0, HUN = 100.0E0 )
REAL HALF
PARAMETER ( HALF = ONE / TWO )
REAL TENTH
PARAMETER ( TENTH = 0.1E0 )
INTEGER MAXTYP
PARAMETER ( MAXTYP = 21 )
PARAMETER ( MAXTYP = 22 )
LOGICAL SRANGE
PARAMETER ( SRANGE = .FALSE. )
LOGICAL SREL
Expand Down Expand Up @@ -656,11 +658,11 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
* ..
* .. Data statements ..
DATA KTYPE / 1, 2, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 8,
$ 8, 8, 9, 9, 9, 9, 9, 10 /
$ 8, 8, 9, 9, 9, 9, 9, 10, 11 /
DATA KMAGN / 1, 1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
$ 2, 3, 1, 1, 1, 2, 3, 1 /
$ 2, 3, 1, 1, 1, 2, 3, 1, 4 /
DATA KMODE / 0, 0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
$ 0, 0, 4, 3, 1, 4, 4, 3 /
$ 0, 0, 4, 3, 1, 4, 4, 3, 0 /
* ..
* .. Executable Statements ..
*
Expand Down Expand Up @@ -787,7 +789,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
* Compute norm
*
GO TO ( 40, 50, 60 )KMAGN( JTYPE )
GO TO ( 40, 50, 60, 65 )KMAGN( JTYPE )
*
40 CONTINUE
ANORM = ONE
Expand All @@ -800,6 +802,15 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
60 CONTINUE
ANORM = RTUNFL*REAL( N )*ULPINV
GO TO 70
*
65 CONTINUE
*
* The smallest norm the eigenvalue drivers scale a matrix up
* to, below which the splitting test of SSTEBZ has to be
* relative to see an off-diagonal entry.
*
ANORM = RTUNFL / SQRT( ULP )
GO TO 70
*
70 CONTINUE
*
Expand Down Expand Up @@ -885,6 +896,20 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
A( I, I-1 ) = A( I-1, I )
END IF
90 CONTINUE
*
ELSE IF( ITYPE.EQ.11 ) THEN
*
* Tridiagonal with equal diagonal entries and an
* off-diagonal entry whose square underflows, but which is
* far above the relative splitting threshold.
*
DO 95 JC = 1, N
A( JC, JC ) = ANORM
95 CONTINUE
DO 96 JC = 1, N - 1
A( JC, JC+1 ) = ANORM*( TENTH*SQRT( ULP ) )
A( JC+1, JC ) = A( JC, JC+1 )
96 CONTINUE
*
ELSE
*
Expand Down Expand Up @@ -1416,6 +1441,13 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*
RESULT( 19 ) = ( TEMP1+TEMP2 ) / MAX( UNFL, TEMP3*ULP )
*
* The eigenvalues of type 22 are closer together than the
* floor the interval above is built with, so its 'V' range
* holds the whole spectrum and cannot match the 'I' range.
*
IF( JTYPE.EQ.22 )
$ RESULT( 19 ) = ZERO
*
* Call SSTEIN to compute eigenvectors corresponding to
* eigenvalues in WA1. (First call SSTEBZ again, to make sure
* it returns these eigenvalues in the correct order.)
Expand Down
2 changes: 1 addition & 1 deletion TESTING/EIG/zchkee.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1861,7 +1861,7 @@ PROGRAM ZCHKEE
* NBMIN = minimum block size
* NX = crossover point
*
MAXTYP = 21
MAXTYP = 22
NTYPES = MIN( MAXTYP, NTYPES )
CALL ALAREQ( C3, NTYPES, DOTYPE, MAXTYP, NIN, NOUT )
CALL XLAENV( 1, 1 )
Expand Down
Loading
Loading