diff --git a/lapack-netlib/SRC/dlaebz.f b/lapack-netlib/SRC/dlaebz.f index 980536ce67..c5021fd76e 100644 --- a/lapack-netlib/SRC/dlaebz.f +++ b/lapack-netlib/SRC/dlaebz.f @@ -63,7 +63,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/lapack-netlib/SRC/dstebz.c b/lapack-netlib/SRC/dstebz.c index 4668bcefcf..d35c33e895 100644 --- a/lapack-netlib/SRC/dstebz.c +++ b/lapack-netlib/SRC/dstebz.c @@ -226,8 +226,6 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; @@ -243,7 +241,6 @@ static integer c__0 = 0; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download DSTEBZ + dependencies */ /* > */ @@ -254,7 +251,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -502,9 +498,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date June 2016 */ - -/* > \ingroup auxOTHERcomputational */ +/* > \ingroup stebz */ /* ===================================================================== */ /* Subroutine */ void dstebz_(char *range, char *order, integer *n, doublereal @@ -518,7 +512,17 @@ f"> */ doublereal d__1, d__2, d__3, d__4, d__5; /* Local variables */ - integer iend, ioff, iout, itmp1, j, jdisc; + integer j, ib, jb, ie, je, nb; + doublereal gl; + integer im, in; + doublereal gu; + integer iw; + doublereal wl, wu; + integer nwl; + doublereal ulp, wlu, wul; + integer nwu; + doublereal tmp1, tmp2; + integer iend, ioff, iout, itmp1, jdisc; extern logical lsame_(char *, char *); integer iinfo; doublereal atoli; @@ -526,39 +530,28 @@ f"> */ doublereal bnorm; integer itmax; doublereal wkill, rtoli, tnorm; - integer ib, jb, ie, je, nb; - doublereal gl; - integer im, in; extern doublereal dlamch_(char *); integer ibegin; - doublereal gu; - integer iw; extern /* Subroutine */ void dlaebz_(integer *, integer *, integer *, integer *, integer *, integer *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, integer *, doublereal *, doublereal *, integer *, integer *, doublereal *, integer *, integer *); - doublereal wl; integer irange, idiscl; - doublereal safemn, wu; + doublereal safemn; integer idumma[1]; + extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); integer idiscu, iorder; logical ncnvrg; doublereal pivmin; logical toofew; - integer nwl; - doublereal ulp, wlu, wul; - integer nwu; - doublereal tmp1, tmp2; -/* -- LAPACK computational routine (version 3.7.0) -- */ +/* -- LAPACK computational routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* June 2016 */ /* ===================================================================== */ @@ -680,10 +673,17 @@ f"> */ /* Computing 2nd power */ d__1 = e[j - 1]; tmp1 = d__1 * d__1; -/* Computing 2nd power */ - d__2 = ulp; - if ((d__1 = d__[j] * d__[j - 1], abs(d__1)) * (d__2 * d__2) + safemn - > tmp1) { +/* 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 = (d__1 = d__[j], abs(d__1)) * ulp * ((d__2 = d__[j - 1], abs( + d__2)) * ulp); + if (tmp1 < tmp2 || tmp1 == 0.) { isplit[*nsplit] = j - 1; ++(*nsplit); work[j - 1] = 0.; @@ -841,8 +841,8 @@ f"> */ if (irange == 1 || wu >= d__[ibegin] - pivmin) { ++nwu; } - if (irange == 1 || (wl < d__[ibegin] - pivmin && wu >= d__[ibegin] - - pivmin)) { + if (irange == 1 || wl < d__[ibegin] - pivmin && wu >= d__[ibegin] + - pivmin) { ++(*m); w[*m] = d__[ibegin]; iblock[*m] = jb; @@ -894,7 +894,7 @@ f"> */ } if (irange > 1) { - if (gu < wl) { + if (gu <= wl) { nwl += in; nwu += in; goto L70; @@ -1007,7 +1007,9 @@ f"> */ } /* L90: */ } - iblock[iw] = 0; + if (iw != 0) { + iblock[iw] = 0; + } /* L100: */ } } @@ -1025,7 +1027,9 @@ f"> */ } /* L110: */ } - iblock[iw] = 0; + if (iw != 0) { + iblock[iw] = 0; + } /* L120: */ } } diff --git a/lapack-netlib/SRC/dstebz.f b/lapack-netlib/SRC/dstebz.f index baf417b0a5..d56f6a9df6 100644 --- a/lapack-netlib/SRC/dstebz.f +++ b/lapack-netlib/SRC/dstebz.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DSTEBZ + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -264,12 +262,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup auxOTHERcomputational +*> \ingroup stebz * * ===================================================================== - SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, + SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, + $ E, $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -419,7 +419,16 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, * 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 @@ -481,7 +490,8 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, IWORK( 5 ) = IL - 1 IWORK( 6 ) = IU * - CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, + CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, + $ E, $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, $ IWORK, W, IBLOCK, IINFO ) * @@ -595,7 +605,7 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, END IF * IF( IRANGE.GT.1 ) THEN - IF( GU.LT.WL ) THEN + IF( GU.LE.WL ) THEN NWL = NWL + IN NWU = NWU + IN GO TO 70 @@ -623,7 +633,8 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, * ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 - CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, + CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, + $ PIVMIN, $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) @@ -698,7 +709,7 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, WKILL = W( JE ) END IF 90 CONTINUE - IBLOCK( IW ) = 0 + IF( IW .NE. 0) IBLOCK( IW ) = 0 100 CONTINUE END IF IF( IDISCU.GT.0 ) THEN @@ -713,7 +724,7 @@ SUBROUTINE DSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, WKILL = W( JE ) END IF 110 CONTINUE - IBLOCK( IW ) = 0 + IF( IW .NE. 0) IBLOCK( IW ) = 0 120 CONTINUE END IF IM = 0 diff --git a/lapack-netlib/SRC/slaebz.f b/lapack-netlib/SRC/slaebz.f index 4eb8e8b7dc..233dbcb683 100644 --- a/lapack-netlib/SRC/slaebz.f +++ b/lapack-netlib/SRC/slaebz.f @@ -63,7 +63,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/lapack-netlib/SRC/sstebz.c b/lapack-netlib/SRC/sstebz.c index f6ca7ddc65..6d0005106a 100644 --- a/lapack-netlib/SRC/sstebz.c +++ b/lapack-netlib/SRC/sstebz.c @@ -220,258 +220,6 @@ typedef struct Namelist Namelist; //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} -/* procedure parameter types for -A and -C++ */ - - -#ifdef __cplusplus -typedef logical (*L_fp)(...); -#else -typedef logical (*L_fp)(); -#endif -#if 0 -static float spow_ui(float x, integer n) { - float pow=1.0; unsigned long int u; - if(n != 0) { - if(n < 0) n = -n, x = 1/x; - for(u = n; ; ) { - if(u & 01) pow *= x; - if(u >>= 1) x *= x; - else break; - } - } - return pow; -} -static double dpow_ui(double x, integer n) { - double pow=1.0; unsigned long int u; - if(n != 0) { - if(n < 0) n = -n, x = 1/x; - for(u = n; ; ) { - if(u & 01) pow *= x; - if(u >>= 1) x *= x; - else break; - } - } - return pow; -} -#ifdef _MSC_VER -static _Fcomplex cpow_ui(complex x, integer n) { - complex pow={1.0,0.0}; unsigned long int u; - if(n != 0) { - if(n < 0) n = -n, x.r = 1/x.r, x.i=1/x.i; - for(u = n; ; ) { - if(u & 01) pow.r *= x.r, pow.i *= x.i; - if(u >>= 1) x.r *= x.r, x.i *= x.i; - else break; - } - } - _Fcomplex p={pow.r, pow.i}; - return p; -} -#else -static _Complex float cpow_ui(_Complex float x, integer n) { - _Complex float pow=1.0; unsigned long int u; - if(n != 0) { - if(n < 0) n = -n, x = 1/x; - for(u = n; ; ) { - if(u & 01) pow *= x; - if(u >>= 1) x *= x; - else break; - } - } - return pow; -} -#endif -#ifdef _MSC_VER -static _Dcomplex zpow_ui(_Dcomplex x, integer n) { - _Dcomplex pow={1.0,0.0}; unsigned long int u; - if(n != 0) { - if(n < 0) n = -n, x._Val[0] = 1/x._Val[0], x._Val[1] =1/x._Val[1]; - for(u = n; ; ) { - if(u & 01) pow._Val[0] *= x._Val[0], pow._Val[1] *= x._Val[1]; - if(u >>= 1) x._Val[0] *= x._Val[0], x._Val[1] *= x._Val[1]; - else break; - } - } - _Dcomplex p = {pow._Val[0], pow._Val[1]}; - return p; -} -#else -static _Complex double zpow_ui(_Complex double x, integer n) { - _Complex double pow=1.0; unsigned long int u; - if(n != 0) { - if(n < 0) n = -n, x = 1/x; - for(u = n; ; ) { - if(u & 01) pow *= x; - if(u >>= 1) x *= x; - else break; - } - } - return pow; -} -#endif -static integer pow_ii(integer x, integer n) { - integer pow; unsigned long int u; - if (n <= 0) { - if (n == 0 || x == 1) pow = 1; - else if (x != -1) pow = x == 0 ? 1/x : 0; - else n = -n; - } - if ((n > 0) || !(n == 0 || x == 1 || x != -1)) { - u = n; - for(pow = 1; ; ) { - if(u & 01) pow *= x; - if(u >>= 1) x *= x; - else break; - } - } - return pow; -} -static integer dmaxloc_(double *w, integer s, integer e, integer *n) -{ - double m; integer i, mi; - for(m=w[s-1], mi=s, i=s+1; i<=e; i++) - if (w[i-1]>m) mi=i ,m=w[i-1]; - return mi-s+1; -} -static integer smaxloc_(float *w, integer s, integer e, integer *n) -{ - float m; integer i, mi; - for(m=w[s-1], mi=s, i=s+1; i<=e; i++) - if (w[i-1]>m) mi=i ,m=w[i-1]; - return mi-s+1; -} -static inline void cdotc_(complex *z, integer *n_, complex *x, integer *incx_, complex *y, integer *incy_) { - integer n = *n_, incx = *incx_, incy = *incy_, i; -#ifdef _MSC_VER - _Fcomplex zdotc = {0.0, 0.0}; - if (incx == 1 && incy == 1) { - for (i=0;i \htmlonly */ /* > Download SSTEBZ + dependencies */ /* > */ @@ -499,7 +246,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -747,9 +493,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date June 2016 */ - -/* > \ingroup auxOTHERcomputational */ +/* > \ingroup stebz */ /* ===================================================================== */ /* Subroutine */ void sstebz_(char *range, char *order, integer *n, real *vl, @@ -762,7 +506,17 @@ f"> */ real r__1, r__2, r__3, r__4, r__5; /* Local variables */ - integer iend, ioff, iout, itmp1, j, jdisc; + integer j, ib, jb, ie, je, nb; + real gl; + integer im, in; + real gu; + integer iw; + real wl, wu; + integer nwl; + real ulp, wlu, wul; + integer nwu; + real tmp1, tmp2; + integer iend, ioff, iout, itmp1, jdisc; extern logical lsame_(char *, char *); integer iinfo; real atoli; @@ -770,19 +524,13 @@ f"> */ real bnorm; integer itmax; real wkill, rtoli, tnorm; - integer ib, jb, ie, je, nb; - real gl; - integer im, in, ibegin; - real gu; - integer iw; - real wl; - integer irange, idiscl; + integer ibegin, irange, idiscl; extern real slamch_(char *); - real safemn, wu; + real safemn; integer idumma[1]; + extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); integer idiscu; extern /* Subroutine */ void slaebz_(integer *, integer *, integer *, integer *, integer *, integer *, real *, real *, real *, real *, @@ -792,16 +540,11 @@ f"> */ logical ncnvrg; real pivmin; logical toofew; - integer nwl; - real ulp, wlu, wul; - integer nwu; - real tmp1, tmp2; -/* -- LAPACK computational routine (version 3.7.0) -- */ +/* -- LAPACK computational routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* June 2016 */ /* ===================================================================== */ @@ -923,10 +666,17 @@ f"> */ /* Computing 2nd power */ r__1 = e[j - 1]; tmp1 = r__1 * r__1; -/* Computing 2nd power */ - r__2 = ulp; - if ((r__1 = d__[j] * d__[j - 1], abs(r__1)) * (r__2 * r__2) + safemn - > tmp1) { +/* 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 = (r__1 = d__[j], abs(r__1)) * ulp * ((r__2 = d__[j - 1], abs( + r__2)) * ulp); + if (tmp1 < tmp2 || tmp1 == 0.f) { isplit[*nsplit] = j - 1; ++(*nsplit); work[j - 1] = 0.f; @@ -975,8 +725,9 @@ f"> */ /* Computing MAX */ r__1 = abs(gl), r__2 = abs(gu); tnorm = f2cmax(r__1,r__2); - gl = gl - tnorm * 2.1f * ulp * *n - pivmin * 4.2000000000000002f; - gu = gu + tnorm * 2.1f * ulp * *n + pivmin * 2.1f; + gl = gl - tnorm * 2.3f * ulp * (real) (*n) - pivmin * + 4.5999999999999996f; + gu = gu + tnorm * 2.3f * ulp * (real) (*n) + pivmin * 2.3f; /* Compute Iteration parameters */ @@ -1085,8 +836,8 @@ f"> */ if (irange == 1 || wu >= d__[ibegin] - pivmin) { ++nwu; } - if (irange == 1 || (wl < d__[ibegin] - pivmin && wu >= d__[ibegin] - - pivmin)) { + if (irange == 1 || wl < d__[ibegin] - pivmin && wu >= d__[ibegin] + - pivmin) { ++(*m); w[*m] = d__[ibegin]; iblock[*m] = jb; @@ -1124,8 +875,8 @@ f"> */ /* Computing MAX */ r__1 = abs(gl), r__2 = abs(gu); bnorm = f2cmax(r__1,r__2); - gl = gl - bnorm * 2.1f * ulp * in - pivmin * 2.1f; - gu = gu + bnorm * 2.1f * ulp * in + pivmin * 2.1f; + gl = gl - bnorm * 2.3f * ulp * (real) in - pivmin * 2.3f; + gu = gu + bnorm * 2.3f * ulp * (real) in + pivmin * 2.3f; /* Compute ATOLI for the current submatrix */ @@ -1138,7 +889,7 @@ f"> */ } if (irange > 1) { - if (gu < wl) { + if (gu <= wl) { nwl += in; nwu += in; goto L70; @@ -1251,7 +1002,9 @@ f"> */ } /* L90: */ } - iblock[iw] = 0; + if (iw != 0) { + iblock[iw] = 0; + } /* L100: */ } } @@ -1269,7 +1022,9 @@ f"> */ } /* L110: */ } - iblock[iw] = 0; + if (iw != 0) { + iblock[iw] = 0; + } /* L120: */ } } diff --git a/lapack-netlib/SRC/sstebz.f b/lapack-netlib/SRC/sstebz.f index 13e7cde43c..ef7e78ee6a 100644 --- a/lapack-netlib/SRC/sstebz.f +++ b/lapack-netlib/SRC/sstebz.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SSTEBZ + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -264,12 +262,14 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup auxOTHERcomputational +*> \ingroup stebz * * ===================================================================== - SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, + SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, + $ E, $ M, NSPLIT, W, IBLOCK, ISPLIT, WORK, IWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK computational routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -292,7 +292,7 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0, $ HALF = 1.0E0 / TWO ) REAL FUDGE, RELFAC - PARAMETER ( FUDGE = 2.1E0, RELFAC = 2.0E0 ) + PARAMETER ( FUDGE = 2.3E0, RELFAC = 2.0E0 ) * .. * .. Local Scalars .. LOGICAL NCNVRG, TOOFEW @@ -418,7 +418,16 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, * 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 @@ -454,8 +463,8 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, GU = MAX( GU, D( N )+TMP1 ) GL = MIN( GL, D( N )-TMP1 ) TNORM = MAX( ABS( GL ), ABS( GU ) ) - GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN - GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN + GL = GL - FUDGE*TNORM*ULP*REAL( N ) - FUDGE*TWO*PIVMIN + GU = GU + FUDGE*TNORM*ULP*REAL( N ) + FUDGE*PIVMIN * * Compute Iteration parameters * @@ -480,7 +489,8 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, IWORK( 5 ) = IL - 1 IWORK( 6 ) = IU * - CALL SLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E, + CALL SLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, + $ E, $ WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, $ IWORK, W, IBLOCK, IINFO ) * @@ -582,8 +592,8 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, GU = MAX( GU, D( IEND )+TMP1 ) GL = MIN( GL, D( IEND )-TMP1 ) BNORM = MAX( ABS( GL ), ABS( GU ) ) - GL = GL - FUDGE*BNORM*ULP*IN - FUDGE*PIVMIN - GU = GU + FUDGE*BNORM*ULP*IN + FUDGE*PIVMIN + GL = GL - FUDGE*BNORM*ULP*REAL( IN ) - FUDGE*PIVMIN + GU = GU + FUDGE*BNORM*ULP*REAL( IN ) + FUDGE*PIVMIN * * Compute ATOLI for the current submatrix * @@ -594,7 +604,7 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, END IF * IF( IRANGE.GT.1 ) THEN - IF( GU.LT.WL ) THEN + IF( GU.LE.WL ) THEN NWL = NWL + IN NWU = NWU + IN GO TO 70 @@ -622,7 +632,8 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, * ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 - CALL SLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, + CALL SLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, + $ PIVMIN, $ D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) @@ -697,7 +708,7 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, WKILL = W( JE ) END IF 90 CONTINUE - IBLOCK( IW ) = 0 + IF( IW .NE. 0) IBLOCK( IW ) = 0 100 CONTINUE END IF IF( IDISCU.GT.0 ) THEN @@ -712,7 +723,7 @@ SUBROUTINE SSTEBZ( RANGE, ORDER, N, VL, VU, IL, IU, ABSTOL, D, E, WKILL = W( JE ) END IF 110 CONTINUE - IBLOCK( IW ) = 0 + IF( IW .NE. 0) IBLOCK( IW ) = 0 120 CONTINUE END IF IM = 0 diff --git a/lapack-netlib/TESTING/EIG/cchkee.F b/lapack-netlib/TESTING/EIG/cchkee.F index eba8022ad4..dccc8ac5f3 100644 --- a/lapack-netlib/TESTING/EIG/cchkee.F +++ b/lapack-netlib/TESTING/EIG/cchkee.F @@ -1885,7 +1885,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/lapack-netlib/TESTING/EIG/cchkst.f b/lapack-netlib/TESTING/EIG/cchkst.f index ed535e91f3..a681a8d111 100644 --- a/lapack-netlib/TESTING/EIG/cchkst.f +++ b/lapack-netlib/TESTING/EIG/cchkst.f @@ -601,6 +601,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, $ LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, $ INFO ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -632,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 @@ -665,18 +668,19 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, EXTERNAL CCOPY, CHET21, CHETRD, CHPT21, CHPTRD, CLACPY, $ CLASET, CLATMR, CLATMS, CPTEQR, CSTEDC, CSTEMR, $ CSTEIN, CSTEQR, CSTT21, CSTT22, CUNGTR, CUPGTR, - $ SCOPY, SLASUM, SSTEBZ, SSTECH, SSTERF, XERBLA + $ SCOPY, SLASUM, SSTEBZ, SSTECH, SSTERF, + $ XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, CONJG, INT, LOG, MAX, MIN, REAL, SQRT * .. * .. 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 .. * @@ -805,7 +809,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 @@ -816,7 +820,16 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, GO TO 70 * 60 CONTINUE - ANORM = RTUNFL*N*ULPINV + 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 @@ -899,6 +912,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 * @@ -1300,8 +1327,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * Do test 17 * - TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) / - $ ( ONE-HALF )**4 + TEMP2 = TWO*( TWO*REAL( N )-ONE )*ULP* + $ ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4 * TEMP1 = ZERO DO 190 J = 1, N @@ -1352,8 +1379,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, IL = 1 IU = N ELSE - IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) - IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) + IL = 1 + INT( REAL( N )*SLARND( 1, ISEED2 ) ) + IU = 1 + INT( REAL( N )*SLARND( 1, ISEED2 ) ) IF( IU.LT.IL ) THEN ITEMP = IU IU = IL @@ -1431,6 +1458,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.) @@ -1604,7 +1638,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * Do test 27 * - TEMP2 = TWO*( TWO*N-ONE )*ULP*( ONE+EIGHT*HALF**2 ) / + TEMP2 = TWO*( TWO*REAL( N )-ONE )*ULP* + $ ( ONE+EIGHT*HALF**2 ) / $ ( ONE-HALF )**4 * TEMP1 = ZERO @@ -1615,8 +1650,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * RESULT( 27 ) = TEMP1 / TEMP2 * - IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) - IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) + IL = 1 + INT( REAL( N )*SLARND( 1, ISEED2 ) ) + IU = 1 + INT( REAL( N )*SLARND( 1, ISEED2 ) ) IF( IU.LT.IL ) THEN ITEMP = IU IU = IL @@ -1646,7 +1681,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * * Do test 28 * - TEMP2 = TWO*( TWO*N-ONE )*ULP* + TEMP2 = TWO*( TWO*REAL( N )-ONE )*ULP* $ ( ONE+EIGHT*HALF**2 ) / ( ONE-HALF )**4 * TEMP1 = ZERO @@ -1675,8 +1710,8 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * IF( CRANGE ) THEN NTEST = 29 - IL = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) - IU = 1 + ( N-1 )*INT( SLARND( 1, ISEED2 ) ) + IL = 1 + INT( REAL( N )*SLARND( 1, ISEED2 ) ) + IU = 1 + INT( REAL( N )*SLARND( 1, ISEED2 ) ) IF( IU.LT.IL ) THEN ITEMP = IU IU = IL diff --git a/lapack-netlib/TESTING/EIG/dchkee.F b/lapack-netlib/TESTING/EIG/dchkee.F index 772f05a74b..31c604259c 100644 --- a/lapack-netlib/TESTING/EIG/dchkee.F +++ b/lapack-netlib/TESTING/EIG/dchkee.F @@ -1884,7 +1884,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/lapack-netlib/TESTING/EIG/dchkst.f b/lapack-netlib/TESTING/EIG/dchkst.f index 95fa26c9d0..4df709e239 100644 --- a/lapack-netlib/TESTING/EIG/dchkst.f +++ b/lapack-netlib/TESTING/EIG/dchkst.f @@ -621,8 +621,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 @@ -663,11 +665,11 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. Data statements .. DATA SMLSZS / 2, 8, 25 / 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 .. * @@ -794,7 +796,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 @@ -807,6 +809,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 * @@ -892,6 +903,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 * @@ -1423,6 +1448,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/lapack-netlib/TESTING/EIG/schkee.F b/lapack-netlib/TESTING/EIG/schkee.F index 3a79ccca0c..872f2d451e 100644 --- a/lapack-netlib/TESTING/EIG/schkee.F +++ b/lapack-netlib/TESTING/EIG/schkee.F @@ -1885,7 +1885,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/lapack-netlib/TESTING/EIG/schkst.f b/lapack-netlib/TESTING/EIG/schkst.f index 7f251eb561..da8ac3ed82 100644 --- a/lapack-netlib/TESTING/EIG/schkst.f +++ b/lapack-netlib/TESTING/EIG/schkst.f @@ -621,8 +621,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 @@ -663,11 +665,11 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. Data statements .. DATA SMLSZS / 2, 8, 25 / 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 .. * @@ -794,7 +796,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 @@ -807,6 +809,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 * @@ -892,6 +903,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 * @@ -1423,6 +1448,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/lapack-netlib/TESTING/EIG/zchkee.F b/lapack-netlib/TESTING/EIG/zchkee.F index dc90708f5e..0491f821b2 100644 --- a/lapack-netlib/TESTING/EIG/zchkee.F +++ b/lapack-netlib/TESTING/EIG/zchkee.F @@ -1885,7 +1885,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/lapack-netlib/TESTING/EIG/zchkst.f b/lapack-netlib/TESTING/EIG/zchkst.f index b77d944388..16c961215e 100644 --- a/lapack-netlib/TESTING/EIG/zchkst.f +++ b/lapack-netlib/TESTING/EIG/zchkst.f @@ -601,6 +601,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ WA1, WA2, WA3, WR, U, LDU, V, VP, TAU, Z, WORK, $ LWORK, RWORK, LRWORK, IWORK, LIWORK, RESULT, $ INFO ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -632,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 @@ -662,21 +665,22 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, EXTERNAL ILAENV, DLAMCH, DLARND, DSXT1 * .. * .. External Subroutines .. - EXTERNAL DCOPY, DLASUM, DSTEBZ, DSTECH, DSTERF, XERBLA, - $ ZCOPY, ZHET21, ZHETRD, ZHPT21, ZHPTRD, ZLACPY, - $ ZLASET, ZLATMR, ZLATMS, ZPTEQR, ZSTEDC, ZSTEMR, - $ ZSTEIN, ZSTEQR, ZSTT21, ZSTT22, ZUNGTR, ZUPGTR + EXTERNAL DCOPY, DLASUM, DSTEBZ, DSTECH, DSTERF, + $ XERBLA, ZCOPY, ZHET21, ZHETRD, ZHPT21, + $ ZHPTRD, ZLACPY, ZLASET, ZLATMR, ZLATMS, ZPTEQR, + $ ZSTEDC, ZSTEMR, ZSTEIN, ZSTEQR, ZSTT21, ZSTT22, + $ ZUNGTR, ZUPGTR * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCONJG, INT, LOG, MAX, MIN, SQRT * .. * .. 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 .. * @@ -805,7 +809,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 @@ -818,6 +822,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 * @@ -899,6 +912,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 * @@ -1352,8 +1379,8 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, IL = 1 IU = N ELSE - IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) - IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) + IL = 1 + INT( N*DLARND( 1, ISEED2 ) ) + IU = 1 + INT( N*DLARND( 1, ISEED2 ) ) IF( IU.LT.IL ) THEN ITEMP = IU IU = IL @@ -1431,6 +1458,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.) @@ -1615,8 +1649,8 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * RESULT( 27 ) = TEMP1 / TEMP2 * - IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) - IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) + IL = 1 + INT( N*DLARND( 1, ISEED2 ) ) + IU = 1 + INT( N*DLARND( 1, ISEED2 ) ) IF( IU.LT.IL ) THEN ITEMP = IU IU = IL @@ -1675,8 +1709,8 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * IF( CRANGE ) THEN NTEST = 29 - IL = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) - IU = 1 + ( N-1 )*INT( DLARND( 1, ISEED2 ) ) + IL = 1 + INT( N*DLARND( 1, ISEED2 ) ) + IU = 1 + INT( N*DLARND( 1, ISEED2 ) ) IF( IU.LT.IL ) THEN ITEMP = IU IU = IL diff --git a/lapack-netlib/TESTING/sep.in b/lapack-netlib/TESTING/sep.in index 99b47c0aa1..60442870f9 100644 --- a/lapack-netlib/TESTING/sep.in +++ b/lapack-netlib/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