SUBROUTINE DLARRD( RANGE, ORDER, N, VL, VU, IL, IU, GERS, $ RELTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, $ M, W, WERR, WL, WU, IBLOCK, INDEXW, $ WORK, IWORK, INFO ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER ORDER, RANGE INTEGER IL, INFO, IU, M, N, NSPLIT DOUBLE PRECISION PIVMIN, RELTOL, VL, VU, WL, WU * .. * .. Array Arguments .. INTEGER IBLOCK( * ), INDEXW( * ), $ ISPLIT( * ), IWORK( * ) DOUBLE PRECISION D( * ), E( * ), E2( * ), $ GERS( * ), W( * ), WERR( * ), WORK( * ) * .. * * Purpose * ======= * * DLARRD computes the eigenvalues of a symmetric tridiagonal * matrix T to suitable accuracy. This is an auxiliary code to be * called from DSTEMR. * The user may ask for all eigenvalues, all eigenvalues * in the half-open interval (VL, VU], or the IL-th through IU-th * eigenvalues. * * To avoid overflow, the matrix must be scaled so that its * largest element is no greater than overflow**(1/2) * * underflow**(1/4) in absolute value, and for greatest * accuracy, it should not be much smaller than that. * * See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal * Matrix", Report CS41, Computer Science Dept., Stanford * University, July 21, 1966. * * Arguments * ========= * * RANGE (input) CHARACTER * = 'A': ("All") all eigenvalues will be found. * = 'V': ("Value") all eigenvalues in the half-open interval * (VL, VU] will be found. * = 'I': ("Index") the IL-th through IU-th eigenvalues (of the * entire matrix) will be found. * * ORDER (input) CHARACTER * = 'B': ("By Block") the eigenvalues will be grouped by * split-off block (see IBLOCK, ISPLIT) and * ordered from smallest to largest within * the block. * = 'E': ("Entire matrix") * the eigenvalues for the entire matrix * will be ordered from smallest to * largest. * * N (input) INTEGER * The order of the tridiagonal matrix T. N >= 0. * * VL (input) DOUBLE PRECISION * VU (input) DOUBLE PRECISION * If RANGE='V', the lower and upper bounds of the interval to * be searched for eigenvalues. Eigenvalues less than or equal * to VL, or greater than VU, will not be returned. VL < VU. * Not referenced if RANGE = 'A' or 'I'. * * IL (input) INTEGER * IU (input) INTEGER * If RANGE='I', the indices (in ascending order) of the * smallest and largest eigenvalues to be returned. * 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. * Not referenced if RANGE = 'A' or 'V'. * * GERS (input) DOUBLE PRECISION array, dimension (2*N) * The N Gerschgorin intervals (the i-th Gerschgorin interval * is (GERS(2*i-1), GERS(2*i)). * * RELTOL (input) DOUBLE PRECISION * The minimum relative width of an interval. When an interval * is narrower than RELTOL times the larger (in * magnitude) endpoint, then it is considered to be * sufficiently small, i.e., converged. Note: this should * always be at least radix*machine epsilon. * * D (input) DOUBLE PRECISION array, dimension (N) * The n diagonal elements of the tridiagonal matrix T. * * E (input) DOUBLE PRECISION array, dimension (N-1) * The (n-1) off-diagonal elements of the tridiagonal matrix T. * * E2 (input) DOUBLE PRECISION array, dimension (N-1) * The (n-1) squared off-diagonal elements of the tridiagonal matrix T. * * PIVMIN (input) DOUBLE PRECISION * The minimum pivot allowed in the Sturm sequence for T. * * NSPLIT (input) INTEGER * The number of diagonal blocks in the matrix T. * 1 <= NSPLIT <= N. * * ISPLIT (input) INTEGER array, dimension (N) * The splitting points, at which T breaks up into submatrices. * The first submatrix consists of rows/columns 1 to ISPLIT(1), * the second of rows/columns ISPLIT(1)+1 through ISPLIT(2), * etc., and the NSPLIT-th consists of rows/columns * ISPLIT(NSPLIT-1)+1 through ISPLIT(NSPLIT)=N. * (Only the first NSPLIT elements will actually be used, but * since the user cannot know a priori what value NSPLIT will * have, N words must be reserved for ISPLIT.) * * M (output) INTEGER * The actual number of eigenvalues found. 0 <= M <= N. * (See also the description of INFO=2,3.) * * W (output) DOUBLE PRECISION array, dimension (N) * On exit, the first M elements of W will contain the * eigenvalue approximations. DLARRD computes an interval * I_j = (a_j, b_j] that includes eigenvalue j. The eigenvalue * approximation is given as the interval midpoint * W(j)= ( a_j + b_j)/2. The corresponding error is bounded by * WERR(j) = abs( a_j - b_j)/2 * * WERR (output) DOUBLE PRECISION array, dimension (N) * The error bound on the corresponding eigenvalue approximation * in W. * * WL (output) DOUBLE PRECISION * WU (output) DOUBLE PRECISION * The interval (WL, WU] contains all the wanted eigenvalues. * If RANGE='V', then WL=VL and WU=VU. * If RANGE='A', then WL and WU are the global Gerschgorin bounds * on the spectrum. * If RANGE='I', then WL and WU are computed by DLAEBZ from the * index range specified. * * IBLOCK (output) INTEGER array, dimension (N) * At each row/column j where E(j) is zero or small, the * matrix T is considered to split into a block diagonal * matrix. On exit, if INFO = 0, IBLOCK(i) specifies to which * block (from 1 to the number of blocks) the eigenvalue W(i) * belongs. (DLARRD may use the remaining N-M elements as * workspace.) * * INDEXW (output) INTEGER array, dimension (N) * The indices of the eigenvalues within each block (submatrix); * for example, INDEXW(i)= j and IBLOCK(i)=k imply that the * i-th eigenvalue W(i) is the j-th eigenvalue in block k. * * WORK (workspace) DOUBLE PRECISION array, dimension (4*N) * * IWORK (workspace) INTEGER array, dimension (3*N) * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * > 0: some or all of the eigenvalues failed to converge or * were not computed: * =1 or 3: Bisection failed to converge for some * eigenvalues; these eigenvalues are flagged by a * negative block number. The effect is that the * eigenvalues may not be as accurate as the * absolute and relative tolerances. This is * generally caused by unexpectedly inaccurate * arithmetic. * =2 or 3: RANGE='I' only: Not all of the eigenvalues * IL:IU were found. * Effect: M < IU+1-IL * Cause: non-monotonic arithmetic, causing the * Sturm sequence to be non-monotonic. * Cure: recalculate, using RANGE='A', and pick * out eigenvalues IL:IU. In some cases, * increasing the PARAMETER "FUDGE" may * make things work. * = 4: RANGE='I', and the Gershgorin interval * initially used was too small. No eigenvalues * were computed. * Probable cause: your machine has sloppy * floating-point arithmetic. * Cure: Increase the PARAMETER "FUDGE", * recompile, and try again. * * Internal Parameters * =================== * * FUDGE DOUBLE PRECISION, default = 2 * A "fudge factor" to widen the Gershgorin intervals. Ideally, * a value of 1 should work, but on machines with sloppy * arithmetic, this needs to be larger. The default for * publicly released versions should be large enough to handle * the worst machine around. Note that this has no effect * on accuracy of the solution. * * Based on contributions by * W. Kahan, University of California, Berkeley, USA * Beresford Parlett, University of California, Berkeley, USA * Jim Demmel, University of California, Berkeley, USA * Inderjit Dhillon, University of Texas, Austin, USA * Osni Marques, LBNL/NERSC, USA * Christof Voemel, University of California, Berkeley, USA * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE, TWO, HALF, FUDGE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, $ TWO = 2.0D0, HALF = ONE/TWO, $ FUDGE = TWO ) INTEGER ALLRNG, VALRNG, INDRNG PARAMETER ( ALLRNG = 1, VALRNG = 2, INDRNG = 3 ) * .. * .. Local Scalars .. LOGICAL NCNVRG, TOOFEW INTEGER I, IB, IBEGIN, IDISCL, IDISCU, IE, IEND, IINFO, $ IM, IN, IOFF, IOUT, IRANGE, ITMAX, ITMP1, $ ITMP2, IW, IWOFF, J, JBLK, JDISC, JE, JEE, NB, $ NWL, NWU DOUBLE PRECISION ATOLI, EPS, GL, GU, RTOLI, SPDIAM, TMP1, TMP2, $ TNORM, UFLOW, WKILL, WLU, WUL * .. * .. Local Arrays .. INTEGER IDUMMA( 1 ) * .. * .. External Functions .. LOGICAL LSAME INTEGER ILAENV DOUBLE PRECISION DLAMCH EXTERNAL LSAME, ILAENV, DLAMCH * .. * .. External Subroutines .. EXTERNAL DLAEBZ * .. * .. Intrinsic Functions .. INTRINSIC ABS, INT, LOG, MAX, MIN * .. * .. Executable Statements .. * INFO = 0 * * Decode RANGE * IF( LSAME( RANGE, 'A' ) ) THEN IRANGE = ALLRNG ELSE IF( LSAME( RANGE, 'V' ) ) THEN IRANGE = VALRNG ELSE IF( LSAME( RANGE, 'I' ) ) THEN IRANGE = INDRNG ELSE IRANGE = 0 END IF * * Check for Errors * IF( IRANGE.LE.0 ) THEN INFO = -1 ELSE IF( .NOT.(LSAME(ORDER,'B').OR.LSAME(ORDER,'E')) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( IRANGE.EQ.VALRNG ) THEN IF( VL.GE.VU ) $ INFO = -5 ELSE IF( IRANGE.EQ.INDRNG .AND. $ ( IL.LT.1 .OR. IL.GT.MAX( 1, N ) ) ) THEN INFO = -6 ELSE IF( IRANGE.EQ.INDRNG .AND. $ ( IU.LT.MIN( N, IL ) .OR. IU.GT.N ) ) THEN INFO = -7 END IF * IF( INFO.NE.0 ) THEN RETURN END IF * Initialize error flags INFO = 0 NCNVRG = .FALSE. TOOFEW = .FALSE. * Quick return if possible M = 0 IF( N.EQ.0 ) RETURN * Simplification: IF( IRANGE.EQ.INDRNG .AND. IL.EQ.1 .AND. IU.EQ.N ) IRANGE = 1 * Get machine constants EPS = DLAMCH( 'P' ) UFLOW = DLAMCH( 'U' ) * Special Case when N=1 * Treat case of 1x1 matrix for quick return IF( N.EQ.1 ) THEN IF( (IRANGE.EQ.ALLRNG).OR. $ ((IRANGE.EQ.VALRNG).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. $ ((IRANGE.EQ.INDRNG).AND.(IL.EQ.1).AND.(IU.EQ.1)) ) THEN M = 1 W(1) = D(1) * The computation error of the eigenvalue is zero WERR(1) = ZERO IBLOCK( 1 ) = 1 INDEXW( 1 ) = 1 ENDIF RETURN END IF * NB is the minimum vector length for vector bisection, or 0 * if only scalar is to be done. NB = ILAENV( 1, 'DSTEBZ', ' ', N, -1, -1, -1 ) IF( NB.LE.1 ) NB = 0 * Find global spectral radius GL = D(1) GU = D(1) DO 5 I = 1,N GL = MIN( GL, GERS( 2*I - 1)) GU = MAX( GU, GERS(2*I) ) 5 CONTINUE * Compute global Gerschgorin bounds and spectral diameter TNORM = MAX( ABS( GL ), ABS( GU ) ) GL = GL - FUDGE*TNORM*EPS*N - FUDGE*TWO*PIVMIN GU = GU + FUDGE*TNORM*EPS*N + FUDGE*TWO*PIVMIN SPDIAM = GU - GL * Input arguments for DLAEBZ: * The relative tolerance. An interval (a,b] lies within * "relative tolerance" if b-a < RELTOL*max(|a|,|b|), RTOLI = RELTOL * Set the absolute tolerance for interval convergence to zero to force * interval convergence based on relative size of the interval. * This is dangerous because intervals might not converge when RELTOL is * small. But at least a very small number should be selected so that for * strongly graded matrices, the code can get relatively accurate * eigenvalues. ATOLI = FUDGE*TWO*UFLOW + FUDGE*TWO*PIVMIN IF( IRANGE.EQ.INDRNG ) THEN * RANGE='I': Compute an interval containing eigenvalues * IL through IU. The initial interval [GL,GU] from the global * Gerschgorin bounds GL and GU is refined by DLAEBZ. ITMAX = INT( ( LOG( TNORM+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 WORK( N+1 ) = GL WORK( N+2 ) = GL WORK( N+3 ) = GU WORK( N+4 ) = GU WORK( N+5 ) = GL WORK( N+6 ) = GU IWORK( 1 ) = -1 IWORK( 2 ) = -1 IWORK( 3 ) = N + 1 IWORK( 4 ) = N + 1 IWORK( 5 ) = IL - 1 IWORK( 6 ) = IU * CALL DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, $ D, E, E2, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT, $ IWORK, W, IBLOCK, IINFO ) IF( IINFO .NE. 0 ) THEN INFO = IINFO RETURN END IF * On exit, output intervals may not be ordered by ascending negcount IF( IWORK( 6 ).EQ.IU ) THEN WL = WORK( N+1 ) WLU = WORK( N+3 ) NWL = IWORK( 1 ) WU = WORK( N+4 ) WUL = WORK( N+2 ) NWU = IWORK( 4 ) ELSE WL = WORK( N+2 ) WLU = WORK( N+4 ) NWL = IWORK( 2 ) WU = WORK( N+3 ) WUL = WORK( N+1 ) NWU = IWORK( 3 ) END IF * On exit, the interval [WL, WLU] contains a value with negcount NWL, * and [WUL, WU] contains a value with negcount NWU. IF( NWL.LT.0 .OR. NWL.GE.N .OR. NWU.LT.1 .OR. NWU.GT.N ) THEN INFO = 4 RETURN END IF ELSEIF( IRANGE.EQ.VALRNG ) THEN WL = VL WU = VU ELSEIF( IRANGE.EQ.ALLRNG ) THEN WL = GL WU = GU ENDIF * Find Eigenvalues -- Loop Over blocks and recompute NWL and NWU. * NWL accumulates the number of eigenvalues .le. WL, * NWU accumulates the number of eigenvalues .le. WU M = 0 IEND = 0 INFO = 0 NWL = 0 NWU = 0 * DO 70 JBLK = 1, NSPLIT IOFF = IEND IBEGIN = IOFF + 1 IEND = ISPLIT( JBLK ) IN = IEND - IOFF * IF( IN.EQ.1 ) THEN * 1x1 block IF( WL.GE.D( IBEGIN )-PIVMIN ) $ NWL = NWL + 1 IF( WU.GE.D( IBEGIN )-PIVMIN ) $ NWU = NWU + 1 IF( IRANGE.EQ.ALLRNG .OR. $ ( WL.LT.D( IBEGIN )-PIVMIN $ .AND. WU.GE. D( IBEGIN )-PIVMIN ) ) THEN M = M + 1 W( M ) = D( IBEGIN ) WERR(M) = ZERO * The gap for a single block doesn't matter for the later * algorithm and is assigned an arbitrary large value IBLOCK( M ) = JBLK INDEXW( M ) = 1 END IF * Disabled 2x2 case because of a failure on the following matrix * RANGE = 'I', IL = IU = 4 * Original Tridiagonal, d = [ * -0.150102010615740E+00 * -0.849897989384260E+00 * -0.128208148052635E-15 * 0.128257718286320E-15 * ]; * e = [ * -0.357171383266986E+00 * -0.180411241501588E-15 * -0.175152352710251E-15 * ]; * * ELSE IF( IN.EQ.2 ) THEN ** 2x2 block * DISC = SQRT( (HALF*(D(IBEGIN)-D(IEND)))**2 + E(IBEGIN)**2 ) * TMP1 = HALF*(D(IBEGIN)+D(IEND)) * L1 = TMP1 - DISC * IF( WL.GE. L1-PIVMIN ) * $ NWL = NWL + 1 * IF( WU.GE. L1-PIVMIN ) * $ NWU = NWU + 1 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L1-PIVMIN .AND. WU.GE. * $ L1-PIVMIN ) ) THEN * M = M + 1 * W( M ) = L1 ** The uncertainty of eigenvalues of a 2x2 matrix is very small * WERR( M ) = EPS * ABS( W( M ) ) * TWO * IBLOCK( M ) = JBLK * INDEXW( M ) = 1 * ENDIF * L2 = TMP1 + DISC * IF( WL.GE. L2-PIVMIN ) * $ NWL = NWL + 1 * IF( WU.GE. L2-PIVMIN ) * $ NWU = NWU + 1 * IF( IRANGE.EQ.ALLRNG .OR. ( WL.LT.L2-PIVMIN .AND. WU.GE. * $ L2-PIVMIN ) ) THEN * M = M + 1 * W( M ) = L2 ** The uncertainty of eigenvalues of a 2x2 matrix is very small * WERR( M ) = EPS * ABS( W( M ) ) * TWO * IBLOCK( M ) = JBLK * INDEXW( M ) = 2 * ENDIF ELSE * General Case - block of size IN >= 2 * Compute local Gerschgorin interval and use it as the initial * interval for DLAEBZ GU = D( IBEGIN ) GL = D( IBEGIN ) TMP1 = ZERO DO 40 J = IBEGIN, IEND GL = MIN( GL, GERS( 2*J - 1)) GU = MAX( GU, GERS(2*J) ) 40 CONTINUE SPDIAM = GU - GL GL = GL - FUDGE*SPDIAM*EPS*IN - FUDGE*PIVMIN GU = GU + FUDGE*SPDIAM*EPS*IN + FUDGE*PIVMIN * IF( IRANGE.GT.1 ) THEN IF( GU.LT.WL ) THEN * the local block contains none of the wanted eigenvalues NWL = NWL + IN NWU = NWU + IN GO TO 70 END IF * refine search interval if possible, only range (WL,WU] matters GL = MAX( GL, WL ) GU = MIN( GU, WU ) IF( GL.GE.GU ) $ GO TO 70 END IF * Find negcount of initial interval boundaries GL and GU WORK( N+1 ) = GL WORK( N+IN+1 ) = GU CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) IF( IINFO .NE. 0 ) THEN INFO = IINFO RETURN END IF * NWL = NWL + IWORK( 1 ) NWU = NWU + IWORK( IN+1 ) IWOFF = M - IWORK( 1 ) * Compute Eigenvalues ITMAX = INT( ( LOG( GU-GL+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN, $ D( IBEGIN ), E( IBEGIN ), E2( IBEGIN ), $ IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT, $ IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO ) IF( IINFO .NE. 0 ) THEN INFO = IINFO RETURN END IF * * Copy eigenvalues into W and IBLOCK * Use -JBLK for block number for unconverged eigenvalues. * Loop over the number of output intervals from DLAEBZ DO 60 J = 1, IOUT * eigenvalue approximation is middle point of interval TMP1 = HALF*( WORK( J+N )+WORK( J+IN+N ) ) * semi length of error interval TMP2 = HALF*ABS( WORK( J+N )-WORK( J+IN+N ) ) IF( J.GT.IOUT-IINFO ) THEN * Flag non-convergence. NCNVRG = .TRUE. IB = -JBLK ELSE IB = JBLK END IF DO 50 JE = IWORK( J ) + 1 + IWOFF, $ IWORK( J+IN ) + IWOFF W( JE ) = TMP1 WERR( JE ) = TMP2 INDEXW( JE ) = JE - IWOFF IBLOCK( JE ) = IB 50 CONTINUE 60 CONTINUE * M = M + IM END IF 70 CONTINUE * If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU * If NWL+1 < IL or NWU > IU, discard extra eigenvalues. IF( IRANGE.EQ.INDRNG ) THEN IDISCL = IL - 1 - NWL IDISCU = NWU - IU * IF( IDISCL.GT.0 ) THEN IM = 0 DO 80 JE = 1, M * Remove some of the smallest eigenvalues from the left so that * at the end IDISCL =0. Move all eigenvalues up to the left. IF( W( JE ).LE.WLU .AND. IDISCL.GT.0 ) THEN IDISCL = IDISCL - 1 ELSE IM = IM + 1 W( IM ) = W( JE ) WERR( IM ) = WERR( JE ) INDEXW( IM ) = INDEXW( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 80 CONTINUE M = IM END IF IF( IDISCU.GT.0 ) THEN * Remove some of the largest eigenvalues from the right so that * at the end IDISCU =0. Move all eigenvalues up to the left. IM=M+1 DO 81 JE = M, 1, -1 IF( W( JE ).GE.WUL .AND. IDISCU.GT.0 ) THEN IDISCU = IDISCU - 1 ELSE IM = IM - 1 W( IM ) = W( JE ) WERR( IM ) = WERR( JE ) INDEXW( IM ) = INDEXW( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 81 CONTINUE JEE = 0 DO 82 JE = IM, M JEE = JEE + 1 W( JEE ) = W( JE ) WERR( JEE ) = WERR( JE ) INDEXW( JEE ) = INDEXW( JE ) IBLOCK( JEE ) = IBLOCK( JE ) 82 CONTINUE M = M-IM+1 END IF IF( IDISCL.GT.0 .OR. IDISCU.GT.0 ) THEN * Code to deal with effects of bad arithmetic. (If N(w) is * monotone non-decreasing, this should never happen.) * Some low eigenvalues to be discarded are not in (WL,WLU], * or high eigenvalues to be discarded are not in (WUL,WU] * so just kill off the smallest IDISCL/largest IDISCU * eigenvalues, by marking the corresponding IBLOCK = 0 IF( IDISCL.GT.0 ) THEN WKILL = WU DO 100 JDISC = 1, IDISCL IW = 0 DO 90 JE = 1, M IF( IBLOCK( JE ).NE.0 .AND. $ ( W( JE ).LT.WKILL .OR. IW.EQ.0 ) ) THEN IW = JE WKILL = W( JE ) END IF 90 CONTINUE IBLOCK( IW ) = 0 100 CONTINUE END IF IF( IDISCU.GT.0 ) THEN WKILL = WL DO 120 JDISC = 1, IDISCU IW = 0 DO 110 JE = 1, M IF( IBLOCK( JE ).NE.0 .AND. $ ( W( JE ).GE.WKILL .OR. IW.EQ.0 ) ) THEN IW = JE WKILL = W( JE ) END IF 110 CONTINUE IBLOCK( IW ) = 0 120 CONTINUE END IF * Now erase all eigenvalues with IBLOCK set to zero IM = 0 DO 130 JE = 1, M IF( IBLOCK( JE ).NE.0 ) THEN IM = IM + 1 W( IM ) = W( JE ) WERR( IM ) = WERR( JE ) INDEXW( IM ) = INDEXW( JE ) IBLOCK( IM ) = IBLOCK( JE ) END IF 130 CONTINUE M = IM END IF IF( IDISCL.LT.0 .OR. IDISCU.LT.0 ) THEN TOOFEW = .TRUE. END IF END IF * IF(( IRANGE.EQ.ALLRNG .AND. M.NE.N ).OR. $ ( IRANGE.EQ.INDRNG .AND. M.NE.IU-IL+1 ) ) THEN TOOFEW = .TRUE. END IF * If ORDER='B', do nothing the eigenvalues are already sorted by * block. * If ORDER='E', sort the eigenvalues from smallest to largest IF( LSAME(ORDER,'E') .AND. NSPLIT.GT.1 ) THEN DO 150 JE = 1, M - 1 IE = 0 TMP1 = W( JE ) DO 140 J = JE + 1, M IF( W( J ).LT.TMP1 ) THEN IE = J TMP1 = W( J ) END IF 140 CONTINUE IF( IE.NE.0 ) THEN TMP2 = WERR( IE ) ITMP1 = IBLOCK( IE ) ITMP2 = INDEXW( IE ) W( IE ) = W( JE ) WERR( IE ) = WERR( JE ) IBLOCK( IE ) = IBLOCK( JE ) INDEXW( IE ) = INDEXW( JE ) W( JE ) = TMP1 WERR( JE ) = TMP2 IBLOCK( JE ) = ITMP1 INDEXW( JE ) = ITMP2 END IF 150 CONTINUE END IF * INFO = 0 IF( NCNVRG ) $ INFO = INFO + 1 IF( TOOFEW ) $ INFO = INFO + 2 RETURN * * End of DLARRD * END