001:       SUBROUTINE SLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
002:      $                   RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
003:      $                   PIVMIN, SPDIAM, TWIST, INFO )
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            IFIRST, ILAST, INFO, N, OFFSET, TWIST
011:       REAL               PIVMIN, RTOL1, RTOL2, SPDIAM
012: *     ..
013: *     .. Array Arguments ..
014:       INTEGER            IWORK( * )
015:       REAL               D( * ), LLD( * ), W( * ),
016:      $                   WERR( * ), WGAP( * ), WORK( * )
017: *     ..
018: *
019: *  Purpose
020: *  =======
021: *
022: *  Given the relatively robust representation(RRR) L D L^T, SLARRB
023: *  does "limited" bisection to refine the eigenvalues of L D L^T,
024: *  W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
025: *  guesses for these eigenvalues are input in W, the corresponding estimate
026: *  of the error in these guesses and their gaps are input in WERR
027: *  and WGAP, respectively. During bisection, intervals
028: *  [left, right] are maintained by storing their mid-points and
029: *  semi-widths in the arrays W and WERR respectively.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  N       (input) INTEGER
035: *          The order of the matrix.
036: *
037: *  D       (input) REAL             array, dimension (N)
038: *          The N diagonal elements of the diagonal matrix D.
039: *
040: *  LLD     (input) REAL             array, dimension (N-1)
041: *          The (N-1) elements L(i)*L(i)*D(i).
042: *
043: *  IFIRST  (input) INTEGER
044: *          The index of the first eigenvalue to be computed.
045: *
046: *  ILAST   (input) INTEGER
047: *          The index of the last eigenvalue to be computed.
048: *
049: *  RTOL1   (input) REAL            
050: *  RTOL2   (input) REAL            
051: *          Tolerance for the convergence of the bisection intervals.
052: *          An interval [LEFT,RIGHT] has converged if
053: *          RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
054: *          where GAP is the (estimated) distance to the nearest
055: *          eigenvalue.
056: *
057: *  OFFSET  (input) INTEGER
058: *          Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
059: *          through ILAST-OFFSET elements of these arrays are to be used.
060: *
061: *  W       (input/output) REAL             array, dimension (N)
062: *          On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
063: *          estimates of the eigenvalues of L D L^T indexed IFIRST throug
064: *          ILAST.
065: *          On output, these estimates are refined.
066: *
067: *  WGAP    (input/output) REAL             array, dimension (N-1)
068: *          On input, the (estimated) gaps between consecutive
069: *          eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
070: *          eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
071: *          then WGAP(IFIRST-OFFSET) must be set to ZERO.
072: *          On output, these gaps are refined.
073: *
074: *  WERR    (input/output) REAL             array, dimension (N)
075: *          On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
076: *          the errors in the estimates of the corresponding elements in W.
077: *          On output, these errors are refined.
078: *
079: *  WORK    (workspace) REAL             array, dimension (2*N)
080: *          Workspace.
081: *
082: *  IWORK   (workspace) INTEGER array, dimension (2*N)
083: *          Workspace.
084: *
085: *  PIVMIN  (input) DOUBLE PRECISION
086: *          The minimum pivot in the Sturm sequence.
087: *
088: *  SPDIAM  (input) DOUBLE PRECISION
089: *          The spectral diameter of the matrix.
090: *
091: *  TWIST   (input) INTEGER
092: *          The twist index for the twisted factorization that is used
093: *          for the negcount.
094: *          TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
095: *          TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
096: *          TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
097: *
098: *  INFO    (output) INTEGER
099: *          Error flag.
100: *
101: *  Further Details
102: *  ===============
103: *
104: *  Based on contributions by
105: *     Beresford Parlett, University of California, Berkeley, USA
106: *     Jim Demmel, University of California, Berkeley, USA
107: *     Inderjit Dhillon, University of Texas, Austin, USA
108: *     Osni Marques, LBNL/NERSC, USA
109: *     Christof Voemel, University of California, Berkeley, USA
110: *
111: *  =====================================================================
112: *
113: *     .. Parameters ..
114:       REAL               ZERO, TWO, HALF
115:       PARAMETER        ( ZERO = 0.0E0, TWO = 2.0E0,
116:      $                   HALF = 0.5E0 )
117:       INTEGER   MAXITR
118: *     ..
119: *     .. Local Scalars ..
120:       INTEGER            I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
121:      $                   OLNINT, PREV, R
122:       REAL               BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
123:      $                   RGAP, RIGHT, TMP, WIDTH
124: *     ..
125: *     .. External Functions ..
126:       INTEGER            SLANEG
127:       EXTERNAL           SLANEG
128: *
129: *     ..
130: *     .. Intrinsic Functions ..
131:       INTRINSIC          ABS, MAX, MIN
132: *     ..
133: *     .. Executable Statements ..
134: *
135:       INFO = 0
136: *
137:       MAXITR = INT( ( LOG( SPDIAM+PIVMIN )-LOG( PIVMIN ) ) /
138:      $           LOG( TWO ) ) + 2
139:       MNWDTH = TWO * PIVMIN
140: *
141:       R = TWIST
142:       IF((R.LT.1).OR.(R.GT.N)) R = N
143: *
144: *     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].
145: *     The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while
146: *     Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )
147: *     for an unconverged interval is set to the index of the next unconverged
148: *     interval, and is -1 or 0 for a converged interval. Thus a linked
149: *     list of unconverged intervals is set up.
150: *
151:       I1 = IFIRST
152: *     The number of unconverged intervals
153:       NINT = 0
154: *     The last unconverged interval found
155:       PREV = 0
156: 
157:       RGAP = WGAP( I1-OFFSET )
158:       DO 75 I = I1, ILAST
159:          K = 2*I
160:          II = I - OFFSET
161:          LEFT = W( II ) - WERR( II )
162:          RIGHT = W( II ) + WERR( II )
163:          LGAP = RGAP
164:          RGAP = WGAP( II )
165:          GAP = MIN( LGAP, RGAP )
166: 
167: *        Make sure that [LEFT,RIGHT] contains the desired eigenvalue
168: *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT
169: *
170: *        Do while( NEGCNT(LEFT).GT.I-1 )
171: *
172:          BACK = WERR( II )
173:  20      CONTINUE
174:          NEGCNT = SLANEG( N, D, LLD, LEFT, PIVMIN, R )
175:          IF( NEGCNT.GT.I-1 ) THEN
176:             LEFT = LEFT - BACK
177:             BACK = TWO*BACK
178:             GO TO 20
179:          END IF
180: *
181: *        Do while( NEGCNT(RIGHT).LT.I )
182: *        Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT
183: *
184:          BACK = WERR( II )
185:  50      CONTINUE
186: 
187:          NEGCNT = SLANEG( N, D, LLD, RIGHT, PIVMIN, R )
188:           IF( NEGCNT.LT.I ) THEN
189:              RIGHT = RIGHT + BACK
190:              BACK = TWO*BACK
191:              GO TO 50
192:           END IF
193:          WIDTH = HALF*ABS( LEFT - RIGHT )
194:          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
195:          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
196:          IF( WIDTH.LE.CVRGD .OR. WIDTH.LE.MNWDTH ) THEN
197: *           This interval has already converged and does not need refinement.
198: *           (Note that the gaps might change through refining the
199: *            eigenvalues, however, they can only get bigger.)
200: *           Remove it from the list.
201:             IWORK( K-1 ) = -1
202: *           Make sure that I1 always points to the first unconverged interval
203:             IF((I.EQ.I1).AND.(I.LT.ILAST)) I1 = I + 1
204:             IF((PREV.GE.I1).AND.(I.LE.ILAST)) IWORK( 2*PREV-1 ) = I + 1
205:          ELSE
206: *           unconverged interval found
207:             PREV = I
208:             NINT = NINT + 1
209:             IWORK( K-1 ) = I + 1
210:             IWORK( K ) = NEGCNT
211:          END IF
212:          WORK( K-1 ) = LEFT
213:          WORK( K ) = RIGHT
214:  75   CONTINUE
215: 
216: *
217: *     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals
218: *     and while (ITER.LT.MAXITR)
219: *
220:       ITER = 0
221:  80   CONTINUE
222:       PREV = I1 - 1
223:       I = I1
224:       OLNINT = NINT
225: 
226:       DO 100 IP = 1, OLNINT
227:          K = 2*I
228:          II = I - OFFSET
229:          RGAP = WGAP( II )
230:          LGAP = RGAP
231:          IF(II.GT.1) LGAP = WGAP( II-1 )
232:          GAP = MIN( LGAP, RGAP )
233:          NEXT = IWORK( K-1 )
234:          LEFT = WORK( K-1 )
235:          RIGHT = WORK( K )
236:          MID = HALF*( LEFT + RIGHT )
237: 
238: *        semiwidth of interval
239:          WIDTH = RIGHT - MID
240:          TMP = MAX( ABS( LEFT ), ABS( RIGHT ) )
241:          CVRGD = MAX(RTOL1*GAP,RTOL2*TMP)
242:          IF( ( WIDTH.LE.CVRGD ) .OR. ( WIDTH.LE.MNWDTH ).OR.
243:      $       ( ITER.EQ.MAXITR ) )THEN
244: *           reduce number of unconverged intervals
245:             NINT = NINT - 1
246: *           Mark interval as converged.
247:             IWORK( K-1 ) = 0
248:             IF( I1.EQ.I ) THEN
249:                I1 = NEXT
250:             ELSE
251: *              Prev holds the last unconverged interval previously examined
252:                IF(PREV.GE.I1) IWORK( 2*PREV-1 ) = NEXT
253:             END IF
254:             I = NEXT
255:             GO TO 100
256:          END IF
257:          PREV = I
258: *
259: *        Perform one bisection step
260: *
261:          NEGCNT = SLANEG( N, D, LLD, MID, PIVMIN, R )
262:          IF( NEGCNT.LE.I-1 ) THEN
263:             WORK( K-1 ) = MID
264:          ELSE
265:             WORK( K ) = MID
266:          END IF
267:          I = NEXT
268:  100  CONTINUE
269:       ITER = ITER + 1
270: *     do another loop if there are still unconverged intervals
271: *     However, in the last iteration, all intervals are accepted
272: *     since this is the best we can do.
273:       IF( ( NINT.GT.0 ).AND.(ITER.LE.MAXITR) ) GO TO 80
274: *
275: *
276: *     At this point, all the intervals have converged
277:       DO 110 I = IFIRST, ILAST
278:          K = 2*I
279:          II = I - OFFSET
280: *        All intervals marked by '0' have been refined.
281:          IF( IWORK( K-1 ).EQ.0 ) THEN
282:             W( II ) = HALF*( WORK( K-1 )+WORK( K ) )
283:             WERR( II ) = WORK( K ) - W( II )
284:          END IF
285:  110  CONTINUE
286: *
287:       DO 111 I = IFIRST+1, ILAST
288:          K = 2*I
289:          II = I - OFFSET
290:          WGAP( II-1 ) = MAX( ZERO,
291:      $                     W(II) - WERR (II) - W( II-1 ) - WERR( II-1 ))
292:  111  CONTINUE
293: 
294:       RETURN
295: *
296: *     End of SLARRB
297: *
298:       END
299: