From 66f0863d77d0e7eb219c89a48e19d43e766a6733 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 21:20:04 -0700 Subject: [PATCH] Drop the absolute floor from the splitting test in xSTEBZ xSTEBZ declares an off-diagonal entry negligible when e(j)**2 < ulp**2 |d(j) d(j+1)| + SAFEMN. The SAFEMN term makes the test absolute: any |e(j)| below sqrt(SAFEMN) (1.5e-154 in double precision) splits the matrix whatever its neighbours are. The drivers built on xSTEBZ (xSTEVX, xSYEVX, xHEEVX, xSPEVX, xHPEVX, xSBEVX, xHBEVX, and xSTEVR/xSYEVR/xHEEVR for RANGE = 'V') scale a small matrix up only to RMIN = sqrt(SAFMIN/ULP), about 1e-146, where that floor corresponds to |e(j)|/|T| < sqrt(ULP) = 1.5e-8. So for a matrix whose largest entry is below sqrt(SAFMIN)/ULP (6.7e-139) the drivers drop off-diagonals of relative size between ULP and 1.5e-8: the eigenvalues are wrong by up to 1.5e-8 |T| and the eigenvectors have residuals of the same size, with INFO = 0, while xSTEV and xSYEV return them to working precision. The test suite never reaches this regime; its smallest matrix scale is RTUNFL*N/ULP, exactly where the floor stops mattering. Use Kahan's relative criterion e(j)**2 < ulp**2 |d(j) d(j+1)|, or e(j) = 0, as xLARRA already does, with the threshold formed as (|d(j)| ulp)(|d(j+1)| ulp) so that it cannot overflow before e(j)**2 does. For a matrix in the normal range this changes the treatment only of an off-diagonal below 1.5e-154 sitting between two diagonal entries below 2e-138, which now stays in its block. Also correct "underflow" to "overflow" in the scaling advice of the xLAEBZ documentation. For T = s tridiag(1, r, 1) with the eigenvalues known in closed form, every DSTEBZ-based driver returns eigenvalue errors of r |T| for r between 1e-14 and 1e-8 at s <= 1e-140 on the parent commit and errors below n ulp |T| on this branch, at every scale. xCHKST gets a matrix type for it: type 22 is a tridiagonal matrix with equal diagonal entries and an off-diagonal entry one tenth of the square root of an ulp of them, at the smallest norm the drivers scale a matrix up to. Its off-diagonal square underflows, so the absolute term of the old test hid it, while the relative test does not. The existing comparison of xSTEBZ against xSTERF, and the residual of the eigenvectors xSTEIN builds from its blocks, both catch the difference: on the parent commit the type fails 7 ratios in single precision and 8 in double at ratios of 1e6 and above, against a threshold of 60. The 'I' against 'V' comparison is skipped for the type, because the interval the test builds is floored at twice the square root of the underflow threshold, which is wider than the whole spectrum here. The complex checkers carry the same type, since they call the real xSTEBZ. A power-of-two twin sweep of DSTEBZ + DSTEIN over 42 exponents from -1070 to 1020 is identical to the parent wherever the parent was correct. The full LAPACK test suite passes: 5446341 LAPACK tests, 0 numerical errors, 0 other errors; the 4440 tests above the parent are the new type. Co-Authored-By: Claude Fable 5.1 --- SRC/dlaebz.f | 2 +- SRC/dstebz.f | 11 ++++++++++- SRC/slaebz.f | 2 +- SRC/sstebz.f | 11 ++++++++++- TESTING/EIG/cchkee.F | 2 +- TESTING/EIG/cchkst.f | 42 +++++++++++++++++++++++++++++++++++++----- TESTING/EIG/dchkee.F | 2 +- TESTING/EIG/dchkst.f | 42 +++++++++++++++++++++++++++++++++++++----- TESTING/EIG/schkee.F | 2 +- TESTING/EIG/schkst.f | 42 +++++++++++++++++++++++++++++++++++++----- TESTING/EIG/zchkee.F | 2 +- TESTING/EIG/zchkst.f | 42 +++++++++++++++++++++++++++++++++++++----- TESTING/sep.in | 4 ++-- 13 files changed, 176 insertions(+), 30 deletions(-) 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