191 SUBROUTINE slarrf( N, D, L, LD, CLSTRT, CLEND,
193 $ spdiam, clgapl, clgapr, pivmin, sigma,
194 $ dplus, lplus, work, info )
202 INTEGER clstrt, clend, info, n
203 REAL clgapl, clgapr, pivmin, sigma, spdiam
206 REAL d( * ), dplus( * ), l( * ), ld( * ),
207 $ lplus( * ), w( * ), wgap( * ), werr( * ), work( * )
213 REAL maxgrowth1, maxgrowth2, one, quart, two
214 parameter( one = 1.0e0, two = 2.0e0,
217 $ maxgrowth2 = 8.e0 )
220 LOGICAL dorrr1, forcer, nofail, sawnan1, sawnan2, tryrrr1
221 INTEGER i, indx, ktry, ktrymax, sleft, sright, shift
222 parameter( ktrymax = 1, sleft = 1, sright = 2 )
223 REAL avgap, bestshift, clwdth, eps, fact, fail,
224 $ fail2, growthbound, ldelta, ldmax, lsigma,
225 $ max1, max2, mingap, oldp, prod, rdelta, rdmax,
226 $ rrr1, rrr2, rsigma, s, smlgrowth, tmp, znm2
242 fact =
REAL(2**ktrymax)
243 eps =
slamch(
'Precision' )
264 clwdth = abs(w(clend)-w(clstrt)) + werr(clend) + werr(clstrt)
265 avgap = clwdth /
REAL(clend-clstrt)
266 mingap = min(clgapl, clgapr)
268 lsigma = min(w( clstrt ),w( clend )) - werr( clstrt )
269 rsigma = max(w( clstrt ),w( clend )) + werr( clend )
272 lsigma = lsigma - abs(lsigma)* two * eps
273 rsigma = rsigma + abs(rsigma)* two * eps
276 ldmax = quart * mingap + two * pivmin
277 rdmax = quart * mingap + two * pivmin
279 ldelta = max(avgap,wgap( clstrt ))/fact
280 rdelta = max(avgap,wgap( clend-1 ))/fact
286 fail =
REAL(n-1)*mingap/(spdiam*eps)
287 fail2 =
REAL(n-1)*mingap/(spdiam*sqrt(eps))
292 growthbound = maxgrowth1*spdiam
298 ldelta = min(ldmax,ldelta)
299 rdelta = min(rdmax,rdelta)
306 dplus( 1 ) = d( 1 ) + s
307 IF(abs(dplus(1)).LT.pivmin)
THEN
313 max1 = abs( dplus( 1 ) )
315 lplus( i ) = ld( i ) / dplus( i )
316 s = s*lplus( i )*l( i ) - lsigma
317 dplus( i+1 ) = d( i+1 ) + s
318 IF(abs(dplus(i+1)).LT.pivmin)
THEN
324 max1 = max( max1,abs(dplus(i+1)) )
326 sawnan1 = sawnan1 .OR.
sisnan( max1 )
329 $ (max1.LE.growthbound .AND. .NOT.sawnan1 ) )
THEN
337 work( 1 ) = d( 1 ) + s
338 IF(abs(work(1)).LT.pivmin)
THEN
344 max2 = abs( work( 1 ) )
346 work( n+i ) = ld( i ) / work( i )
347 s = s*work( n+i )*l( i ) - rsigma
348 work( i+1 ) = d( i+1 ) + s
349 IF(abs(work(i+1)).LT.pivmin)
THEN
355 max2 = max( max2,abs(work(i+1)) )
357 sawnan2 = sawnan2 .OR.
sisnan( max2 )
360 $ (max2.LE.growthbound .AND. .NOT.sawnan2 ) )
THEN
368 IF(sawnan1.AND.sawnan2)
THEN
372 IF( .NOT.sawnan1 )
THEN
374 IF(max1.LE.smlgrowth)
THEN
379 IF( .NOT.sawnan2 )
THEN
380 IF(sawnan1 .OR. max2.LE.max1) indx = 2
381 IF(max2.LE.smlgrowth)
THEN
393 IF((clwdth.LT.mingap/
REAL(128)) .AND.
394 $ (min(max1,max2).LT.fail2)
395 $ .AND.(.NOT.sawnan1).AND.(.NOT.sawnan2))
THEN
401 IF( tryrrr1 .AND. dorrr1 )
THEN
403 tmp = abs( dplus( n ) )
408 IF( prod .LE. eps )
THEN
410 $ ((dplus(i+1)*work(n+i+1))/(dplus(i)*work(n+i)))*oldp
412 prod = prod*abs(work(n+i))
415 znm2 = znm2 + prod**2
416 tmp = max( tmp, abs( dplus( i ) * prod ))
418 rrr1 = tmp/( spdiam * sqrt( znm2 ) )
419 IF (rrr1.LE.maxgrowth2)
THEN
424 ELSE IF(indx.EQ.2)
THEN
425 tmp = abs( work( n ) )
430 IF( prod .LE. eps )
THEN
431 prod = ((work(i+1)*lplus(i+1))/(work(i)*lplus(i)))*oldp
433 prod = prod*abs(lplus(i))
436 znm2 = znm2 + prod**2
437 tmp = max( tmp, abs( work( i ) * prod ))
439 rrr2 = tmp/( spdiam * sqrt( znm2 ) )
440 IF (rrr2.LE.maxgrowth2)
THEN
450 IF (ktry.LT.ktrymax)
THEN
453 lsigma = max( lsigma - ldelta,
455 rsigma = min( rsigma + rdelta,
457 ldelta = two * ldelta
458 rdelta = two * rdelta
464 IF((smlgrowth.LT.fail).OR.nofail)
THEN
476 IF (shift.EQ.sleft)
THEN
477 ELSEIF (shift.EQ.sright)
THEN
479 CALL
scopy( n, work, 1, dplus, 1 )
480 CALL
scopy( n-1, work(n+1), 1, lplus, 1 )
subroutine slarrf(N, D, L, LD, CLSTRT, CLEND, W, WGAP, WERR, SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA, DPLUS, LPLUS, WORK, INFO)
SLARRF finds a new relatively robust representation such that at least one of the eigenvalues is rela...
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
real function slamch(CMACH)
SLAMCH
logical function sisnan(SIN)
SISNAN tests input for NaN.