SUBROUTINE DLARRE2A( RANGE, N, VL, VU, IL, IU, D, E, E2, $ RTOL1, RTOL2, SPLTOL, NSPLIT, ISPLIT, $ M, DOL, DOU, NEEDIL, NEEDIU, $ W, WERR, WGAP, IBLOCK, INDEXW, GERS, $ SDIAM, PIVMIN, WORK, IWORK, MINRGP, INFO ) * * -- ScaLAPACK auxiliary routine (version 2.0) -- * Univ. of Tennessee, Univ. of California Berkeley, Univ of Colorado Denver * July 4, 2010 * IMPLICIT NONE * * .. Scalar Arguments .. CHARACTER RANGE INTEGER DOL, DOU, IL, INFO, IU, M, N, NSPLIT, $ NEEDIL, NEEDIU DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, SPLTOL, VL, VU * .. * .. Array Arguments .. INTEGER IBLOCK( * ), ISPLIT( * ), IWORK( * ), $ INDEXW( * ) DOUBLE PRECISION D( * ), E( * ), E2( * ), GERS( * ), $ SDIAM( * ), W( * ),WERR( * ), $ WGAP( * ), WORK( * ) * * Purpose * ======= * * To find the desired eigenvalues of a given real symmetric * tridiagonal matrix T, DLARRE2 sets any "small" off-diagonal * elements to zero, and for each unreduced block T_i, it finds * (a) a suitable shift at one end of the block's spectrum, * (b) the base representation, T_i - sigma_i I = L_i D_i L_i^T, and * (c) eigenvalues of each L_i D_i L_i^T. * * NOTE: * The algorithm obtains a crude picture of all the wanted eigenvalues * (as selected by RANGE). However, to reduce work and improve scalability, * only the eigenvalues DOL to DOU are refined. Furthermore, if the matrix * splits into blocks, RRRs for blocks that do not contain eigenvalues * from DOL to DOU are skipped. * The DQDS algorithm (subroutine DLASQ2) is not used, unlike in * the sequential case. Instead, eigenvalues are computed in parallel to some * figures using bisection. * * 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. * * N (input) INTEGER * The order of the matrix. N > 0. * * VL (input/output) DOUBLE PRECISION * VU (input/output) DOUBLE PRECISION * If RANGE='V', the lower and upper bounds for the eigenvalues. * Eigenvalues less than or equal to VL, or greater than VU, * will not be returned. VL < VU. * If RANGE='I' or ='A', DLARRE2A computes bounds on the desired * part of the spectrum. * * 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. * * D (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the N diagonal elements of the tridiagonal * matrix T. * On exit, the N diagonal elements of the diagonal * matrices D_i. * * E (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the first (N-1) entries contain the subdiagonal * elements of the tridiagonal matrix T; E(N) need not be set. * On exit, E contains the subdiagonal elements of the unit * bidiagonal matrices L_i. The entries E( ISPLIT( I ) ), * 1 <= I <= NSPLIT, contain the base points sigma_i on output. * * E2 (input/output) DOUBLE PRECISION array, dimension (N) * On entry, the first (N-1) entries contain the SQUARES of the * subdiagonal elements of the tridiagonal matrix T; * E2(N) need not be set. * On exit, the entries E2( ISPLIT( I ) ), * 1 <= I <= NSPLIT, have been set to zero * * RTOL1 (input) DOUBLE PRECISION * RTOL2 (input) DOUBLE PRECISION * Parameters for bisection. * An interval [LEFT,RIGHT] has converged if * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) * * SPLTOL (input) DOUBLE PRECISION * The threshold for splitting. * * NSPLIT (output) INTEGER * The number of blocks T splits into. 1 <= NSPLIT <= N. * * ISPLIT (output) INTEGER array, dimension (N) * The splitting points, at which T breaks up into blocks. * The first block 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. * * M (output) INTEGER * The total number of eigenvalues (of all L_i D_i L_i^T) * found. * * DOL (input) INTEGER * DOU (input) INTEGER * If the user wants to work on only a selected part of the * representation tree, he can specify an index range DOL:DOU. * Otherwise, the setting DOL=1, DOU=N should be applied. * Note that DOL and DOU refer to the order in which the eigenvalues * are stored in W. * * NEEDIL (output) INTEGER * NEEDIU (output) INTEGER * The indices of the leftmost and rightmost eigenvalues * of the root node RRR which are * needed to accurately compute the relevant part of the * representation tree. * * W (output) DOUBLE PRECISION array, dimension (N) * The first M elements contain the eigenvalues. The * eigenvalues of each of the blocks, L_i D_i L_i^T, are * sorted in ascending order ( DLARRE2A may use the * remaining N-M elements as workspace). * Note that immediately after exiting this routine, only * the eigenvalues from position DOL:DOU in W are * reliable on this processor * because the eigenvalue computation is done in parallel. * * WERR (output) DOUBLE PRECISION array, dimension (N) * The error bound on the corresponding eigenvalue in W. * Note that immediately after exiting this routine, only * the uncertainties from position DOL:DOU in WERR are * reliable on this processor * because the eigenvalue computation is done in parallel. * * WGAP (output) DOUBLE PRECISION array, dimension (N) * The separation from the right neighbor eigenvalue in W. * The gap is only with respect to the eigenvalues of the same block * as each block has its own representation tree. * Exception: at the right end of a block we store the left gap * Note that immediately after exiting this routine, only * the gaps from position DOL:DOU in WGAP are * reliable on this processor * because the eigenvalue computation is done in parallel. * * IBLOCK (output) INTEGER array, dimension (N) * The indices of the blocks (submatrices) associated with the * corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue * W(i) belongs to the first block from the top, =2 if W(i) * belongs to the second block, etc. * * INDEXW (output) INTEGER array, dimension (N) * The indices of the eigenvalues within each block (submatrix); * for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the * i-th eigenvalue W(i) is the 10-th eigenvalue in block 2 * * GERS (output) DOUBLE PRECISION array, dimension (2*N) * The N Gerschgorin intervals (the i-th Gerschgorin interval * is (GERS(2*i-1), GERS(2*i)). * * PIVMIN (output) DOUBLE PRECISION * The minimum pivot in the sturm sequence for T. * * WORK (workspace) DOUBLE PRECISION array, dimension (6*N) * Workspace. * * IWORK (workspace) INTEGER array, dimension (5*N) * Workspace. * * MINRGP (input) DOUBLE PRECISION * The minimum relativ gap threshold to decide whether an eigenvalue * or a cluster boundary is reached. * * INFO (output) INTEGER * = 0: successful exit * > 0: A problem occured in DLARRE2A. * < 0: One of the called subroutines signaled an internal problem. * Needs inspection of the corresponding parameter IINFO * for further information. * * =-1: Problem in DLARRD2. * = 2: No base representation could be found in MAXTRY iterations. * Increasing MAXTRY and recompilation might be a remedy. * =-3: Problem in DLARRB2 when computing the refined root * representation * =-4: Problem in DLARRB2 when preforming bisection on the * desired part of the spectrum. * = -9 Problem: M < DOU-DOL+1, that is the code found fewer * eigenvalues than it was supposed to * * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION FAC, FOUR, FOURTH, FUDGE, HALF, HNDRD, $ MAXGROWTH, ONE, PERT, TWO, ZERO PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, $ TWO = 2.0D0, FOUR=4.0D0, $ HNDRD = 100.0D0, $ PERT = 8.0D0, $ HALF = ONE/TWO, FOURTH = ONE/FOUR, FAC= HALF, $ MAXGROWTH = 64.0D0, FUDGE = 2.0D0 ) INTEGER MAXTRY PARAMETER ( MAXTRY = 6 ) * .. * .. Local Scalars .. LOGICAL NOREP, RNDPRT, USEDQD INTEGER CNT, CNT1, CNT2, I, IBEGIN, IDUM, IEND, IINFO, $ IN, INDL, INDU, IRANGE, J, JBLK, MB, MM, $ MYINDL, MYINDU, MYWBEG, MYWEND, WBEGIN, WEND DOUBLE PRECISION AVGAP, BSRTOL, CLWDTH, DMAX, DPIVOT, EABS, $ EMAX, EOLD, EPS, GL, GU, ISLEFT, ISRGHT, $ LGPVMN, LGSPDM, RTL, S1, S2, SAFMIN, SGNDEF, $ SIGMA, SPDIAM, TAU, TMP, TMP1, TMP2 * .. * .. Local Arrays .. INTEGER ISEED( 4 ) * .. * .. External Functions .. LOGICAL LSAME DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH, LSAME * .. * .. External Subroutines .. EXTERNAL DCOPY, DLARNV, DLARRA, DLARRB2, $ DLARRC, DLARRD2 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Executable Statements .. * INFO = 0 * Dis-/Enable a small random perturbation of the root representation RNDPRT = .TRUE. * * Decode RANGE * IF( LSAME( RANGE, 'A' ) ) THEN IRANGE = 1 ELSE IF( LSAME( RANGE, 'V' ) ) THEN IRANGE = 2 ELSE IF( LSAME( RANGE, 'I' ) ) THEN IRANGE = 3 END IF M = 0 * Get machine constants SAFMIN = DLAMCH( 'S' ) EPS = DLAMCH( 'P' ) * Set parameters RTL = SQRT(EPS) BSRTOL = 1.0D-1 * Treat case of 1x1 matrix for quick return IF( N.EQ.1 ) THEN IF( (IRANGE.EQ.1).OR. $ ((IRANGE.EQ.2).AND.(D(1).GT.VL).AND.(D(1).LE.VU)).OR. $ ((IRANGE.EQ.3).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 WGAP(1) = ZERO IBLOCK( 1 ) = 1 INDEXW( 1 ) = 1 GERS(1) = D( 1 ) GERS(2) = D( 1 ) ENDIF * store the shift for the initial RRR, which is zero in this case E(1) = ZERO RETURN END IF * General case: tridiagonal matrix of order > 1 * Init WERR, WGAP. DO 1 I =1,N WERR(I) = ZERO 1 CONTINUE DO 2 I =1,N WGAP(I) = ZERO 2 CONTINUE * Compute Gerschgorin intervals and spectral diameter. * Compute maximum off-diagonal entry and pivmin. GL = D(1) GU = D(1) EOLD = ZERO EMAX = ZERO E(N) = ZERO DO 5 I = 1,N EABS = ABS( E(I) ) IF( EABS .GE. EMAX ) THEN EMAX = EABS END IF TMP = EABS + EOLD EOLD = EABS TMP1 = D(I) - TMP TMP2 = D(I) + TMP GL = MIN( GL, TMP1 ) GU = MAX( GU, TMP2 ) GERS( 2*I-1) = TMP1 GERS( 2*I ) = TMP2 5 CONTINUE * The minimum pivot allowed in the sturm sequence for T PIVMIN = SAFMIN * MAX( ONE, EMAX**2 ) * Compute spectral diameter. The Gerschgorin bounds give an * estimate that is wrong by at most a factor of SQRT(2) SPDIAM = GU - GL * Compute splitting points CALL DLARRA( N, D, E, E2, SPLTOL, SPDIAM, $ NSPLIT, ISPLIT, IINFO ) IF( IRANGE.EQ.1 ) THEN * Set interval [VL,VU] that contains all eigenvalues VL = GL VU = GU ENDIF * We call DLARRD2 to find crude approximations to the eigenvalues * in the desired range. In case IRANGE = 3, we also obtain the * interval (VL,VU] that contains all the wanted eigenvalues. * An interval [LEFT,RIGHT] has converged if * RIGHT-LEFT.LT.RTOL*MAX(ABS(LEFT),ABS(RIGHT)) * DLARRD2 needs a WORK of size 4*N, IWORK of size 3*N CALL DLARRD2( RANGE, 'B', N, VL, VU, IL, IU, GERS, $ BSRTOL, D, E, E2, PIVMIN, NSPLIT, ISPLIT, $ MM, W, WERR, VL, VU, IBLOCK, INDEXW, $ WORK, IWORK, DOL, DOU, IINFO ) IF( IINFO.NE.0 ) THEN INFO = -1 RETURN ENDIF * Make sure that the entries M+1 to N in W, WERR, IBLOCK, INDEXW are 0 DO 14 I = MM+1,N W( I ) = ZERO WERR( I ) = ZERO IBLOCK( I ) = 0 INDEXW( I ) = 0 14 CONTINUE *** * Loop over unreduced blocks IBEGIN = 1 WBEGIN = 1 DO 170 JBLK = 1, NSPLIT IEND = ISPLIT( JBLK ) IN = IEND - IBEGIN + 1 * 1 X 1 block IF( IN.EQ.1 ) THEN IF( (IRANGE.EQ.1).OR.( (IRANGE.EQ.2).AND. $ ( D( IBEGIN ).GT.VL ).AND.( D( IBEGIN ).LE.VU ) ) $ .OR. ( (IRANGE.EQ.3).AND.(IBLOCK(WBEGIN).EQ.JBLK)) $ ) 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 WGAP(M) = ZERO IBLOCK( M ) = JBLK INDEXW( M ) = 1 WBEGIN = WBEGIN + 1 ENDIF * E( IEND ) holds the shift for the initial RRR E( IEND ) = ZERO IBEGIN = IEND + 1 GO TO 170 END IF * * Blocks of size larger than 1x1 * * E( IEND ) will hold the shift for the initial RRR, for now set it =0 E( IEND ) = ZERO IF( ( IRANGE.EQ.1 ) .OR. $ ((IRANGE.EQ.3).AND.(IL.EQ.1.AND.IU.EQ.N)) ) THEN * MB = number of eigenvalues to compute MB = IN WEND = WBEGIN + MB - 1 INDL = 1 INDU = IN ELSE * Count the number of eigenvalues in the current block. MB = 0 DO 20 I = WBEGIN,MM IF( IBLOCK(I).EQ.JBLK ) THEN MB = MB+1 ELSE GOTO 21 ENDIF 20 CONTINUE 21 CONTINUE IF( MB.EQ.0) THEN * No eigenvalue in the current block lies in the desired range * E( IEND ) holds the shift for the initial RRR E( IEND ) = ZERO IBEGIN = IEND + 1 GO TO 170 ENDIF * WEND = WBEGIN + MB - 1 * Find local index of the first and last desired evalue. INDL = INDEXW(WBEGIN) INDU = INDEXW( WEND ) ENDIF * IF( (WEND.LT.DOL).OR.(WBEGIN.GT.DOU) ) THEN * if this subblock contains no desired eigenvalues, * skip the computation of this representation tree IBEGIN = IEND + 1 WBEGIN = WEND + 1 M = M + MB GO TO 170 END IF * IF(.NOT. ( IRANGE.EQ.1 ) ) THEN * At this point, the sequential code decides * whether dqds or bisection is more efficient. * Note: in the parallel code, we do not use dqds. * However, we do not change the shift strategy * if USEDQD is TRUE, then the same shift is used as for * the sequential code when it uses dqds. * USEDQD = ( MB .GT. FAC*IN ) * * Calculate gaps for the current block * In later stages, when representations for individual * eigenvalues are different, we use SIGMA = E( IEND ). SIGMA = ZERO DO 30 I = WBEGIN, WEND - 1 WGAP( I ) = MAX( ZERO, $ W(I+1)-WERR(I+1) - (W(I)+WERR(I)) ) 30 CONTINUE WGAP( WEND ) = MAX( ZERO, $ VU - SIGMA - (W( WEND )+WERR( WEND ))) ENDIF * * Find local outer bounds GL,GU for the block GL = D(IBEGIN) GU = D(IBEGIN) DO 15 I = IBEGIN , IEND GL = MIN( GERS( 2*I-1 ), GL ) GU = MAX( GERS( 2*I ), GU ) 15 CONTINUE SPDIAM = GU - GL * Save local spectral diameter for later use SDIAM(JBLK) = SPDIAM * Find approximations to the extremal eigenvalues of the block CALL DLARRK( IN, 1, GL, GU, D(IBEGIN), $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) IF( IINFO.NE.0 ) THEN INFO = -1 RETURN ENDIF ISLEFT = MAX(GL, TMP - TMP1 $ - HNDRD * EPS* ABS(TMP - TMP1)) CALL DLARRK( IN, IN, GL, GU, D(IBEGIN), $ E2(IBEGIN), PIVMIN, RTL, TMP, TMP1, IINFO ) IF( IINFO.NE.0 ) THEN INFO = -1 RETURN ENDIF ISRGHT = MIN(GU, TMP + TMP1 $ + HNDRD * EPS * ABS(TMP + TMP1)) IF( ( IRANGE.EQ.1 ).OR.USEDQD ) THEN * Case of DQDS shift * Improve the estimate of the spectral diameter SPDIAM = ISRGHT - ISLEFT ELSE * Case of bisection * Find approximations to the wanted extremal eigenvalues ISLEFT = MAX(GL, W(WBEGIN) - WERR(WBEGIN) $ - HNDRD * EPS*ABS(W(WBEGIN)- WERR(WBEGIN) )) ISRGHT = MIN(GU,W(WEND) + WERR(WEND) $ + HNDRD * EPS * ABS(W(WEND)+ WERR(WEND))) ENDIF * Decide whether the base representation for the current block * L_JBLK D_JBLK L_JBLK^T = T_JBLK - sigma_JBLK I * should be on the left or the right end of the current block. * The strategy is to shift to the end which is "more populated" IF( IRANGE.EQ.1 ) THEN * If all the eigenvalues have to be computed, we use dqd USEDQD = .TRUE. * INDL is the local index of the first eigenvalue to compute INDL = 1 INDU = IN * MB = number of eigenvalues to compute MB = IN WEND = WBEGIN + MB - 1 * Define 1/4 and 3/4 points of the spectrum S1 = ISLEFT + FOURTH * SPDIAM S2 = ISRGHT - FOURTH * SPDIAM ELSE * DLARRD2 has computed IBLOCK and INDEXW for each eigenvalue * approximation. * choose sigma IF( USEDQD ) THEN S1 = ISLEFT + FOURTH * SPDIAM S2 = ISRGHT - FOURTH * SPDIAM ELSE TMP = MIN(ISRGHT,VU) - MAX(ISLEFT,VL) S1 = MAX(ISLEFT,VL) + FOURTH * TMP S2 = MIN(ISRGHT,VU) - FOURTH * TMP ENDIF ENDIF * Compute the negcount at the 1/4 and 3/4 points IF(MB.GT.2) THEN CALL DLARRC( 'T', IN, S1, S2, D(IBEGIN), $ E(IBEGIN), PIVMIN, CNT, CNT1, CNT2, IINFO) ENDIF IF(MB.LE.2) THEN SIGMA = GL SGNDEF = ONE ELSEIF( CNT1 - INDL .GE. INDU - CNT2 ) THEN IF( IRANGE.EQ.1 ) THEN SIGMA = MAX(ISLEFT,GL) ELSEIF( USEDQD ) THEN * use Gerschgorin bound as shift to get pos def matrix SIGMA = ISLEFT ELSE * use approximation of the first desired eigenvalue of the * block as shift SIGMA = MAX(ISLEFT,VL) ENDIF SGNDEF = ONE ELSE IF( IRANGE.EQ.1 ) THEN SIGMA = MIN(ISRGHT,GU) ELSEIF( USEDQD ) THEN * use Gerschgorin bound as shift to get neg def matrix * for dqds SIGMA = ISRGHT ELSE * use approximation of the first desired eigenvalue of the * block as shift SIGMA = MIN(ISRGHT,VU) ENDIF SGNDEF = -ONE ENDIF * An initial SIGMA has been chosen that will be used for computing * T - SIGMA I = L D L^T * Define the increment TAU of the shift in case the initial shift * needs to be refined to obtain a factorization with not too much * element growth. IF( USEDQD ) THEN TAU = SPDIAM*EPS*N + TWO*PIVMIN TAU = MAX(TAU,EPS*ABS(SIGMA)) ELSE IF(MB.GT.1) THEN CLWDTH = W(WEND) + WERR(WEND) - W(WBEGIN) - WERR(WBEGIN) AVGAP = ABS(CLWDTH / DBLE(WEND-WBEGIN)) IF( SGNDEF.EQ.ONE ) THEN TAU = HALF*MAX(WGAP(WBEGIN),AVGAP) TAU = MAX(TAU,WERR(WBEGIN)) ELSE TAU = HALF*MAX(WGAP(WEND-1),AVGAP) TAU = MAX(TAU,WERR(WEND)) ENDIF ELSE TAU = WERR(WBEGIN) ENDIF ENDIF * DO 80 IDUM = 1, MAXTRY * Compute L D L^T factorization of tridiagonal matrix T - sigma I. * Store D in WORK(1:IN), L in WORK(IN+1:2*IN), and reciprocals of * pivots in WORK(2*IN+1:3*IN) DPIVOT = D( IBEGIN ) - SIGMA WORK( 1 ) = DPIVOT DMAX = ABS( WORK(1) ) J = IBEGIN DO 70 I = 1, IN - 1 WORK( 2*IN+I ) = ONE / WORK( I ) TMP = E( J )*WORK( 2*IN+I ) WORK( IN+I ) = TMP DPIVOT = ( D( J+1 )-SIGMA ) - TMP*E( J ) WORK( I+1 ) = DPIVOT DMAX = MAX( DMAX, ABS(DPIVOT) ) J = J + 1 70 CONTINUE * check for element growth IF( DMAX .GT. MAXGROWTH*SPDIAM ) THEN NOREP = .TRUE. ELSE NOREP = .FALSE. ENDIF IF(NOREP) THEN * Note that in the case of IRANGE=1, we use the Gerschgorin * shift which makes the matrix definite. So we should end up * here really only in the case of IRANGE = 2,3 IF( IDUM.EQ.MAXTRY-1 ) THEN IF( SGNDEF.EQ.ONE ) THEN * The fudged Gerschgorin shift should succeed SIGMA = $ GL - FUDGE*SPDIAM*EPS*N - FUDGE*TWO*PIVMIN ELSE SIGMA = $ GU + FUDGE*SPDIAM*EPS*N + FUDGE*TWO*PIVMIN END IF ELSE SIGMA = SIGMA - SGNDEF * TAU TAU = TWO * TAU END IF ELSE * an initial RRR is found GO TO 83 END IF 80 CONTINUE * if the program reaches this point, no base representation could be * found in MAXTRY iterations. INFO = 2 RETURN 83 CONTINUE * At this point, we have found an initial base representation * T - SIGMA I = L D L^T with not too much element growth. * Store the shift. E( IEND ) = SIGMA * Store D and L. CALL DCOPY( IN, WORK, 1, D( IBEGIN ), 1 ) CALL DCOPY( IN-1, WORK( IN+1 ), 1, E( IBEGIN ), 1 ) IF(RNDPRT .AND. MB.GT.1 ) THEN * * Perturb each entry of the base representation by a small * (but random) relative amount to overcome difficulties with * glued matrices. * DO 122 I = 1, 4 ISEED( I ) = 1 122 CONTINUE CALL DLARNV(2, ISEED, 2*IN-1, WORK(1)) DO 125 I = 1,IN-1 D(IBEGIN+I-1) = D(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(2*I-1)) E(IBEGIN+I-1) = E(IBEGIN+I-1)*(ONE+EPS*PERT*WORK(2*I)) 125 CONTINUE D(IEND) = D(IEND)*(ONE+EPS*PERT*WORK(2*IN-1)) * ENDIF * * Compute the required eigenvalues of L D L' by bisection * Shift the eigenvalue approximations * according to the shift of their representation. DO 134 J=WBEGIN,WEND W(J) = W(J) - SIGMA WERR(J) = WERR(J) + ABS(W(J)) * EPS 134 CONTINUE * call DLARRB2 to reduce eigenvalue error of the approximations * from DLARRD2 DO 135 I = IBEGIN, IEND-1 WORK( I ) = D( I ) * E( I )**2 135 CONTINUE * use bisection to find EV from INDL to INDU INDL = INDEXW( WBEGIN ) INDU = INDEXW( WEND ) * * Indicate that the current block contains eigenvalues that * are potentially needed later. * NEEDIL = MIN(NEEDIL,WBEGIN) NEEDIU = MAX(NEEDIU,WEND) * * For the parallel distributed case, only compute * those eigenvalues that have to be computed as indicated by DOL, DOU * MYWBEG = MAX(WBEGIN,DOL) MYWEND = MIN(WEND,DOU) * IF(MYWBEG.GT.WBEGIN) THEN * This is the leftmost block containing wanted eigenvalues * as well as unwanted ones. To save on communication, * check if NEEDIL can be increased even further: * on the left end, only the eigenvalues of the cluster * including MYWBEG are needed DO 136 I = WBEGIN, MYWBEG-1 IF ( WGAP(I).GE.MINRGP*ABS(W(I)) ) THEN NEEDIL = MAX(I+1,NEEDIL) ENDIF 136 CONTINUE ENDIF IF(MYWEND.LT.WEND) THEN * This is the rightmost block containing wanted eigenvalues * as well as unwanted ones. To save on communication, * Check if NEEDIU can be decreased even further. DO 137 I = MYWEND,WEND-1 IF ( WGAP(I).GE.MINRGP*ABS(W(I)) ) THEN NEEDIU = MIN(I,NEEDIU) GOTO 138 ENDIF 137 CONTINUE 138 CONTINUE ENDIF * * Only compute eigenvalues from MYINDL to MYINDU * instead of INDL to INDU * MYINDL = INDEXW( MYWBEG ) MYINDU = INDEXW( MYWEND ) * LGPVMN = LOG( PIVMIN ) LGSPDM = LOG( SPDIAM + PIVMIN ) CALL DLARRB2(IN, D(IBEGIN), WORK(IBEGIN), $ MYINDL, MYINDU, RTOL1, RTOL2, MYINDL-1, $ W(MYWBEG), WGAP(MYWBEG), WERR(MYWBEG), $ WORK( 2*N+1 ), IWORK, PIVMIN, $ LGPVMN, LGSPDM, IN, IINFO ) IF( IINFO .NE. 0 ) THEN INFO = -4 RETURN END IF * DLARRB2 computes all gaps correctly except for the last one * Record distance to VU/GU WGAP( WEND ) = MAX( ZERO, $ ( VU-SIGMA ) - ( W( WEND ) + WERR( WEND ) ) ) DO 140 I = INDL, INDU M = M + 1 IBLOCK(M) = JBLK INDEXW(M) = I 140 CONTINUE * * proceed with next block IBEGIN = IEND + 1 WBEGIN = WEND + 1 170 CONTINUE * IF (M.LT.DOU-DOL+1) THEN INFO = -9 ENDIF RETURN * * end of DLARRE2A * END