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