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 c62aa2065..6ffe7bea7 100644 --- a/TESTING/EIG/dchkee.F +++ b/TESTING/EIG/dchkee.F @@ -1939,7 +1939,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 e9c560ee9..98e3bf65e 100644 --- a/TESTING/EIG/schkee.F +++ b/TESTING/EIG/schkee.F @@ -1940,7 +1940,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 )