diff --git a/SRC/dlaebz.f b/SRC/dlaebz.f index 24a611643..57908c349 100644 --- a/SRC/dlaebz.f +++ b/SRC/dlaebz.f @@ -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 diff --git a/SRC/dstebz.f b/SRC/dstebz.f index 9b228a9e3..d56f6a9df 100644 --- a/SRC/dstebz.f +++ b/SRC/dstebz.f @@ -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 diff --git a/SRC/slaebz.f b/SRC/slaebz.f index 58c913ef1..fc5071426 100644 --- a/SRC/slaebz.f +++ b/SRC/slaebz.f @@ -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 diff --git a/SRC/sstebz.f b/SRC/sstebz.f index b4a0bce00..ef7e78ee6 100644 --- a/SRC/sstebz.f +++ b/SRC/sstebz.f @@ -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 diff --git a/TESTING/EIG/cchkee.F b/TESTING/EIG/cchkee.F index 6b2bca670..77aa17a80 100644 --- a/TESTING/EIG/cchkee.F +++ b/TESTING/EIG/cchkee.F @@ -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 ) diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 790c7f5f1..e53248b61 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -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 @@ -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 .. * @@ -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 @@ -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 * @@ -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 * @@ -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.) diff --git a/TESTING/EIG/dchkee.F b/TESTING/EIG/dchkee.F index 07bb24e9f..b3f894350 100644 --- a/TESTING/EIG/dchkee.F +++ b/TESTING/EIG/dchkee.F @@ -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 ) diff --git a/TESTING/EIG/dchkst.f b/TESTING/EIG/dchkst.f index 012aa95d4..6048eb0df 100644 --- a/TESTING/EIG/dchkst.f +++ b/TESTING/EIG/dchkst.f @@ -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 @@ -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 .. * @@ -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 @@ -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 * @@ -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 * @@ -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.) diff --git a/TESTING/EIG/schkee.F b/TESTING/EIG/schkee.F index 4fa913e30..a7f774eb2 100644 --- a/TESTING/EIG/schkee.F +++ b/TESTING/EIG/schkee.F @@ -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 ) diff --git a/TESTING/EIG/schkst.f b/TESTING/EIG/schkst.f index c5c7f57c3..a3f1366d9 100644 --- a/TESTING/EIG/schkst.f +++ b/TESTING/EIG/schkst.f @@ -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 @@ -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 .. * @@ -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 @@ -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 * @@ -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 * @@ -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.) diff --git a/TESTING/EIG/zchkee.F b/TESTING/EIG/zchkee.F index 3dc9ab4d8..dd5fd7464 100644 --- a/TESTING/EIG/zchkee.F +++ b/TESTING/EIG/zchkee.F @@ -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 ) diff --git a/TESTING/EIG/zchkst.f b/TESTING/EIG/zchkst.f index 4335e15f0..11ac55961 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -633,8 +633,10 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ CONE = ( 1.0D+0, 0.0D+0 ) ) DOUBLE PRECISION HALF PARAMETER ( HALF = ONE / TWO ) + DOUBLE PRECISION TENTH + PARAMETER ( TENTH = 0.1D0 ) INTEGER MAXTYP - PARAMETER ( MAXTYP = 21 ) + PARAMETER ( MAXTYP = 22 ) LOGICAL CRANGE PARAMETER ( CRANGE = .FALSE. ) LOGICAL CREL @@ -673,11 +675,11 @@ SUBROUTINE ZCHKST( 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 .. * @@ -806,7 +808,7 @@ SUBROUTINE ZCHKST( 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 @@ -819,6 +821,15 @@ SUBROUTINE ZCHKST( 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 * @@ -900,6 +911,20 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, A( I, I-1 ) = DCONJG( 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 * @@ -1432,6 +1457,13 @@ SUBROUTINE ZCHKST( 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 ZSTEIN to compute eigenvectors corresponding to * eigenvalues in WA1. (First call DSTEBZ again, to make sure * it returns these eigenvalues in the correct order.) diff --git a/TESTING/sep.in b/TESTING/sep.in index 99b47c0aa..60442870f 100644 --- a/TESTING/sep.in +++ b/TESTING/sep.in @@ -10,6 +10,6 @@ T Put T to test the LAPACK routines T Put T to test the driver routines T Put T to test the error exits 1 Code to interpret the seed -SEP 20 -1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 21 +SEP 21 +1 2 3 4 5 6 7 8 10 11 12 13 14 15 16 17 18 19 20 21 22