From 0c58a6a05088a5a882a502bf1007775420653e35 Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Tue, 27 Oct 2026 08:59:39 +0100 Subject: [PATCH 1/4] Drop the absolute floor from the splitting test in ?STEBZ (Reference-LAPACK PR 1385) --- lapack-netlib/SRC/dlaebz.f | 2 +- lapack-netlib/SRC/dstebz.c | 66 ++++---- lapack-netlib/SRC/dstebz.f | 31 ++-- lapack-netlib/SRC/slaebz.f | 2 +- lapack-netlib/SRC/sstebz.c | 327 +++++-------------------------------- lapack-netlib/SRC/sstebz.f | 41 +++-- 6 files changed, 125 insertions(+), 344 deletions(-) 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 From de3bbd6ee119aacf6a584dc46d6ca210d87bbdaf Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Tue, 27 Oct 2026 09:01:57 +0100 Subject: [PATCH 2/4] Drop the absolute floor from the splitting test in ?STEBZ (Reference-LAPACK PR 1385) --- lapack-netlib/TESTING/sep.in | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From 99b04f0ff68c2320458f0d12c02e3de5787fd490 Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Tue, 27 Oct 2026 09:03:05 +0100 Subject: [PATCH 3/4] Drop the absolute floor from the splitting test in ?STEBZ (Reference-LAPACK PR 1385) --- lapack-netlib/TESTING/EIG/cchkee.F | 29 ++++++++++-- lapack-netlib/TESTING/EIG/cchkst.f | 71 ++++++++++++++++++++++-------- lapack-netlib/TESTING/EIG/dchkee.F | 68 +++++++++++++++++++++------- lapack-netlib/TESTING/EIG/dchkst.f | 46 ++++++++++++++++--- lapack-netlib/TESTING/EIG/schkee.F | 68 +++++++++++++++++++++------- lapack-netlib/TESTING/EIG/schkst.f | 46 ++++++++++++++++--- lapack-netlib/TESTING/EIG/zchkee.F | 29 ++++++++++-- lapack-netlib/TESTING/EIG/zchkst.f | 66 ++++++++++++++++++++------- 8 files changed, 337 insertions(+), 86 deletions(-) diff --git a/lapack-netlib/TESTING/EIG/cchkee.F b/lapack-netlib/TESTING/EIG/cchkee.F index a1f82400a1..dccc8ac5f3 100644 --- a/lapack-netlib/TESTING/EIG/cchkee.F +++ b/lapack-netlib/TESTING/EIG/cchkee.F @@ -1080,8 +1080,9 @@ PROGRAM CCHKEE REAL EPS, S1, S2, THRESH, THRSHN * .. * .. Local Arrays .. - LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ) - INTEGER IOLDSD( 4 ), ISEED( 4 ), IWORK( LIWORK ), + LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ), DOTYP1( 10 ) + INTEGER IOLDSD( 4 ), ISEED( 4 ), ISEED2( 4 ), + $ IWORK( LIWORK ), $ KVAL( MAXIN ), MVAL( MAXIN ), MXBVAL( MAXIN ), $ NBCOL( MAXIN ), NBMIN( MAXIN ), NBVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), @@ -1109,6 +1110,7 @@ PROGRAM CCHKEE $ CCHKGG, CCHKGK, CCHKGL, CCHKHS, CCHKST, CCKCSD, $ CCKGLM, CCKGQR, CCKGSV, CCKLSE, CDRGES, CDRGEV, $ CDRGSX, CDRGVX, CDRVBD, CDRVES, CDRVEV, CDRVST, + $ CDRVMR3, $ CDRVSX, CDRVVX, CERRBD, CERRED, CERRGG, CERRHS, $ CERRST, ILAVER, XLAENV, CDRGES3, CDRGEV3, $ CCHKST2STG, CDRVST2STG, CCHKHB2STG, CDRVSG2STG @@ -1136,6 +1138,7 @@ PROGRAM CCHKEE * .. * .. Data statements .. DATA INTSTR / '0123456789' / + DATA DOTYP1 / 10*.TRUE. / DATA IOLDSD / 0, 0, 0, 1 / * .. * .. Allocate memory dynamically .. @@ -1882,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 ) @@ -2053,6 +2056,26 @@ PROGRAM CCHKEE $ A( 1, 6 ), DR( 1, 1 ), DR( 1, 2 ), $ DR( 1, 3 ), WORK, LWORK, RWORK, IWORK, NOUT, $ INFO ) +* +* CDRVMR3 uses a seed of its own, taken from ISEED at the first +* parameter set, so that the matrices generated by CCHKBD and +* CDRVBD in the following parameter sets are the same as +* without it. +* + IF( TSTDRV ) THEN + IF( I.EQ.1 .OR. NEWSD.EQ.0 ) THEN + ISEED2( 1 ) = ISEED( 1 ) + ISEED2( 2 ) = ISEED( 2 ) + ISEED2( 3 ) = ISEED( 3 ) + ISEED2( 4 ) = ISEED( 4 ) + END IF + CALL CDRVMR3( NN, MVAL, NVAL, 10, DOTYP1, ISEED2, + $ THRESH, A( 1, 1 ), NMAX, A( 1, 2 ), NMAX, + $ A( 1, 3 ), NMAX, A( 1, 4 ), A( 1, 5 ), + $ A( 1, 6 ), DR( 1, 1 ), DR( 1, 2 ), + $ DR( 1, 3 ), WORK, LWORK, RWORK, IWORK, NOUT, + $ INFO ) + END IF 330 CONTINUE * ELSE IF( LSAMEN( 3, C3, 'CEV' ) ) THEN diff --git a/lapack-netlib/TESTING/EIG/cchkst.f b/lapack-netlib/TESTING/EIG/cchkst.f index ed535e91f3..a02a88c288 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, + $ XER_REPLACE * .. * .. 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 .. * @@ -719,7 +723,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XERBLA( 'CCHKST', -INFO ) + CALL XER_REPLACE( 'CCHKST', -INFO ) RETURN END IF * @@ -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 20c4cb47c2..31c604259c 100644 --- a/lapack-netlib/TESTING/EIG/dchkee.F +++ b/lapack-netlib/TESTING/EIG/dchkee.F @@ -775,7 +775,10 @@ *> printed about each test for which the test ratio is greater *> than or equal to the threshold. *> -*> lines 3-EOF: +*> line 3: TSTERR, LOGICAL +*> Flag indicating whether or not to test the error exits. +*> +*> lines 4-EOF: *> *> Input for testing the eigencondition routines consists of a set of *> specially constructed test cases and their solutions. The data @@ -1083,8 +1086,9 @@ PROGRAM DCHKEE DOUBLE PRECISION EPS, S1, S2, THRESH, THRSHN * .. * .. Local Arrays .. - LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ) - INTEGER IOLDSD( 4 ), ISEED( 4 ), IWORK( LIWORK ), + LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ), DOTYP1( 10 ) + INTEGER IOLDSD( 4 ), ISEED( 4 ), ISEED2( 4 ), + $ IWORK( LIWORK ), $ KVAL( MAXIN ), MVAL( MAXIN ), MXBVAL( MAXIN ), $ NBCOL( MAXIN ), NBMIN( MAXIN ), NBVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), @@ -1106,12 +1110,12 @@ PROGRAM DCHKEE * .. * .. External Subroutines .. EXTERNAL ALAREQ, DCHKBB, DCHKBD, DCHKBK, DCHKBL, DCHKEC, - $ DCHKGG, DCHKGK, DCHKGL, DCHKHS, DCHKSB, DCHKST, - $ DCKCSD, DCKGLM, DCKGQR, DCKGSV, DCKLSE, DDRGES, - $ DDRGEV, DDRGSX, DDRGVX, DDRVBD, DDRVES, DDRVEV, - $ DDRVSG, DDRVST, DDRVSX, DDRVVX, DERRBD, - $ DERRED, DERRGG, DERRHS, DERRST, ILAVER, XLAENV, - $ DDRGES3, DDRGEV3, + $ DCHKGG, DCHKGK, DCHKGL, DCHKHS, DCHKST, DCKCSD, + $ DCKGLM, DCKGQR, DCKGSV, DCKLSE, DDRGES, DDRGEV, + $ DDRGSX, DDRGVX, DDRVBD, DDRVES, DDRVEV, DDRVST, + $ DDRVMR3, + $ DDRVSX, DDRVVX, DERRBD, DERRED, DERRGG, DERRHS, + $ DERRST, ILAVER, XLAENV, DDRGES3, DDRGEV3, $ DCHKST2STG, DDRVST2STG, DCHKSB2STG, DDRVSG2STG * .. * .. Intrinsic Functions .. @@ -1137,18 +1141,31 @@ PROGRAM DCHKEE * .. * .. Data statements .. DATA INTSTR / '0123456789' / + DATA DOTYP1 / 10*.TRUE. / DATA IOLDSD / 0, 0, 0, 1 / * .. * .. Allocate memory dynamically .. * ALLOCATE ( A(NMAX*NMAX,NEED), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( B(NMAX*NMAX,5), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( C(NCMAX*NCMAX,NCMAX*NCMAX), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( WORK(LWORK), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF * .. * .. Executable Statements .. * @@ -1264,13 +1281,13 @@ PROGRAM DCHKEE * DEC: Eigencondition estimation * READ( NIN, FMT = * )THRESH + READ( NIN, FMT = * )TSTERR CALL XLAENV( 1, 1 ) CALL XLAENV( 12, 11 ) CALL XLAENV( 13, 2 ) CALL XLAENV( 14, 0 ) CALL XLAENV( 15, 2 ) CALL XLAENV( 16, 2 ) - TSTERR = .TRUE. CALL DCHKEC( THRESH, TSTERR, NIN, NOUT ) GO TO 10 ELSE @@ -1731,7 +1748,7 @@ PROGRAM DCHKEE * IF( FATAL ) THEN WRITE( NOUT, FMT = 9999 ) - STOP + STOP 1 END IF * * Read the input lines indicating the test path and its parameters. @@ -1867,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 ) @@ -2032,6 +2049,25 @@ PROGRAM DCHKEE $ A( 1, 3 ), NMAX, A( 1, 4 ), A( 1, 5 ), $ A( 1, 6 ), D( 1, 1 ), D( 1, 2 ), D( 1, 3 ), $ WORK, LWORK, IWORK, NOUT, INFO ) +* +* DDRVMR3 uses a seed of its own, taken from ISEED at the first +* parameter set, so that the matrices generated by DCHKBD and +* DDRVBD in the following parameter sets are the same as +* without it. +* + IF( TSTDRV ) THEN + IF( I.EQ.1 .OR. NEWSD.EQ.0 ) THEN + ISEED2( 1 ) = ISEED( 1 ) + ISEED2( 2 ) = ISEED( 2 ) + ISEED2( 3 ) = ISEED( 3 ) + ISEED2( 4 ) = ISEED( 4 ) + END IF + CALL DDRVMR3( NN, MVAL, NVAL, 10, DOTYP1, ISEED2, + $ THRESH, A( 1, 1 ), NMAX, A( 1, 2 ), NMAX, + $ A( 1, 3 ), NMAX, A( 1, 4 ), A( 1, 5 ), + $ A( 1, 6 ), D( 1, 1 ), D( 1, 2 ), D( 1, 3 ), + $ WORK, LWORK, IWORK, NOUT, INFO ) + END IF 330 CONTINUE * ELSE IF( LSAMEN( 3, C3, 'DEV' ) ) THEN diff --git a/lapack-netlib/TESTING/EIG/dchkst.f b/lapack-netlib/TESTING/EIG/dchkst.f index 95fa26c9d0..5b59f609c2 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 @@ -655,7 +657,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD, DSTEBZ, $ DSTEBR, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR, $ DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, - $ XERBLA, XLAENV + $ XER_REPLACE, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT @@ -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 .. * @@ -710,7 +712,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XERBLA( 'DCHKST', -INFO ) + CALL XER_REPLACE( 'DCHKST', -INFO ) RETURN END IF * @@ -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 fa780b695c..872f2d451e 100644 --- a/lapack-netlib/TESTING/EIG/schkee.F +++ b/lapack-netlib/TESTING/EIG/schkee.F @@ -775,7 +775,10 @@ *> printed about each test for which the test ratio is greater *> than or equal to the threshold. *> -*> lines 3-EOF: +*> line 3: TSTERR, LOGICAL +*> Flag indicating whether or not to test the error exits. +*> +*> lines 4-EOF: *> *> Input for testing the eigencondition routines consists of a set of *> specially constructed test cases and their solutions. The data @@ -1083,8 +1086,9 @@ PROGRAM SCHKEE REAL EPS, S1, S2, THRESH, THRSHN * .. * .. Local Arrays .. - LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ) - INTEGER IOLDSD( 4 ), ISEED( 4 ), IWORK( LIWORK ), + LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ), DOTYP1( 10 ) + INTEGER IOLDSD( 4 ), ISEED( 4 ), ISEED2( 4 ), + $ IWORK( LIWORK ), $ KVAL( MAXIN ), MVAL( MAXIN ), MXBVAL( MAXIN ), $ NBCOL( MAXIN ), NBMIN( MAXIN ), NBVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), @@ -1106,12 +1110,12 @@ PROGRAM SCHKEE * .. * .. External Subroutines .. EXTERNAL ALAREQ, SCHKBB, SCHKBD, SCHKBK, SCHKBL, SCHKEC, - $ SCHKGG, SCHKGK, SCHKGL, SCHKHS, SCHKSB, SCHKST, - $ SCKCSD, SCKGLM, SCKGQR, SCKGSV, SCKLSE, SDRGES, - $ SDRGEV, SDRGSX, SDRGVX, SDRVBD, SDRVES, SDRVEV, - $ SDRVSG, SDRVST, SDRVSX, SDRVVX, SERRBD, - $ SERRED, SERRGG, SERRHS, SERRST, ILAVER, XLAENV, - $ SDRGES3, SDRGEV3, + $ SCHKGG, SCHKGK, SCHKGL, SCHKHS, SCHKST, SCKCSD, + $ SCKGLM, SCKGQR, SCKGSV, SCKLSE, SDRGES, SDRGEV, + $ SDRGSX, SDRGVX, SDRVBD, SDRVES, SDRVEV, SDRVST, + $ SDRVMR3, + $ SDRVSX, SDRVVX, SERRBD, SERRED, SERRGG, SERRHS, + $ SERRST, ILAVER, XLAENV, SDRGES3, SDRGEV3, $ SCHKST2STG, SDRVST2STG, SCHKSB2STG, SDRVSG2STG * .. * .. Intrinsic Functions .. @@ -1137,18 +1141,31 @@ PROGRAM SCHKEE * .. * .. Data statements .. DATA INTSTR / '0123456789' / + DATA DOTYP1 / 10*.TRUE. / DATA IOLDSD / 0, 0, 0, 1 / * .. * .. Allocate memory dynamically .. * ALLOCATE ( A(NMAX*NMAX,NEED), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( B(NMAX*NMAX,5), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( C(NCMAX*NCMAX,NCMAX*NCMAX), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF ALLOCATE ( WORK(LWORK), STAT = AllocateStatus ) - IF (AllocateStatus /= 0) STOP "*** Not enough memory ***" + IF (AllocateStatus /= 0) THEN + WRITE( NOUT, FMT = * )'*** Not enough memory ***' + STOP 1 + END IF * .. * .. Executable Statements .. * @@ -1265,13 +1282,13 @@ PROGRAM SCHKEE * SEC: Eigencondition estimation * READ( NIN, FMT = * )THRESH + READ( NIN, FMT = * )TSTERR CALL XLAENV( 1, 1 ) CALL XLAENV( 12, 11 ) CALL XLAENV( 13, 2 ) CALL XLAENV( 14, 0 ) CALL XLAENV( 15, 2 ) CALL XLAENV( 16, 2 ) - TSTERR = .TRUE. CALL SCHKEC( THRESH, TSTERR, NIN, NOUT ) GO TO 10 ELSE @@ -1732,7 +1749,7 @@ PROGRAM SCHKEE * IF( FATAL ) THEN WRITE( NOUT, FMT = 9999 ) - STOP + STOP 1 END IF * * Read the input lines indicating the test path and its parameters. @@ -1868,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 ) @@ -2033,6 +2050,25 @@ PROGRAM SCHKEE $ A( 1, 3 ), NMAX, A( 1, 4 ), A( 1, 5 ), $ A( 1, 6 ), D( 1, 1 ), D( 1, 2 ), D( 1, 3 ), $ WORK, LWORK, IWORK, NOUT, INFO ) +* +* SDRVMR3 uses a seed of its own, taken from ISEED at the first +* parameter set, so that the matrices generated by SCHKBD and +* SDRVBD in the following parameter sets are the same as +* without it. +* + IF( TSTDRV ) THEN + IF( I.EQ.1 .OR. NEWSD.EQ.0 ) THEN + ISEED2( 1 ) = ISEED( 1 ) + ISEED2( 2 ) = ISEED( 2 ) + ISEED2( 3 ) = ISEED( 3 ) + ISEED2( 4 ) = ISEED( 4 ) + END IF + CALL SDRVMR3( NN, MVAL, NVAL, 10, DOTYP1, ISEED2, + $ THRESH, A( 1, 1 ), NMAX, A( 1, 2 ), NMAX, + $ A( 1, 3 ), NMAX, A( 1, 4 ), A( 1, 5 ), + $ A( 1, 6 ), D( 1, 1 ), D( 1, 2 ), D( 1, 3 ), + $ WORK, LWORK, IWORK, NOUT, INFO ) + END IF 330 CONTINUE * ELSE IF( LSAMEN( 3, C3, 'SEV' ) ) THEN diff --git a/lapack-netlib/TESTING/EIG/schkst.f b/lapack-netlib/TESTING/EIG/schkst.f index 7f251eb561..a21f307c69 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 @@ -655,7 +657,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ SOPGTR, SORGTR, SPTEQR, SSPT21, SSPTRD, SSTEBZ, $ SSTEBR, SSTECH, SSTEDC, SSTEMR, SSTEIN, SSTEQR, $ SSTERF, SSTT21, SSTT22, SSYT21, SSYTRD, - $ XERBLA, XLAENV + $ XER_REPLACE, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN, REAL, SQRT @@ -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 .. * @@ -710,7 +712,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XERBLA( 'SCHKST', -INFO ) + CALL XER_REPLACE( 'SCHKST', -INFO ) RETURN END IF * @@ -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 b00a412723..0491f821b2 100644 --- a/lapack-netlib/TESTING/EIG/zchkee.F +++ b/lapack-netlib/TESTING/EIG/zchkee.F @@ -1080,8 +1080,9 @@ PROGRAM ZCHKEE DOUBLE PRECISION EPS, S1, S2, THRESH, THRSHN * .. * .. Local Arrays .. - LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ) - INTEGER IOLDSD( 4 ), ISEED( 4 ), IWORK( LIWORK ), + LOGICAL DOTYPE( MAXT ), LOGWRK( NMAX ), DOTYP1( 10 ) + INTEGER IOLDSD( 4 ), ISEED( 4 ), ISEED2( 4 ), + $ IWORK( LIWORK ), $ KVAL( MAXIN ), MVAL( MAXIN ), MXBVAL( MAXIN ), $ NBCOL( MAXIN ), NBMIN( MAXIN ), NBVAL( MAXIN ), $ NSVAL( MAXIN ), NVAL( MAXIN ), NXVAL( MAXIN ), @@ -1109,6 +1110,7 @@ PROGRAM ZCHKEE $ ZCHKEC, ZCHKGG, ZCHKGK, ZCHKGL, ZCHKHS, ZCHKST, $ ZCKCSD, ZCKGLM, ZCKGQR, ZCKGSV, ZCKLSE, ZDRGES, $ ZDRGEV, ZDRGSX, ZDRGVX, ZDRVBD, ZDRVES, ZDRVEV, + $ ZDRVMR3, $ ZDRVST, ZDRVSX, ZDRVVX, ZERRBD, ZERRED, ZERRGG, $ ZERRHS, ZERRST, ILAVER, ZDRGES3, ZDRGEV3, $ ZCHKST2STG, ZDRVST2STG, ZCHKHB2STG, ZDRVSG2STG @@ -1136,6 +1138,7 @@ PROGRAM ZCHKEE * .. * .. Data statements .. DATA INTSTR / '0123456789' / + DATA DOTYP1 / 10*.TRUE. / DATA IOLDSD / 0, 0, 0, 1 / * .. * .. Allocate memory dynamically .. @@ -1882,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 ) @@ -2053,6 +2056,26 @@ PROGRAM ZCHKEE $ A( 1, 6 ), DR( 1, 1 ), DR( 1, 2 ), $ DR( 1, 3 ), WORK, LWORK, RWORK, IWORK, NOUT, $ INFO ) +* +* ZDRVMR3 uses a seed of its own, taken from ISEED at the first +* parameter set, so that the matrices generated by ZCHKBD and +* ZDRVBD in the following parameter sets are the same as +* without it. +* + IF( TSTDRV ) THEN + IF( I.EQ.1 .OR. NEWSD.EQ.0 ) THEN + ISEED2( 1 ) = ISEED( 1 ) + ISEED2( 2 ) = ISEED( 2 ) + ISEED2( 3 ) = ISEED( 3 ) + ISEED2( 4 ) = ISEED( 4 ) + END IF + CALL ZDRVMR3( NN, MVAL, NVAL, 10, DOTYP1, ISEED2, + $ THRESH, A( 1, 1 ), NMAX, A( 1, 2 ), NMAX, + $ A( 1, 3 ), NMAX, A( 1, 4 ), A( 1, 5 ), + $ A( 1, 6 ), DR( 1, 1 ), DR( 1, 2 ), + $ DR( 1, 3 ), WORK, LWORK, RWORK, IWORK, NOUT, + $ INFO ) + END IF 330 CONTINUE * ELSE IF( LSAMEN( 3, C3, 'ZEV' ) ) THEN diff --git a/lapack-netlib/TESTING/EIG/zchkst.f b/lapack-netlib/TESTING/EIG/zchkst.f index b77d944388..597e93729f 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, + $ XER_REPLACE, 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 .. * @@ -719,7 +723,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XERBLA( 'ZCHKST', -INFO ) + CALL XER_REPLACE( 'ZCHKST', -INFO ) RETURN END IF * @@ -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 From 67169c859c8ceee57a1c577c01e1298eda82638a Mon Sep 17 00:00:00 2001 From: Martin Kroeker Date: Wed, 28 Oct 2026 19:05:57 +0100 Subject: [PATCH 4/4] Rename XER_REPLACE back to XERBLA for now --- lapack-netlib/TESTING/EIG/cchkst.f | 4 ++-- lapack-netlib/TESTING/EIG/dchkst.f | 4 ++-- lapack-netlib/TESTING/EIG/schkst.f | 4 ++-- lapack-netlib/TESTING/EIG/zchkst.f | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/lapack-netlib/TESTING/EIG/cchkst.f b/lapack-netlib/TESTING/EIG/cchkst.f index a02a88c288..a681a8d111 100644 --- a/lapack-netlib/TESTING/EIG/cchkst.f +++ b/lapack-netlib/TESTING/EIG/cchkst.f @@ -669,7 +669,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ CLASET, CLATMR, CLATMS, CPTEQR, CSTEDC, CSTEMR, $ CSTEIN, CSTEQR, CSTT21, CSTT22, CUNGTR, CUPGTR, $ SCOPY, SLASUM, SSTEBZ, SSTECH, SSTERF, - $ XER_REPLACE + $ XERBLA * .. * .. Intrinsic Functions .. INTRINSIC ABS, CONJG, INT, LOG, MAX, MIN, REAL, SQRT @@ -723,7 +723,7 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XER_REPLACE( 'CCHKST', -INFO ) + CALL XERBLA( 'CCHKST', -INFO ) RETURN END IF * diff --git a/lapack-netlib/TESTING/EIG/dchkst.f b/lapack-netlib/TESTING/EIG/dchkst.f index 5b59f609c2..4df709e239 100644 --- a/lapack-netlib/TESTING/EIG/dchkst.f +++ b/lapack-netlib/TESTING/EIG/dchkst.f @@ -657,7 +657,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ DOPGTR, DORGTR, DPTEQR, DSPT21, DSPTRD, DSTEBZ, $ DSTEBR, DSTECH, DSTEDC, DSTEMR, DSTEIN, DSTEQR, $ DSTERF, DSTT21, DSTT22, DSYT21, DSYTRD, - $ XER_REPLACE, XLAENV + $ XERBLA, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, LOG, MAX, MIN, SQRT @@ -712,7 +712,7 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XER_REPLACE( 'DCHKST', -INFO ) + CALL XERBLA( 'DCHKST', -INFO ) RETURN END IF * diff --git a/lapack-netlib/TESTING/EIG/schkst.f b/lapack-netlib/TESTING/EIG/schkst.f index a21f307c69..da8ac3ed82 100644 --- a/lapack-netlib/TESTING/EIG/schkst.f +++ b/lapack-netlib/TESTING/EIG/schkst.f @@ -657,7 +657,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, $ SOPGTR, SORGTR, SPTEQR, SSPT21, SSPTRD, SSTEBZ, $ SSTEBR, SSTECH, SSTEDC, SSTEMR, SSTEIN, SSTEQR, $ SSTERF, SSTT21, SSTT22, SSYT21, SSYTRD, - $ XER_REPLACE, XLAENV + $ XERBLA, XLAENV * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN, REAL, SQRT @@ -712,7 +712,7 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XER_REPLACE( 'SCHKST', -INFO ) + CALL XERBLA( 'SCHKST', -INFO ) RETURN END IF * diff --git a/lapack-netlib/TESTING/EIG/zchkst.f b/lapack-netlib/TESTING/EIG/zchkst.f index 597e93729f..16c961215e 100644 --- a/lapack-netlib/TESTING/EIG/zchkst.f +++ b/lapack-netlib/TESTING/EIG/zchkst.f @@ -666,7 +666,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, * .. * .. External Subroutines .. EXTERNAL DCOPY, DLASUM, DSTEBZ, DSTECH, DSTERF, - $ XER_REPLACE, ZCOPY, ZHET21, ZHETRD, ZHPT21, + $ XERBLA, ZCOPY, ZHET21, ZHETRD, ZHPT21, $ ZHPTRD, ZLACPY, ZLASET, ZLATMR, ZLATMS, ZPTEQR, $ ZSTEDC, ZSTEMR, ZSTEIN, ZSTEQR, ZSTT21, ZSTT22, $ ZUNGTR, ZUPGTR @@ -723,7 +723,7 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, END IF * IF( INFO.NE.0 ) THEN - CALL XER_REPLACE( 'ZCHKST', -INFO ) + CALL XERBLA( 'ZCHKST', -INFO ) RETURN END IF *