From fcc595819e08dab21c47302990fa9b2cc02f2112 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Mon, 7 Sep 2026 20:27:42 -0700 Subject: [PATCH] Bound the bracket widening in xLARRB and xLARRJ so that a NaN matrix cannot hang xSTEMR Before bisecting, xLARRB widens the initial interval [W-WERR, W+WERR] until the Sturm count of the representation agrees with the eigenvalue index, moving the endpoint by BACK and doubling BACK at every step. In exact arithmetic that ends quickly, because the count is 0 below the spectrum and N above it. A NaN pivot is never counted, so with a NaN in the representation the count can stay short of the index and the loop runs forever: xSTEMR, and through it xSTEVR, xSYEVR and xHEEVR with RANGE = 'A', never return for a symmetric tridiagonal matrix with a NaN in one of its first two rows. The loop also stalls when WERR is exactly zero, which xLARRE produces for an eigenvalue that dqds returns as exactly zero. xLARRJ, which refines the eigenvalues by bisection on the matrix itself when only eigenvalues are wanted, has the same two widening loops with the same step and the same hang: xSTEMR with JOBZ = 'N' never returns for the same matrices. Start BACK at no less than the minimum interval width 2*PIVMIN, and stop widening with INFO = 1 once BACK has overflowed, when no further step can move the endpoint; a finite matrix never gets there, since BACK passes the spectral diameter within MAXITR doublings. The xLARRB callers already propagate a nonzero INFO (xLARRV returns -1, xLARRE -4, xSTEMR 21 or 11) and the drivers then fall back to xSTEBZ and xSTEIN; the one xLARRB call in xLARRV that ignored INFO now checks it like the others, in all four precisions. xSTEMR never looked at the INFO of xLARRJ; it now returns INFO = 3X for a failure there, next to the documented 1X and 2X codes of xLARRE and xLARRV. xCHKST gets the case as a regression test: after the size and type loops it calls xSTEMR with RANGE = 'A' on a tridiagonal matrix whose first diagonal entry is a NaN, once with JOBZ = 'V' and once with JOBZ = 'N', so that both refinement paths are covered. Only a negative INFO is reported as a failure; the test is that the call returns at all. On the parent all four xeigtst binaries stop responding there. Over 100 NaN cases (DSTEMR, SSTEMR, DSYEVR, DSTEVR, ZHEEVR; five NaN positions; n = 3, 5, 8, 20) the parent hangs in 60, this branch returns in all of them; with JOBZ = 'N' (xSTEMR in all four precisions, five NaN positions, n = 3, 5, 8, 20) the parent hangs in 48 of 80, this branch returns INFO = 31 in every case that hung. 200 finite matrices give bit-identical eigenvalues from DSTEMR on both, with JOBZ = 'V' and with JOBZ = 'N'. The full LAPACK test suite passes: 5441941 LAPACK tests, 0 numerical errors, 0 other errors, 40 tests more than the parent from the new calls. Co-Authored-By: Claude Fable 5.1 --- SRC/clarrv.f | 4 ++++ SRC/cstemr.f | 11 +++++++--- SRC/dlarrb.f | 38 +++++++++++++++++++++++++-------- SRC/dlarrj.f | 37 +++++++++++++++++++++++++------- SRC/dlarrv.f | 4 ++++ SRC/dstemr.f | 11 +++++++--- SRC/slarrb.f | 38 +++++++++++++++++++++++++-------- SRC/slarrj.f | 37 +++++++++++++++++++++++++------- SRC/slarrv.f | 4 ++++ SRC/sstemr.f | 11 +++++++--- SRC/zlarrv.f | 4 ++++ SRC/zstemr.f | 11 +++++++--- TESTING/EIG/cchkst.f | 51 ++++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/dchkst.f | 50 +++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/schkst.f | 50 +++++++++++++++++++++++++++++++++++++++++++ TESTING/EIG/zchkst.f | 51 ++++++++++++++++++++++++++++++++++++++++++++ 16 files changed, 366 insertions(+), 46 deletions(-) diff --git a/SRC/clarrv.f b/SRC/clarrv.f index dac574d75..ed9ae14fb 100644 --- a/SRC/clarrv.f +++ b/SRC/clarrv.f @@ -702,6 +702,10 @@ SUBROUTINE CLARRV( N, VL, VU, D, L, PIVMIN, $ WERR(WBEGIN),WORK( INDWRK ), $ IWORK( IINDWK ), PIVMIN, SPDIAM, $ IN, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = -1 + RETURN + END IF 55 CONTINUE * IF((WBEGIN+NEWLST-1.LT.DOL).OR. diff --git a/SRC/cstemr.f b/SRC/cstemr.f index 15dc809d3..8ff54a2d4 100644 --- a/SRC/cstemr.f +++ b/SRC/cstemr.f @@ -304,10 +304,11 @@ *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = 1X, internal error in SLARRE, -*> if INFO = 2X, internal error in CLARRV. +*> if INFO = 2X, internal error in CLARRV, +*> if INFO = 3X, internal error in SLARRJ. *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is -*> the nonzero error code returned by SLARRE or -*> CLARRV, respectively. +*> the nonzero error code returned by SLARRE, CLARRV +*> or SLARRJ, respectively. *> \endverbatim * * Authors: @@ -752,6 +753,10 @@ SUBROUTINE CSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, $ WORK( INDERR+WBEGIN-1 ), $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, $ TNRM, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = 30 + ABS( IINFO ) + RETURN + END IF IBEGIN = IEND + 1 WBEGIN = WEND + 1 39 CONTINUE diff --git a/SRC/dlarrb.f b/SRC/dlarrb.f index 38555f818..c3a923aa6 100644 --- a/SRC/dlarrb.f +++ b/SRC/dlarrb.f @@ -165,7 +165,11 @@ *> \param[out] INFO *> \verbatim *> INFO is INTEGER -*> Error flag. +*> = 0: successful exit +*> = 1: the interval around one of the eigenvalues could not +*> be widened to contain it: the Sturm counts of the +*> representation are inconsistent, as they are when it +*> contains a NaN or an Inf. *> \endverbatim * * Authors: @@ -274,10 +278,23 @@ SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1, * * Do while( NEGCNT(LEFT).GT.I-1 ) * - BACK = WERR( II ) +* The interval is widened by BACK, doubled at every step, until +* the Sturm count agrees with the index. In exact arithmetic +* this terminates because the count is 0 below the spectrum and +* N above it. A NaN or an Inf in the representation is never +* counted, so the count can stay short of the index; stop once +* BACK has overflowed, when no further step can move the +* endpoint. BACK starts at no less than the minimum interval +* width so that a zero WERR still makes progress. +* + BACK = MAX( WERR( II ), MNWDTH ) 20 CONTINUE NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R ) IF( NEGCNT.GT.I-1 ) THEN + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF LEFT = LEFT - BACK BACK = TWO*BACK GO TO 20 @@ -286,15 +303,18 @@ SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1, * Do while( NEGCNT(RIGHT).LT.I ) * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT * - BACK = WERR( II ) + BACK = MAX( WERR( II ), MNWDTH ) 50 CONTINUE - NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R ) - IF( NEGCNT.LT.I ) THEN - RIGHT = RIGHT + BACK - BACK = TWO*BACK - GO TO 50 - END IF + IF( NEGCNT.LT.I ) THEN + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF + RIGHT = RIGHT + BACK + BACK = TWO*BACK + GO TO 50 + END IF WIDTH = HALF*ABS( LEFT - RIGHT ) TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) CVRGD = MAX(RTOL1*GAP,RTOL2*TMP) diff --git a/SRC/dlarrj.f b/SRC/dlarrj.f index af3110f6e..5a8a43a81 100644 --- a/SRC/dlarrj.f +++ b/SRC/dlarrj.f @@ -137,7 +137,11 @@ *> \param[out] INFO *> \verbatim *> INFO is INTEGER -*> Error flag. +*> = 0: successful exit +*> = 1: the interval around one of the eigenvalues could not +*> be widened to contain it: the Sturm counts of the +*> matrix are inconsistent, as they are when it contains +*> a NaN or an Inf. *> \endverbatim * * Authors: @@ -190,7 +194,7 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, * .. Local Scalars .. INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT, $ OLNINT, P, PREV, SAVI1 - DOUBLE PRECISION DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH + DOUBLE PRECISION BACK, DPLUS, LEFT, MID, RIGHT, S, TMP, WIDTH * * .. * .. Intrinsic Functions .. @@ -249,7 +253,16 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, * * Do while( CNT(LEFT).GT.I-1 ) * - FAC = ONE +* The interval is widened by BACK, doubled at every step, +* until the Sturm count agrees with the index. In exact +* arithmetic this terminates because the count is 0 below the +* spectrum and N above it. A NaN or an Inf in the matrix is +* never counted, so the count can stay short of the index; +* stop once BACK has overflowed, when no further step can move +* the endpoint. BACK starts at no less than the minimum +* interval width so that a zero WERR still makes progress. +* + BACK = MAX( WERR( II ), TWO*PIVMIN ) 20 CONTINUE CNT = 0 S = LEFT @@ -260,14 +273,18 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, IF( DPLUS.LT.ZERO ) CNT = CNT + 1 30 CONTINUE IF( CNT.GT.I-1 ) THEN - LEFT = LEFT - WERR( II )*FAC - FAC = TWO*FAC + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF + LEFT = LEFT - BACK + BACK = TWO*BACK GO TO 20 END IF * * Do while( CNT(RIGHT).LT.I ) * - FAC = ONE + BACK = MAX( WERR( II ), TWO*PIVMIN ) 50 CONTINUE CNT = 0 S = RIGHT @@ -278,8 +295,12 @@ SUBROUTINE DLARRJ( N, D, E2, IFIRST, ILAST, IF( DPLUS.LT.ZERO ) CNT = CNT + 1 60 CONTINUE IF( CNT.LT.I ) THEN - RIGHT = RIGHT + WERR( II )*FAC - FAC = TWO*FAC + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF + RIGHT = RIGHT + BACK + BACK = TWO*BACK GO TO 50 END IF NINT = NINT + 1 diff --git a/SRC/dlarrv.f b/SRC/dlarrv.f index 4cf354ab7..19fe245d4 100644 --- a/SRC/dlarrv.f +++ b/SRC/dlarrv.f @@ -698,6 +698,10 @@ SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN, $ WERR(WBEGIN),WORK( INDWRK ), $ IWORK( IINDWK ), PIVMIN, SPDIAM, $ IN, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = -1 + RETURN + END IF 55 CONTINUE * IF((WBEGIN+NEWLST-1.LT.DOL).OR. diff --git a/SRC/dstemr.f b/SRC/dstemr.f index 71b32ebd3..070a318de 100644 --- a/SRC/dstemr.f +++ b/SRC/dstemr.f @@ -287,10 +287,11 @@ *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = 1X, internal error in DLARRE, -*> if INFO = 2X, internal error in DLARRV. +*> if INFO = 2X, internal error in DLARRV, +*> if INFO = 3X, internal error in DLARRJ. *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is -*> the nonzero error code returned by DLARRE or -*> DLARRV, respectively. +*> the nonzero error code returned by DLARRE, DLARRV +*> or DLARRJ, respectively. *> \endverbatim * * Authors: @@ -735,6 +736,10 @@ SUBROUTINE DSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, $ WORK( INDERR+WBEGIN-1 ), $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, $ TNRM, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = 30 + ABS( IINFO ) + RETURN + END IF IBEGIN = IEND + 1 WBEGIN = WEND + 1 39 CONTINUE diff --git a/SRC/slarrb.f b/SRC/slarrb.f index 6e43f126f..efa9e0410 100644 --- a/SRC/slarrb.f +++ b/SRC/slarrb.f @@ -165,7 +165,11 @@ *> \param[out] INFO *> \verbatim *> INFO is INTEGER -*> Error flag. +*> = 0: successful exit +*> = 1: the interval around one of the eigenvalues could not +*> be widened to contain it: the Sturm counts of the +*> representation are inconsistent, as they are when it +*> contains a NaN or an Inf. *> \endverbatim * * Authors: @@ -274,10 +278,23 @@ SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1, * * Do while( NEGCNT(LEFT).GT.I-1 ) * - BACK = WERR( II ) +* The interval is widened by BACK, doubled at every step, until +* the Sturm count agrees with the index. In exact arithmetic +* this terminates because the count is 0 below the spectrum and +* N above it. A NaN or an Inf in the representation is never +* counted, so the count can stay short of the index; stop once +* BACK has overflowed, when no further step can move the +* endpoint. BACK starts at no less than the minimum interval +* width so that a zero WERR still makes progress. +* + BACK = MAX( WERR( II ), MNWDTH ) 20 CONTINUE NEGCNT = SLANEG( N, D, LLD, LEFT, PIVMIN, R ) IF( NEGCNT.GT.I-1 ) THEN + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF LEFT = LEFT - BACK BACK = TWO*BACK GO TO 20 @@ -286,15 +303,18 @@ SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1, * Do while( NEGCNT(RIGHT).LT.I ) * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT * - BACK = WERR( II ) + BACK = MAX( WERR( II ), MNWDTH ) 50 CONTINUE - NEGCNT = SLANEG( N, D, LLD, RIGHT, PIVMIN, R ) - IF( NEGCNT.LT.I ) THEN - RIGHT = RIGHT + BACK - BACK = TWO*BACK - GO TO 50 - END IF + IF( NEGCNT.LT.I ) THEN + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF + RIGHT = RIGHT + BACK + BACK = TWO*BACK + GO TO 50 + END IF WIDTH = HALF*ABS( LEFT - RIGHT ) TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) CVRGD = MAX(RTOL1*GAP,RTOL2*TMP) diff --git a/SRC/slarrj.f b/SRC/slarrj.f index 6cfb2eec8..5fc94ca5f 100644 --- a/SRC/slarrj.f +++ b/SRC/slarrj.f @@ -137,7 +137,11 @@ *> \param[out] INFO *> \verbatim *> INFO is INTEGER -*> Error flag. +*> = 0: successful exit +*> = 1: the interval around one of the eigenvalues could not +*> be widened to contain it: the Sturm counts of the +*> matrix are inconsistent, as they are when it contains +*> a NaN or an Inf. *> \endverbatim * * Authors: @@ -190,7 +194,7 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST, * .. Local Scalars .. INTEGER CNT, I, I1, I2, II, ITER, J, K, NEXT, NINT, $ OLNINT, P, PREV, SAVI1 - REAL DPLUS, FAC, LEFT, MID, RIGHT, S, TMP, WIDTH + REAL BACK, DPLUS, LEFT, MID, RIGHT, S, TMP, WIDTH * * .. * .. Intrinsic Functions .. @@ -249,7 +253,16 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST, * * Do while( CNT(LEFT).GT.I-1 ) * - FAC = ONE +* The interval is widened by BACK, doubled at every step, +* until the Sturm count agrees with the index. In exact +* arithmetic this terminates because the count is 0 below the +* spectrum and N above it. A NaN or an Inf in the matrix is +* never counted, so the count can stay short of the index; +* stop once BACK has overflowed, when no further step can move +* the endpoint. BACK starts at no less than the minimum +* interval width so that a zero WERR still makes progress. +* + BACK = MAX( WERR( II ), TWO*PIVMIN ) 20 CONTINUE CNT = 0 S = LEFT @@ -260,14 +273,18 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST, IF( DPLUS.LT.ZERO ) CNT = CNT + 1 30 CONTINUE IF( CNT.GT.I-1 ) THEN - LEFT = LEFT - WERR( II )*FAC - FAC = TWO*FAC + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF + LEFT = LEFT - BACK + BACK = TWO*BACK GO TO 20 END IF * * Do while( CNT(RIGHT).LT.I ) * - FAC = ONE + BACK = MAX( WERR( II ), TWO*PIVMIN ) 50 CONTINUE CNT = 0 S = RIGHT @@ -278,8 +295,12 @@ SUBROUTINE SLARRJ( N, D, E2, IFIRST, ILAST, IF( DPLUS.LT.ZERO ) CNT = CNT + 1 60 CONTINUE IF( CNT.LT.I ) THEN - RIGHT = RIGHT + WERR( II )*FAC - FAC = TWO*FAC + IF( .NOT.( BACK.LT.TWO*BACK ) ) THEN + INFO = 1 + RETURN + END IF + RIGHT = RIGHT + BACK + BACK = TWO*BACK GO TO 50 END IF NINT = NINT + 1 diff --git a/SRC/slarrv.f b/SRC/slarrv.f index 59ad49770..2734fb0d9 100644 --- a/SRC/slarrv.f +++ b/SRC/slarrv.f @@ -698,6 +698,10 @@ SUBROUTINE SLARRV( N, VL, VU, D, L, PIVMIN, $ WERR(WBEGIN),WORK( INDWRK ), $ IWORK( IINDWK ), PIVMIN, SPDIAM, $ IN, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = -1 + RETURN + END IF 55 CONTINUE * IF((WBEGIN+NEWLST-1.LT.DOL).OR. diff --git a/SRC/sstemr.f b/SRC/sstemr.f index 8fe8292a5..df8f6e59c 100644 --- a/SRC/sstemr.f +++ b/SRC/sstemr.f @@ -287,10 +287,11 @@ *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = 1X, internal error in SLARRE, -*> if INFO = 2X, internal error in SLARRV. +*> if INFO = 2X, internal error in SLARRV, +*> if INFO = 3X, internal error in SLARRJ. *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is -*> the nonzero error code returned by SLARRE or -*> SLARRV, respectively. +*> the nonzero error code returned by SLARRE, SLARRV +*> or SLARRJ, respectively. *> \endverbatim * * Authors: @@ -732,6 +733,10 @@ SUBROUTINE SSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, $ WORK( INDERR+WBEGIN-1 ), $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, $ TNRM, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = 30 + ABS( IINFO ) + RETURN + END IF IBEGIN = IEND + 1 WBEGIN = WEND + 1 39 CONTINUE diff --git a/SRC/zlarrv.f b/SRC/zlarrv.f index 67c1e186a..7bcb99a35 100644 --- a/SRC/zlarrv.f +++ b/SRC/zlarrv.f @@ -702,6 +702,10 @@ SUBROUTINE ZLARRV( N, VL, VU, D, L, PIVMIN, $ WERR(WBEGIN),WORK( INDWRK ), $ IWORK( IINDWK ), PIVMIN, SPDIAM, $ IN, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = -1 + RETURN + END IF 55 CONTINUE * IF((WBEGIN+NEWLST-1.LT.DOL).OR. diff --git a/SRC/zstemr.f b/SRC/zstemr.f index f8fc18ca1..0ceae2308 100644 --- a/SRC/zstemr.f +++ b/SRC/zstemr.f @@ -304,10 +304,11 @@ *> = 0: successful exit *> < 0: if INFO = -i, the i-th argument had an illegal value *> > 0: if INFO = 1X, internal error in DLARRE, -*> if INFO = 2X, internal error in ZLARRV. +*> if INFO = 2X, internal error in ZLARRV, +*> if INFO = 3X, internal error in DLARRJ. *> Here, the digit X = ABS( IINFO ) < 10, where IINFO is -*> the nonzero error code returned by DLARRE or -*> ZLARRV, respectively. +*> the nonzero error code returned by DLARRE, ZLARRV +*> or DLARRJ, respectively. *> \endverbatim * * Authors: @@ -751,6 +752,10 @@ SUBROUTINE ZSTEMR( JOBZ, RANGE, N, D, E, VL, VU, IL, IU, $ WORK( INDERR+WBEGIN-1 ), $ WORK( INDWRK ), IWORK( IINDWK ), PIVMIN, $ TNRM, IINFO ) + IF( IINFO.NE.0 ) THEN + INFO = 30 + ABS( IINFO ) + RETURN + END IF IBEGIN = IEND + 1 WBEGIN = WEND + 1 39 CONTINUE diff --git a/TESTING/EIG/cchkst.f b/TESTING/EIG/cchkst.f index 790c7f5f1..f150cebf8 100644 --- a/TESTING/EIG/cchkst.f +++ b/TESTING/EIG/cchkst.f @@ -641,12 +641,15 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, PARAMETER ( CREL = .FALSE. ) * .. * .. Local Scalars .. + CHARACTER JOBZ + INTEGER JOBNUM LOGICAL BADNN, TRYRAC INTEGER I, IINFO, IL, IMODE, INDE, INDRWK, ITEMP, $ ITYPE, IU, J, JC, JR, JSIZE, JTYPE, LGN, $ LIWEDC, LOG2UI, LRWEDC, LWEDC, M, M2, M3, $ MTYPES, N, NAP, NBLOCK, NERRS, NMATS, NMAX, $ NSPLIT, NTEST, NTESTT + REAL RNAN, RONE REAL ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU @@ -1954,10 +1957,58 @@ SUBROUTINE CCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 310 CONTINUE * * Summary +* +* +* CSTEMR must return for a tridiagonal matrix that contains a NaN +* instead of widening a bracket forever. A NaN pivot is never +* counted, so the Sturm count of such a matrix never reaches the +* index of the wanted eigenvalue, and the widening loops of +* SLARRB, which refines the eigenvalues of a representation, and +* of SLARRJ, which refines them on the matrix itself when no +* eigenvectors are wanted, can only be ended by a bound on the +* step. Any nonnegative INFO is accepted: the test is that the +* call returns at all. +* + IF( NMAX.GE.3 .AND. LRWORK.GE.18*NMAX .AND. + $ LIWORK.GE.12*NMAX .AND. + $ ILAENV( 10, 'CSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'CSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN + N = NMAX + RONE = ONE + RNAN = SQRT( -RONE ) + DO 330 JOBNUM = 1, 2 + IF( JOBNUM.EQ.1 ) THEN + JOBZ = 'V' + ELSE + JOBZ = 'N' + END IF + DO 320 J = 1, N + SD( J ) = REAL( J ) + SE( J ) = ONE / REAL( J+1 ) + 320 CONTINUE + SD( 1 ) = RNAN + SE( N ) = ZERO + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL CSTEMR( JOBZ, 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, RWORK, + $ LRWORK, IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.LT.0 ) THEN + WRITE( NOUNIT, FMT = 9986 )JOBZ, IINFO + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + END IF * CALL SLASUM( 'CST', NOUNIT, NERRS, NTESTT ) RETURN * + 9986 FORMAT( ' CCHKST: CSTEMR( ', A1, ', A ) on a matrix with a', + $ ' NaN returned INFO=', I6 ) 9999 FORMAT( ' CCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * diff --git a/TESTING/EIG/dchkst.f b/TESTING/EIG/dchkst.f index 012aa95d4..4b4599b3d 100644 --- a/TESTING/EIG/dchkst.f +++ b/TESTING/EIG/dchkst.f @@ -625,11 +625,14 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, PARAMETER ( SREL = .FALSE. ) * .. * .. Local Scalars .. + CHARACTER JOBZ + INTEGER JOBNUM LOGICAL BADNN, TRYRAC INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC, $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, $ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS, $ NMATS, NMAX, NSPLIT, NTEST, NTESTT + DOUBLE PRECISION RNAN, RONE DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU @@ -1932,10 +1935,57 @@ SUBROUTINE DCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 310 CONTINUE * * Summary +* +* +* DSTEMR must return for a tridiagonal matrix that contains a NaN +* instead of widening a bracket forever. A NaN pivot is never +* counted, so the Sturm count of such a matrix never reaches the +* index of the wanted eigenvalue, and the widening loops of +* DLARRB, which refines the eigenvalues of a representation, and +* of DLARRJ, which refines them on the matrix itself when no +* eigenvectors are wanted, can only be ended by a bound on the +* step. Any nonnegative INFO is accepted: the test is that the +* call returns at all. +* + IF( NMAX.GE.3 .AND. LWORK.GE.18*NMAX .AND. LIWORK.GE.12*NMAX .AND. + $ ILAENV( 10, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'DSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN + N = NMAX + RONE = ONE + RNAN = SQRT( -RONE ) + DO 330 JOBNUM = 1, 2 + IF( JOBNUM.EQ.1 ) THEN + JOBZ = 'V' + ELSE + JOBZ = 'N' + END IF + DO 320 J = 1, N + SD( J ) = DBLE( J ) + SE( J ) = ONE / DBLE( J+1 ) + 320 CONTINUE + SD( 1 ) = RNAN + SE( N ) = ZERO + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL DSTEMR( JOBZ, 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, WORK, + $ LWORK, IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.LT.0 ) THEN + WRITE( NOUNIT, FMT = 9986 )JOBZ, IINFO + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + END IF * CALL DLASUM( 'DST', NOUNIT, NERRS, NTESTT ) RETURN * + 9986 FORMAT( ' DCHKST: DSTEMR( ', A1, ', A ) on a matrix with a', + $ ' NaN returned INFO=', I6 ) 9999 FORMAT( ' DCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * diff --git a/TESTING/EIG/schkst.f b/TESTING/EIG/schkst.f index c5c7f57c3..d0356f2c2 100644 --- a/TESTING/EIG/schkst.f +++ b/TESTING/EIG/schkst.f @@ -625,11 +625,14 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, PARAMETER ( SREL = .FALSE. ) * .. * .. Local Scalars .. + CHARACTER JOBZ + INTEGER JOBNUM LOGICAL BADNN, TRYRAC INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, J, JC, $ JR, JSIZE, JTYPE, LGN, LIWEDC, LOG2UI, LWEDC, $ M, M2, M3, MTYPES, N, NAP, NBLOCK, NERRS, $ NMATS, NMAX, NSPLIT, NTEST, NTESTT + REAL RNAN, RONE REAL ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU @@ -1932,10 +1935,57 @@ SUBROUTINE SCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 310 CONTINUE * * Summary +* +* +* SSTEMR must return for a tridiagonal matrix that contains a NaN +* instead of widening a bracket forever. A NaN pivot is never +* counted, so the Sturm count of such a matrix never reaches the +* index of the wanted eigenvalue, and the widening loops of +* SLARRB, which refines the eigenvalues of a representation, and +* of SLARRJ, which refines them on the matrix itself when no +* eigenvectors are wanted, can only be ended by a bound on the +* step. Any nonnegative INFO is accepted: the test is that the +* call returns at all. +* + IF( NMAX.GE.3 .AND. LWORK.GE.18*NMAX .AND. LIWORK.GE.12*NMAX .AND. + $ ILAENV( 10, 'SSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'SSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN + N = NMAX + RONE = ONE + RNAN = SQRT( -RONE ) + DO 330 JOBNUM = 1, 2 + IF( JOBNUM.EQ.1 ) THEN + JOBZ = 'V' + ELSE + JOBZ = 'N' + END IF + DO 320 J = 1, N + SD( J ) = REAL( J ) + SE( J ) = ONE / REAL( J+1 ) + 320 CONTINUE + SD( 1 ) = RNAN + SE( N ) = ZERO + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL SSTEMR( JOBZ, 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, WORK, + $ LWORK, IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.LT.0 ) THEN + WRITE( NOUNIT, FMT = 9986 )JOBZ, IINFO + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + END IF * CALL SLASUM( 'SST', NOUNIT, NERRS, NTESTT ) RETURN * + 9986 FORMAT( ' SCHKST: SSTEMR( ', A1, ', A ) on a matrix with a', + $ ' NaN returned INFO=', I6 ) 9999 FORMAT( ' SCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) * diff --git a/TESTING/EIG/zchkst.f b/TESTING/EIG/zchkst.f index 4335e15f0..c43a6496c 100644 --- a/TESTING/EIG/zchkst.f +++ b/TESTING/EIG/zchkst.f @@ -641,12 +641,15 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, PARAMETER ( CREL = .FALSE. ) * .. * .. Local Scalars .. + CHARACTER JOBZ + INTEGER JOBNUM LOGICAL BADNN, TRYRAC INTEGER I, IINFO, IL, IMODE, INDE, INDRWK, ITEMP, $ ITYPE, IU, J, JC, JR, JSIZE, JTYPE, LGN, $ LIWEDC, LOG2UI, LRWEDC, LWEDC, M, M2, M3, $ MTYPES, N, NAP, NBLOCK, NERRS, NMATS, NMAX, $ NSPLIT, NTEST, NTESTT + DOUBLE PRECISION RNAN, RONE DOUBLE PRECISION ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL, $ RTUNFL, TEMP1, TEMP2, TEMP3, TEMP4, ULP, $ ULPINV, UNFL, VL, VU @@ -1953,10 +1956,58 @@ SUBROUTINE ZCHKST( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, 310 CONTINUE * * Summary +* +* +* ZSTEMR must return for a tridiagonal matrix that contains a NaN +* instead of widening a bracket forever. A NaN pivot is never +* counted, so the Sturm count of such a matrix never reaches the +* index of the wanted eigenvalue, and the widening loops of +* DLARRB, which refines the eigenvalues of a representation, and +* of DLARRJ, which refines them on the matrix itself when no +* eigenvectors are wanted, can only be ended by a bound on the +* step. Any nonnegative INFO is accepted: the test is that the +* call returns at all. +* + IF( NMAX.GE.3 .AND. LRWORK.GE.18*NMAX .AND. + $ LIWORK.GE.12*NMAX .AND. + $ ILAENV( 10, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 .AND. + $ ILAENV( 11, 'ZSTEMR', 'VA', 1, 0, 0, 0 ).EQ.1 ) THEN + N = NMAX + RONE = ONE + RNAN = SQRT( -RONE ) + DO 330 JOBNUM = 1, 2 + IF( JOBNUM.EQ.1 ) THEN + JOBZ = 'V' + ELSE + JOBZ = 'N' + END IF + DO 320 J = 1, N + SD( J ) = DBLE( J ) + SE( J ) = ONE / DBLE( J+1 ) + 320 CONTINUE + SD( 1 ) = RNAN + SE( N ) = ZERO + VL = ZERO + VU = ZERO + IL = 0 + IU = 0 + TRYRAC = .TRUE. + CALL ZSTEMR( JOBZ, 'A', N, SD, SE, VL, VU, IL, IU, M, + $ WR, Z, LDU, N, IWORK( 1 ), TRYRAC, RWORK, + $ LRWORK, IWORK( 2*N+1 ), LIWORK-2*N, IINFO ) + IF( IINFO.LT.0 ) THEN + WRITE( NOUNIT, FMT = 9986 )JOBZ, IINFO + NERRS = NERRS + 1 + END IF + NTESTT = NTESTT + 1 + 330 CONTINUE + END IF * CALL DLASUM( 'ZST', NOUNIT, NERRS, NTESTT ) RETURN * + 9986 FORMAT( ' ZCHKST: ZSTEMR( ', A1, ', A ) on a matrix with a', + $ ' NaN returned INFO=', I6 ) 9999 FORMAT( ' ZCHKST: ', A, ' returned INFO=', I6, '.', / 9X, 'N=', $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' ) *