SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1, $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK, $ PIVMIN, SPDIAM, TWIST, INFO ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM * .. * .. Array Arguments .. INTEGER IWORK( * ) DOUBLE PRECISION D( * ), LLD( * ), W( * ), $ WERR( * ), WGAP( * ), WORK( * ) * .. * * Purpose * ======= * * Given the relatively robust representation(RRR) L D L^T, DLARRB * does "limited" bisection to refine the eigenvalues of L D L^T, * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial * guesses for these eigenvalues are input in W, the corresponding estimate * of the error in these guesses and their gaps are input in WERR * and WGAP, respectively. During bisection, intervals * [left, right] are maintained by storing their mid-points and * semi-widths in the arrays W and WERR respectively. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix. * * D (input) DOUBLE PRECISION array, dimension (N) * The N diagonal elements of the diagonal matrix D. * * LLD (input) DOUBLE PRECISION array, dimension (N-1) * The (N-1) elements L(i)*L(i)*D(i). * * IFIRST (input) INTEGER * The index of the first eigenvalue to be computed. * * ILAST (input) INTEGER * The index of the last eigenvalue to be computed. * * RTOL1 (input) DOUBLE PRECISION * RTOL2 (input) DOUBLE PRECISION * Tolerance for the convergence of the bisection intervals. * An interval [LEFT,RIGHT] has converged if * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) * where GAP is the (estimated) distance to the nearest * eigenvalue. * * OFFSET (input) INTEGER * Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET * through ILAST-OFFSET elements of these arrays are to be used. * * W (input/output) DOUBLE PRECISION array, dimension (N) * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are * estimates of the eigenvalues of L D L^T indexed IFIRST throug * ILAST. * On output, these estimates are refined. * * WGAP (input/output) DOUBLE PRECISION array, dimension (N-1) * On input, the (estimated) gaps between consecutive * eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between * eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST * then WGAP(IFIRST-OFFSET) must be set to ZERO. * On output, these gaps are refined. * * WERR (input/output) DOUBLE PRECISION array, dimension (N) * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are * the errors in the estimates of the corresponding elements in W. * On output, these errors are refined. * * WORK (workspace) DOUBLE PRECISION array, dimension (2*N) * Workspace. * * IWORK (workspace) INTEGER array, dimension (2*N) * Workspace. * * PIVMIN (input) DOUBLE PRECISION * The minimum pivot in the Sturm sequence. * * SPDIAM (input) DOUBLE PRECISION * The spectral diameter of the matrix. * * TWIST (input) INTEGER * The twist index for the twisted factorization that is used * for the negcount. * TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T * TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T * TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r) * * INFO (output) INTEGER * Error flag. * * Further Details * =============== * * Based on contributions by * 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, TWO, HALF PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0, $ HALF = 0.5D0 ) INTEGER MAXITR * .. * .. Local Scalars .. INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT, $ OLNINT, PREV, R DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH, $ RGAP, RIGHT, TMP, WIDTH * .. * .. External Functions .. INTEGER DLANEG EXTERNAL DLANEG * * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Executable Statements .. * INFO = 0 * MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) / $ LOG( TWO ) ) + 2 MNWDTH = TWO * PIVMIN * R = TWIST IF((R.LT.1).OR.(R.GT.N)) R = N * * Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ]. * The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while * Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 ) * for an unconverged interval is set to the index of the next unconverged * interval, and is -1 or 0 for a converged interval. Thus a linked * list of unconverged intervals is set up. * I1 = IFIRST * The number of unconverged intervals NINT = 0 * The last unconverged interval found PREV = 0 RGAP = WGAP( I1-OFFSET ) DO 75 I = I1, ILAST K = 2*I II = I - OFFSET LEFT = W( II ) - WERR( II ) RIGHT = W( II ) + WERR( II ) LGAP = RGAP RGAP = WGAP( II ) GAP = MIN( LGAP, RGAP ) * Make sure that [LEFT,RIGHT] contains the desired eigenvalue * Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT * * Do while( NEGCNT(LEFT).GT.I-1 ) * BACK = WERR( II ) 20 CONTINUE NEGCNT = DLANEG( N, D, LLD, LEFT, PIVMIN, R ) IF( NEGCNT.GT.I-1 ) THEN LEFT = LEFT - BACK BACK = TWO*BACK GO TO 20 END IF * * Do while( NEGCNT(RIGHT).LT.I ) * Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT * BACK = WERR( II ) 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 WIDTH = HALF*ABS( LEFT - RIGHT ) TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) CVRGD = MAX(RTOL1*GAP,RTOL2*TMP) IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN * This interval has already converged and does not need refinement. * (Note that the gaps might change through refining the * eigenvalues, however, they can only get bigger.) * Remove it from the list. IWORK( K-1 ) = -1 * Make sure that I1 always points to the first unconverged interval IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1 IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1 ELSE * unconverged interval found PREV = I NINT = NINT + 1 IWORK( K-1 ) = I + 1 IWORK( K ) = NEGCNT END IF WORK( K-1 ) = LEFT WORK( K ) = RIGHT 75 CONTINUE * * Do while( NINT.GT.0 ), i.e. there are still unconverged intervals * and while (ITER.LT.MAXITR) * ITER = 0 80 CONTINUE PREV = I1 - 1 I = I1 OLNINT = NINT DO 100 IP = 1, OLNINT K = 2*I II = I - OFFSET RGAP = WGAP( II ) LGAP = RGAP IF(II.GT.1) LGAP = WGAP( II-1 ) GAP = MIN( LGAP, RGAP ) NEXT = IWORK( K-1 ) LEFT = WORK( K-1 ) RIGHT = WORK( K ) MID = HALF*( LEFT + RIGHT ) * semiwidth of interval WIDTH = RIGHT - MID TMP = MAX( ABS( LEFT ), ABS( RIGHT ) ) CVRGD = MAX(RTOL1*GAP,RTOL2*TMP) IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR. $ ( ITER.EQ.MAXITR ) )THEN * reduce number of unconverged intervals NINT = NINT - 1 * Mark interval as converged. IWORK( K-1 ) = 0 IF( I1.EQ.I ) THEN I1 = NEXT ELSE * Prev holds the last unconverged interval previously examined IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT END IF I = NEXT GO TO 100 END IF PREV = I * * Perform one bisection step * NEGCNT = DLANEG( N, D, LLD, MID, PIVMIN, R ) IF( NEGCNT.LE.I-1 ) THEN WORK( K-1 ) = MID ELSE WORK( K ) = MID END IF I = NEXT 100 CONTINUE ITER = ITER + 1 * do another loop if there are still unconverged intervals * However, in the last iteration, all intervals are accepted * since this is the best we can do. IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80 * * * At this point, all the intervals have converged DO 110 I = IFIRST, ILAST K = 2*I II = I - OFFSET * All intervals marked by '0' have been refined. IF( IWORK( K-1 ).EQ.0 ) THEN W( II ) = HALF*( WORK( K-1 )+WORK( K ) ) WERR( II ) = WORK( K ) - W( II ) END IF 110 CONTINUE * DO 111 I = IFIRST+1, ILAST K = 2*I II = I - OFFSET WGAP( II-1 ) = MAX( ZERO, $ W(II) - WERR (II) - W( II-1 ) - WERR( II-1 )) 111 CONTINUE RETURN * * End of DLARRB * END