191 SUBROUTINE dlarrf( N, D, L, LD, CLSTRT, CLEND,
193 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
194 $ DPLUS, LPLUS, WORK, INFO )
202 INTEGER clstrt, clend, info, n
203 DOUBLE PRECISION clgapl, clgapr, pivmin, sigma, spdiam
206 DOUBLE PRECISION d( * ), dplus( * ), l( * ), ld( * ),
207 $ lplus( * ), w( * ), wgap( * ), werr( * ), work( * )
213 DOUBLE PRECISION four, maxgrowth1, maxgrowth2, one, quart, two,
215 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
216 $ four = 4.0d0, quart = 0.25d0,
218 $ maxgrowth2 = 8.d0 )
221 LOGICAL dorrr1, forcer, nofail, sawnan1, sawnan2, tryrrr1
222 INTEGER i, indx, ktry, ktrymax, sleft, sright, shift
223 parameter( ktrymax = 1, sleft = 1, sright = 2 )
224 DOUBLE PRECISION avgap, bestshift, clwdth, eps, fact, fail,
225 $ fail2, growthbound, ldelta, ldmax, lsigma,
226 $ max1, max2, mingap, oldp, prod, rdelta, rdmax,
227 $ rrr1, rrr2, rsigma, s, smlgrowth, tmp, znm2
231 DOUBLE PRECISION dlamch
232 EXTERNAL disnan, dlamch
243 fact = dble(2**ktrymax)
244 eps = dlamch(
'Precision' )
265 clwdth = abs(w(clend)-w(clstrt)) + werr(clend) + werr(clstrt)
266 avgap = clwdth / dble(clend-clstrt)
267 mingap = min(clgapl, clgapr)
269 lsigma = min(w( clstrt ),w( clend )) - werr( clstrt )
270 rsigma = max(w( clstrt ),w( clend )) + werr( clend )
273 lsigma = lsigma - abs(lsigma)* four * eps
274 rsigma = rsigma + abs(rsigma)* four * eps
277 ldmax = quart * mingap + two * pivmin
278 rdmax = quart * mingap + two * pivmin
280 ldelta = max(avgap,wgap( clstrt ))/fact
281 rdelta = max(avgap,wgap( clend-1 ))/fact
287 fail = dble(n-1)*mingap/(spdiam*eps)
288 fail2 = dble(n-1)*mingap/(spdiam*sqrt(eps))
293 growthbound = maxgrowth1*spdiam
299 ldelta = min(ldmax,ldelta)
300 rdelta = min(rdmax,rdelta)
307 dplus( 1 ) = d( 1 ) + s
308 IF(abs(dplus(1)).LT.pivmin)
THEN 314 max1 = abs( dplus( 1 ) )
316 lplus( i ) = ld( i ) / dplus( i )
317 s = s*lplus( i )*l( i ) - lsigma
318 dplus( i+1 ) = d( i+1 ) + s
319 IF(abs(dplus(i+1)).LT.pivmin)
THEN 325 max1 = max( max1,abs(dplus(i+1)) )
327 sawnan1 = sawnan1 .OR. disnan( max1 )
330 $ (max1.LE.growthbound .AND. .NOT.sawnan1 ) )
THEN 338 work( 1 ) = d( 1 ) + s
339 IF(abs(work(1)).LT.pivmin)
THEN 345 max2 = abs( work( 1 ) )
347 work( n+i ) = ld( i ) / work( i )
348 s = s*work( n+i )*l( i ) - rsigma
349 work( i+1 ) = d( i+1 ) + s
350 IF(abs(work(i+1)).LT.pivmin)
THEN 356 max2 = max( max2,abs(work(i+1)) )
358 sawnan2 = sawnan2 .OR. disnan( max2 )
361 $ (max2.LE.growthbound .AND. .NOT.sawnan2 ) )
THEN 369 IF(sawnan1.AND.sawnan2)
THEN 373 IF( .NOT.sawnan1 )
THEN 375 IF(max1.LE.smlgrowth)
THEN 380 IF( .NOT.sawnan2 )
THEN 381 IF(sawnan1 .OR. max2.LE.max1) indx = 2
382 IF(max2.LE.smlgrowth)
THEN 394 IF((clwdth.LT.mingap/dble(128)) .AND.
395 $ (min(max1,max2).LT.fail2)
396 $ .AND.(.NOT.sawnan1).AND.(.NOT.sawnan2))
THEN 402 IF( tryrrr1 .AND. dorrr1 )
THEN 404 tmp = abs( dplus( n ) )
409 IF( prod .LE. eps )
THEN 411 $ ((dplus(i+1)*work(n+i+1))/(dplus(i)*work(n+i)))*oldp
413 prod = prod*abs(work(n+i))
416 znm2 = znm2 + prod**2
417 tmp = max( tmp, abs( dplus( i ) * prod ))
419 rrr1 = tmp/( spdiam * sqrt( znm2 ) )
420 IF (rrr1.LE.maxgrowth2)
THEN 425 ELSE IF(indx.EQ.2)
THEN 426 tmp = abs( work( n ) )
431 IF( prod .LE. eps )
THEN 432 prod = ((work(i+1)*lplus(i+1))/(work(i)*lplus(i)))*oldp
434 prod = prod*abs(lplus(i))
437 znm2 = znm2 + prod**2
438 tmp = max( tmp, abs( work( i ) * prod ))
440 rrr2 = tmp/( spdiam * sqrt( znm2 ) )
441 IF (rrr2.LE.maxgrowth2)
THEN 451 IF (ktry.LT.ktrymax)
THEN 454 lsigma = max( lsigma - ldelta,
456 rsigma = min( rsigma + rdelta,
458 ldelta = two * ldelta
459 rdelta = two * rdelta
465 IF((smlgrowth.LT.fail).OR.nofail)
THEN 477 IF (shift.EQ.sleft)
THEN 478 ELSEIF (shift.EQ.sright)
THEN 480 CALL dcopy( n, work, 1, dplus, 1 )
481 CALL dcopy( n-1, work(n+1), 1, lplus, 1 )