1 SUBROUTINE DLARRB( N, D, LLD, IFIRST, ILAST, RTOL1,
2 $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
3 $ PIVMIN, SPDIAM, TWIST, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
12 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * )
16 DOUBLE PRECISION D( * ), LLD( * ), W( * ),
17 $ WERR( * ), WGAP( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * Given the relatively robust representation(RRR) L D L^T, DLARRB
24 * does "limited" bisection to refine the eigenvalues of L D L^T,
25 * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
26 * guesses for these eigenvalues are input in W, the corresponding estimate
27 * of the error in these guesses and their gaps are input in WERR
28 * and WGAP, respectively. During bisection, intervals
29 * [left, right] are maintained by storing their mid-points and
30 * semi-widths in the arrays W and WERR respectively.
31 *
32 * Arguments
33 * =========
34 *
35 * N (input) INTEGER
36 * The order of the matrix.
37 *
38 * D (input) DOUBLE PRECISION array, dimension (N)
39 * The N diagonal elements of the diagonal matrix D.
40 *
41 * LLD (input) DOUBLE PRECISION array, dimension (N-1)
42 * The (N-1) elements L(i)*L(i)*D(i).
43 *
44 * IFIRST (input) INTEGER
45 * The index of the first eigenvalue to be computed.
46 *
47 * ILAST (input) INTEGER
48 * The index of the last eigenvalue to be computed.
49 *
50 * RTOL1 (input) DOUBLE PRECISION
51 * RTOL2 (input) DOUBLE PRECISION
52 * Tolerance for the convergence of the bisection intervals.
53 * An interval [LEFT,RIGHT] has converged if
54 * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
55 * where GAP is the (estimated) distance to the nearest
56 * eigenvalue.
57 *
58 * OFFSET (input) INTEGER
59 * Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
60 * through ILAST-OFFSET elements of these arrays are to be used.
61 *
62 * W (input/output) DOUBLE PRECISION array, dimension (N)
63 * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
64 * estimates of the eigenvalues of L D L^T indexed IFIRST throug
65 * ILAST.
66 * On output, these estimates are refined.
67 *
68 * WGAP (input/output) DOUBLE PRECISION array, dimension (N-1)
69 * On input, the (estimated) gaps between consecutive
70 * eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
71 * eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
72 * then WGAP(IFIRST-OFFSET) must be set to ZERO.
73 * On output, these gaps are refined.
74 *
75 * WERR (input/output) DOUBLE PRECISION array, dimension (N)
76 * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
77 * the errors in the estimates of the corresponding elements in W.
78 * On output, these errors are refined.
79 *
80 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
81 * Workspace.
82 *
83 * IWORK (workspace) INTEGER array, dimension (2*N)
84 * Workspace.
85 *
86 * PIVMIN (input) DOUBLE PRECISION
87 * The minimum pivot in the Sturm sequence.
88 *
89 * SPDIAM (input) DOUBLE PRECISION
90 * The spectral diameter of the matrix.
91 *
92 * TWIST (input) INTEGER
93 * The twist index for the twisted factorization that is used
94 * for the negcount.
95 * TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
96 * TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
97 * TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
98 *
99 * 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 DOUBLE PRECISION ZERO, TWO, HALF
116 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
117 $ HALF = 0.5D0 )
118 INTEGER MAXITR
119 * ..
120 * .. Local Scalars ..
121 INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
122 $ OLNINT, PREV, R
123 DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
124 $ RGAP, RIGHT, TMP, WIDTH
125 * ..
126 * .. External Functions ..
127 INTEGER DLANEG
128 EXTERNAL DLANEG
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 = DLANEG( 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 = DLANEG( 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 = DLANEG( 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 DLARRB
298 *
299 END
2 $ RTOL2, OFFSET, W, WGAP, WERR, WORK, IWORK,
3 $ PIVMIN, SPDIAM, TWIST, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.2) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER IFIRST, ILAST, INFO, N, OFFSET, TWIST
12 DOUBLE PRECISION PIVMIN, RTOL1, RTOL2, SPDIAM
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * )
16 DOUBLE PRECISION D( * ), LLD( * ), W( * ),
17 $ WERR( * ), WGAP( * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * Given the relatively robust representation(RRR) L D L^T, DLARRB
24 * does "limited" bisection to refine the eigenvalues of L D L^T,
25 * W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial
26 * guesses for these eigenvalues are input in W, the corresponding estimate
27 * of the error in these guesses and their gaps are input in WERR
28 * and WGAP, respectively. During bisection, intervals
29 * [left, right] are maintained by storing their mid-points and
30 * semi-widths in the arrays W and WERR respectively.
31 *
32 * Arguments
33 * =========
34 *
35 * N (input) INTEGER
36 * The order of the matrix.
37 *
38 * D (input) DOUBLE PRECISION array, dimension (N)
39 * The N diagonal elements of the diagonal matrix D.
40 *
41 * LLD (input) DOUBLE PRECISION array, dimension (N-1)
42 * The (N-1) elements L(i)*L(i)*D(i).
43 *
44 * IFIRST (input) INTEGER
45 * The index of the first eigenvalue to be computed.
46 *
47 * ILAST (input) INTEGER
48 * The index of the last eigenvalue to be computed.
49 *
50 * RTOL1 (input) DOUBLE PRECISION
51 * RTOL2 (input) DOUBLE PRECISION
52 * Tolerance for the convergence of the bisection intervals.
53 * An interval [LEFT,RIGHT] has converged if
54 * RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )
55 * where GAP is the (estimated) distance to the nearest
56 * eigenvalue.
57 *
58 * OFFSET (input) INTEGER
59 * Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET
60 * through ILAST-OFFSET elements of these arrays are to be used.
61 *
62 * W (input/output) DOUBLE PRECISION array, dimension (N)
63 * On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are
64 * estimates of the eigenvalues of L D L^T indexed IFIRST throug
65 * ILAST.
66 * On output, these estimates are refined.
67 *
68 * WGAP (input/output) DOUBLE PRECISION array, dimension (N-1)
69 * On input, the (estimated) gaps between consecutive
70 * eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between
71 * eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST
72 * then WGAP(IFIRST-OFFSET) must be set to ZERO.
73 * On output, these gaps are refined.
74 *
75 * WERR (input/output) DOUBLE PRECISION array, dimension (N)
76 * On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are
77 * the errors in the estimates of the corresponding elements in W.
78 * On output, these errors are refined.
79 *
80 * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
81 * Workspace.
82 *
83 * IWORK (workspace) INTEGER array, dimension (2*N)
84 * Workspace.
85 *
86 * PIVMIN (input) DOUBLE PRECISION
87 * The minimum pivot in the Sturm sequence.
88 *
89 * SPDIAM (input) DOUBLE PRECISION
90 * The spectral diameter of the matrix.
91 *
92 * TWIST (input) INTEGER
93 * The twist index for the twisted factorization that is used
94 * for the negcount.
95 * TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T
96 * TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T
97 * TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)
98 *
99 * 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 DOUBLE PRECISION ZERO, TWO, HALF
116 PARAMETER ( ZERO = 0.0D0, TWO = 2.0D0,
117 $ HALF = 0.5D0 )
118 INTEGER MAXITR
119 * ..
120 * .. Local Scalars ..
121 INTEGER I, I1, II, IP, ITER, K, NEGCNT, NEXT, NINT,
122 $ OLNINT, PREV, R
123 DOUBLE PRECISION BACK, CVRGD, GAP, LEFT, LGAP, MID, MNWDTH,
124 $ RGAP, RIGHT, TMP, WIDTH
125 * ..
126 * .. External Functions ..
127 INTEGER DLANEG
128 EXTERNAL DLANEG
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 = DLANEG( 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 = DLANEG( 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 = DLANEG( 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 DLARRB
298 *
299 END