From b364dad0b36219ccac8097ab69c3f72f4267652d Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 21:19:21 -0700 Subject: [PATCH] Scale ABSTOL with the matrix in xSTEVX and xSTEVR xSTEVX and xSTEVR scale the tridiagonal matrix into [RMIN, RMAX] before calling xSTEBZ but, unlike the other 28 drivers with the same scaling block, pass the caller's ABSTOL unscaled. A positive ABSTOL is then applied to a matrix multiplied by SIGMA = RMAX / |T| whenever |T| > RMAX = SAFMIN**(-1/4), which is 1.6e77 in double precision and 3e9 in single, so bisection stops as soon as an interval is narrower than ABSTOL / SIGMA in the caller's units, a width that can exceed the norm of the scaled matrix. SSTEVX on tridiag(1e12, 5e11, 1e12) with ABSTOL = 1e-6 |T| returns eigenvalues with relative errors of 2e-4, 3e-2 at 1e14, and DSTEVX at |T| = 1e100 with ABSTOL = 1e-10 |T| returns 0.4, all with INFO = 0, while SSYEVX and DSYEVX honour the tolerance. Below RMIN the tolerance is only applied more tightly than requested. Scale a positive ABSTOL with the matrix and pass ABSTLL = MAX( MIN( ABSTOL, |T| )*SIGMA, SAFMIN ) to xSTEBZ. The plain product of xSYEVX is not enough: for the documented "most accurate" ABSTOL = 2*SAFMIN it underflows to zero once |T| > RMAX/EPS, and xSTEBZ reads a zero tolerance as a request for its default EPS*|T|, so DSTEVX on D = (1e100, 1), E = 5e49 would return 1.1e83 for the eigenvalue 0.75. Every positive tolerance below PIVMIN >= SAFMIN is equivalent in xLAEBZ, so the floor changes nothing that did not underflow; a tolerance above |T| is met by the Gershgorin interval, so the cap changes nothing that could not overflow when a tiny matrix is scaled up. Nothing changes for ABSTOL <= 0 or for a matrix the driver does not scale. xDRVST gets a matrix type for it: type 19 is a tridiagonal matrix with equal diagonal entries and half that off the diagonal, scaled to the "large" magnitude of the existing types, which is above the norm the drivers scale down from, and it is the only type that asks for a positive ABSTOL, ten ulps of the matrix norm. Every other type passes 2*SAFMIN, which is inert in the bisection and so cannot show the defect, and the tridiagonal drivers were tested on diagonal matrices alone, where each block is one by one and the tolerance never matters. On the parent commit the type fails 18 ratios in single precision and 40 in double. With the change xSTEVX and xSYEVX agree to the last digit on the matrices above at every scale, and ABSTOL = 2*SAFMIN gives the same eigenvalues as on master. The full LAPACK test suite passes: 5449821 LAPACK tests, 0 numerical errors, 0 other errors; the 7920 tests above the parent are the new type. Co-Authored-By: Claude Fable 5.1 --- SRC/dstevr.f | 14 ++++++++++++-- SRC/dstevx.f | 14 ++++++++++++-- SRC/sstevr.f | 14 ++++++++++++-- SRC/sstevx.f | 14 ++++++++++++-- TESTING/EIG/dchkee.F | 2 +- TESTING/EIG/ddrvst.f | 37 ++++++++++++++++++++++++++++++------- TESTING/EIG/schkee.F | 2 +- TESTING/EIG/sdrvst.f | 37 ++++++++++++++++++++++++++++++------- 8 files changed, 110 insertions(+), 24 deletions(-) diff --git a/SRC/dstevr.f b/SRC/dstevr.f index d75afcb37..64e09c9fd 100644 --- a/SRC/dstevr.f +++ b/SRC/dstevr.f @@ -329,7 +329,8 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP, $ INDIWO, ISCALE, ITMP1, J, JJ, LIWMIN, LWMIN, $ NSPLIT - DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, + DOUBLE PRECISION ABSTLL, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, + $ SMLNUM, $ TMP1, TNRM, VLL, VUU * .. * .. External Functions .. @@ -440,6 +441,7 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, * Scale matrix to allowable range, if necessary. * ISCALE = 0 + ABSTLL = ABSTOL IF( VALEIG ) THEN VLL = VL VUU = VU @@ -456,6 +458,14 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, IF( ISCALE.EQ.1 ) THEN CALL DSCAL( N, SIGMA, D, 1 ) CALL DSCAL( N-1, SIGMA, E( 1 ), 1 ) +* +* Scale a positive ABSTOL with the matrix. Keep the product +* positive, as DSTEBZ takes ABSTOL <= 0 as a request for its +* default tolerance EPS*|T|, and cap ABSTOL at |T| so that the +* product cannot overflow. +* + IF( ABSTOL.GT.ZERO ) + $ ABSTLL = MAX( MIN( ABSTOL, TNRM )*SIGMA, SAFMIN ) IF( VALEIG ) THEN VLL = VL*SIGMA VUU = VU*SIGMA @@ -522,7 +532,7 @@ SUBROUTINE DSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ORDER = 'E' END IF - CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, + CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, D, E, $ M, $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK, $ IWORK( INDIWO ), INFO ) diff --git a/SRC/dstevx.f b/SRC/dstevx.f index e18c03265..4c8b61ff0 100644 --- a/SRC/dstevx.f +++ b/SRC/dstevx.f @@ -253,7 +253,8 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, CHARACTER ORDER INTEGER I, IINFO, IMAX, INDISP, INDIWO, INDWRK, $ ISCALE, ITMP1, J, JJ, NSPLIT - DOUBLE PRECISION BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, + DOUBLE PRECISION ABSTLL, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, + $ SMLNUM, $ TMP1, TNRM, VLL, VUU * .. * .. External Functions .. @@ -340,6 +341,7 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, * Scale matrix to allowable range, if necessary. * ISCALE = 0 + ABSTLL = ABSTOL IF( VALEIG ) THEN VLL = VL VUU = VU @@ -358,6 +360,14 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, IF( ISCALE.EQ.1 ) THEN CALL DSCAL( N, SIGMA, D, 1 ) CALL DSCAL( N-1, SIGMA, E( 1 ), 1 ) +* +* Scale a positive ABSTOL with the matrix. Keep the product +* positive, as DSTEBZ takes ABSTOL <= 0 as a request for its +* default tolerance EPS*|T|, and cap ABSTOL at |T| so that the +* product cannot overflow. +* + IF( ABSTOL.GT.ZERO ) + $ ABSTLL = MAX( MIN( ABSTOL, TNRM )*SIGMA, SAFMIN ) IF( VALEIG ) THEN VLL = VL*SIGMA VUU = VU*SIGMA @@ -406,7 +416,7 @@ SUBROUTINE DSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, INDWRK = 1 INDISP = 1 + N INDIWO = INDISP + N - CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, + CALL DSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, D, E, $ M, $ NSPLIT, W, IWORK( 1 ), IWORK( INDISP ), $ WORK( INDWRK ), IWORK( INDIWO ), IINFO ) diff --git a/SRC/sstevr.f b/SRC/sstevr.f index 0b72c2f0e..750a00429 100644 --- a/SRC/sstevr.f +++ b/SRC/sstevr.f @@ -330,7 +330,8 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, CHARACTER ORDER INTEGER I, IEEEOK, IMAX, INDIBL, INDIFL, INDISP, $ INDIWO, ISCALE, J, JJ, LIWMIN, LWMIN, NSPLIT - REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, + REAL ABSTLL, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, + $ SMLNUM, $ TMP1, TNRM, VLL, VUU * .. * .. External Functions .. @@ -442,6 +443,7 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, * Scale matrix to allowable range, if necessary. * ISCALE = 0 + ABSTLL = ABSTOL IF( VALEIG ) THEN VLL = VL VUU = VU @@ -458,6 +460,14 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, IF( ISCALE.EQ.1 ) THEN CALL SSCAL( N, SIGMA, D, 1 ) CALL SSCAL( N-1, SIGMA, E( 1 ), 1 ) +* +* Scale a positive ABSTOL with the matrix. Keep the product +* positive, as SSTEBZ takes ABSTOL <= 0 as a request for its +* default tolerance EPS*|T|, and cap ABSTOL at |T| so that the +* product cannot overflow. +* + IF( ABSTOL.GT.ZERO ) + $ ABSTLL = MAX( MIN( ABSTOL, TNRM )*SIGMA, SAFMIN ) IF( VALEIG ) THEN VLL = VL*SIGMA VUU = VU*SIGMA @@ -524,7 +534,7 @@ SUBROUTINE SSTEVR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ORDER = 'E' END IF - CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, + CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, D, E, $ M, $ NSPLIT, W, IWORK( INDIBL ), IWORK( INDISP ), WORK, $ IWORK( INDIWO ), INFO ) diff --git a/SRC/sstevx.f b/SRC/sstevx.f index f84e740ef..461ff899d 100644 --- a/SRC/sstevx.f +++ b/SRC/sstevx.f @@ -253,7 +253,8 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, CHARACTER ORDER INTEGER I, IINFO, IMAX, INDISP, INDIWO, INDWRK, $ ISCALE, ITMP1, J, JJ, NSPLIT - REAL BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, SMLNUM, + REAL ABSTLL, BIGNUM, EPS, RMAX, RMIN, SAFMIN, SIGMA, + $ SMLNUM, $ TMP1, TNRM, VLL, VUU * .. * .. External Functions .. @@ -340,6 +341,7 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, * Scale matrix to allowable range, if necessary. * ISCALE = 0 + ABSTLL = ABSTOL IF ( VALEIG ) THEN VLL = VL VUU = VU @@ -358,6 +360,14 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, IF( ISCALE.EQ.1 ) THEN CALL SSCAL( N, SIGMA, D, 1 ) CALL SSCAL( N-1, SIGMA, E( 1 ), 1 ) +* +* Scale a positive ABSTOL with the matrix. Keep the product +* positive, as SSTEBZ takes ABSTOL <= 0 as a request for its +* default tolerance EPS*|T|, and cap ABSTOL at |T| so that the +* product cannot overflow. +* + IF( ABSTOL.GT.ZERO ) + $ ABSTLL = MAX( MIN( ABSTOL, TNRM )*SIGMA, SAFMIN ) IF( VALEIG ) THEN VLL = VL*SIGMA VUU = VU*SIGMA @@ -406,7 +416,7 @@ SUBROUTINE SSTEVX( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, INDWRK = 1 INDISP = 1 + N INDIWO = INDISP + N - CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTOL, D, E, + CALL SSTEBZ( RANGE, ORDER, N, VLL, VUU, IL, IU, ABSTLL, D, E, $ M, $ NSPLIT, W, IWORK( 1 ), IWORK( INDISP ), $ WORK( INDWRK ), IWORK( INDIWO ), IINFO ) diff --git a/TESTING/EIG/dchkee.F b/TESTING/EIG/dchkee.F index 07bb24e9f..937b99ebd 100644 --- a/TESTING/EIG/dchkee.F +++ b/TESTING/EIG/dchkee.F @@ -1924,7 +1924,7 @@ PROGRAM DCHKEE $ A( 1, 3 ), D( 1, 12 ), A( 1, 4 ), WORK, $ LWORK, IWORK, LIWORK, RESULT, INFO ) ELSE - CALL DDRVST( NN, NVAL, 18, DOTYPE, ISEED, THRESH, NOUT, + CALL DDRVST( NN, NVAL, 19, DOTYPE, ISEED, THRESH, NOUT, $ A( 1, 1 ), NMAX, D( 1, 3 ), D( 1, 4 ), $ D( 1, 5 ), D( 1, 6 ), D( 1, 8 ), D( 1, 9 ), $ D( 1, 10 ), D( 1, 11 ), A( 1, 2 ), NMAX, diff --git a/TESTING/EIG/ddrvst.f b/TESTING/EIG/ddrvst.f index fbd3ff0f7..c0b3f0e96 100644 --- a/TESTING/EIG/ddrvst.f +++ b/TESTING/EIG/ddrvst.f @@ -479,7 +479,7 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, DOUBLE PRECISION HALF PARAMETER ( HALF = 0.5D0 ) INTEGER MAXTYP - PARAMETER ( MAXTYP = 18 ) + PARAMETER ( MAXTYP = 19 ) * .. * .. Local Scalars .. LOGICAL BADNN @@ -519,11 +519,11 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT * .. * .. Data statements .. - DATA KTYPE / 1, 2, 5*4, 5*5, 3*8, 3*9 / + DATA KTYPE / 1, 2, 5*4, 5*5, 3*8, 3*9, 10 / DATA KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, - $ 2, 3, 1, 2, 3 / + $ 2, 3, 1, 2, 3, 2 / DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, - $ 0, 0, 4, 4, 4 / + $ 0, 0, 4, 4, 4, 0 / * .. * .. Executable Statements .. * @@ -741,6 +741,22 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, A( I, J ) = U( IROW, J ) 90 CONTINUE 100 CONTINUE +* + ELSE IF( ITYPE.EQ.10 ) THEN +* +* Tridiagonal with equal diagonal entries and half that +* off the diagonal, so that the absolute tolerance below +* is a meaningful fraction of the spread of the +* eigenvalues. +* + DO 105 JCOL = 1, N + A( JCOL, JCOL ) = ANORM + 105 CONTINUE + DO 106 JCOL = 1, N - 1 + A( JCOL, JCOL+1 ) = ANORM / TWO + A( JCOL+1, JCOL ) = A( JCOL, JCOL+1 ) + 106 CONTINUE +* ELSE IINFO = 1 END IF @@ -755,6 +771,13 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 110 CONTINUE * ABSTOL = UNFL + UNFL +* +* A tolerance at or below the smallest pivot is inert in the +* bisection, so type 19 asks for one that the routine has to +* scale along with the matrix. +* + IF( JTYPE.EQ.19 ) + $ ABSTOL = TEN*ULP*ANORM IF( N.LE.1 ) THEN IL = 1 IU = N @@ -770,7 +793,7 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * 3) If matrix is tridiagonal, call DSTEV and DSTEVX. * - IF( JTYPE.LE.7 ) THEN + IF( JTYPE.LE.7 .OR. JTYPE.EQ.19 ) THEN NTEST = 1 DO 120 I = 1, N D1( I ) = A( I, I ) @@ -2051,7 +2074,7 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * 6) Call DSBEV and DSBEVX. * - IF( JTYPE.LE.7 ) THEN + IF( JTYPE.LE.7 .OR. JTYPE.EQ.19 ) THEN KD = 1 ELSE IF( JTYPE.GE.8 .AND. JTYPE.LE.15 ) THEN KD = MAX( N-1, 0 ) @@ -2555,7 +2578,7 @@ SUBROUTINE DDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * 9) Call DSBEVD. * - IF( JTYPE.LE.7 ) THEN + IF( JTYPE.LE.7 .OR. JTYPE.EQ.19 ) THEN KD = 1 ELSE IF( JTYPE.GE.8 .AND. JTYPE.LE.15 ) THEN KD = MAX( N-1, 0 ) diff --git a/TESTING/EIG/schkee.F b/TESTING/EIG/schkee.F index 4fa913e30..5c6f0a9fb 100644 --- a/TESTING/EIG/schkee.F +++ b/TESTING/EIG/schkee.F @@ -1925,7 +1925,7 @@ PROGRAM SCHKEE $ A( 1, 3 ), D( 1, 12 ), A( 1, 4 ), WORK, $ LWORK, IWORK, LIWORK, RESULT, INFO ) ELSE - CALL SDRVST( NN, NVAL, 18, DOTYPE, ISEED, THRESH, + CALL SDRVST( NN, NVAL, 19, DOTYPE, ISEED, THRESH, $ NOUT, A( 1, 1 ), NMAX, D( 1, 3 ), D( 1, 4 ), $ D( 1, 5 ), D( 1, 6 ), D( 1, 8 ), D( 1, 9 ), $ D( 1, 10 ), D( 1, 11), A( 1, 2 ), NMAX, diff --git a/TESTING/EIG/sdrvst.f b/TESTING/EIG/sdrvst.f index 002f86a49..dc8eb6034 100644 --- a/TESTING/EIG/sdrvst.f +++ b/TESTING/EIG/sdrvst.f @@ -479,7 +479,7 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, REAL HALF PARAMETER ( HALF = 0.5E0 ) INTEGER MAXTYP - PARAMETER ( MAXTYP = 18 ) + PARAMETER ( MAXTYP = 19 ) * .. * .. Local Scalars .. LOGICAL BADNN @@ -519,11 +519,11 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, INTRINSIC ABS, INT, LOG, MAX, MIN, REAL, SQRT * .. * .. Data statements .. - DATA KTYPE / 1, 2, 5*4, 5*5, 3*8, 3*9 / + DATA KTYPE / 1, 2, 5*4, 5*5, 3*8, 3*9, 10 / DATA KMAGN / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1, - $ 2, 3, 1, 2, 3 / + $ 2, 3, 1, 2, 3, 2 / DATA KMODE / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0, - $ 0, 0, 4, 4, 4 / + $ 0, 0, 4, 4, 4, 0 / * .. * .. Executable Statements .. * @@ -741,6 +741,22 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, A( I, J ) = U( IROW, J ) 90 CONTINUE 100 CONTINUE +* + ELSE IF( ITYPE.EQ.10 ) THEN +* +* Tridiagonal with equal diagonal entries and half that +* off the diagonal, so that the absolute tolerance below +* is a meaningful fraction of the spread of the +* eigenvalues. +* + DO 105 JCOL = 1, N + A( JCOL, JCOL ) = ANORM + 105 CONTINUE + DO 106 JCOL = 1, N - 1 + A( JCOL, JCOL+1 ) = ANORM / TWO + A( JCOL+1, JCOL ) = A( JCOL, JCOL+1 ) + 106 CONTINUE +* ELSE IINFO = 1 END IF @@ -755,6 +771,13 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 110 CONTINUE * ABSTOL = UNFL + UNFL +* +* A tolerance at or below the smallest pivot is inert in the +* bisection, so type 19 asks for one that the routine has to +* scale along with the matrix. +* + IF( JTYPE.EQ.19 ) + $ ABSTOL = TEN*ULP*ANORM IF( N.LE.1 ) THEN IL = 1 IU = N @@ -770,7 +793,7 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * 3) If matrix is tridiagonal, call SSTEV and SSTEVX. * - IF( JTYPE.LE.7 ) THEN + IF( JTYPE.LE.7 .OR. JTYPE.EQ.19 ) THEN NTEST = 1 DO 120 I = 1, N D1( I ) = A( I, I ) @@ -2051,7 +2074,7 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * 6) Call SSBEV and SSBEVX. * - IF( JTYPE.LE.7 ) THEN + IF( JTYPE.LE.7 .OR. JTYPE.EQ.19 ) THEN KD = 1 ELSE IF( JTYPE.GE.8 .AND. JTYPE.LE.15 ) THEN KD = MAX( N-1, 0 ) @@ -2555,7 +2578,7 @@ SUBROUTINE SDRVST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * 9) Call SSBEVD. * - IF( JTYPE.LE.7 ) THEN + IF( JTYPE.LE.7 .OR. JTYPE.EQ.19 ) THEN KD = 1 ELSE IF( JTYPE.GE.8 .AND. JTYPE.LE.15 ) THEN KD = MAX( N-1, 0 )