187 SUBROUTINE slarrf( N, D, L, LD, CLSTRT, CLEND,
189 $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
190 $ DPLUS, LPLUS, WORK, INFO )
197 INTEGER CLSTRT, CLEND, INFO, N
198 REAL CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
201 REAL D( * ), DPLUS( * ), L( * ), LD( * ),
202 $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
208 REAL MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO
209 PARAMETER ( ONE = 1.0e0, two = 2.0e0,
212 $ maxgrowth2 = 8.e0 )
215 LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
216 INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
217 PARAMETER ( KTRYMAX = 1, sleft = 1, sright = 2 )
218 REAL AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
219 $ fail2, growthbound, ldelta, ldmax, lsigma,
220 $ max1, max2, mingap, oldp, prod, rdelta, rdmax,
221 $ rrr1, rrr2, rsigma, s, smlgrowth, tmp, znm2
226 EXTERNAL SISNAN, SLAMCH
244 fact = real(2**ktrymax)
245 eps = slamch(
'Precision' )
267 clwdth = abs(w(clend)-w(clstrt)) + werr(clend) + werr(clstrt)
268 avgap = clwdth / real(clend-clstrt)
269 mingap = min(clgapl, clgapr)
271 lsigma = min(w( clstrt ),w( clend )) - werr( clstrt )
272 rsigma = max(w( clstrt ),w( clend )) + werr( clend )
275 lsigma = lsigma - abs(lsigma)* two * eps
276 rsigma = rsigma + abs(rsigma)* two * eps
279 ldmax = quart * mingap + two * pivmin
280 rdmax = quart * mingap + two * pivmin
282 ldelta = max(avgap,wgap( clstrt ))/fact
283 rdelta = max(avgap,wgap( clend-1 ))/fact
289 fail = real(n-1)*mingap/(spdiam*eps)
290 fail2 = real(n-1)*mingap/(spdiam*sqrt(eps))
295 growthbound = maxgrowth1*spdiam
301 ldelta = min(ldmax,ldelta)
302 rdelta = min(rdmax,rdelta)
309 dplus( 1 ) = d( 1 ) + s
310 IF(abs(dplus(1)).LT.pivmin)
THEN
316 max1 = abs( dplus( 1 ) )
318 lplus( i ) = ld( i ) / dplus( i )
319 s = s*lplus( i )*l( i ) - lsigma
320 dplus( i+1 ) = d( i+1 ) + s
321 IF(abs(dplus(i+1)).LT.pivmin)
THEN
327 max1 = max( max1,abs(dplus(i+1)) )
329 sawnan1 = sawnan1 .OR. sisnan( max1 )
332 $ (max1.LE.growthbound .AND. .NOT.sawnan1 ) )
THEN
340 work( 1 ) = d( 1 ) + s
341 IF(abs(work(1)).LT.pivmin)
THEN
347 max2 = abs( work( 1 ) )
349 work( n+i ) = ld( i ) / work( i )
350 s = s*work( n+i )*l( i ) - rsigma
351 work( i+1 ) = d( i+1 ) + s
352 IF(abs(work(i+1)).LT.pivmin)
THEN
358 max2 = max( max2,abs(work(i+1)) )
360 sawnan2 = sawnan2 .OR. sisnan( max2 )
363 $ (max2.LE.growthbound .AND. .NOT.sawnan2 ) )
THEN
371 IF(sawnan1.AND.sawnan2)
THEN
375 IF( .NOT.sawnan1 )
THEN
377 IF(max1.LE.smlgrowth)
THEN
382 IF( .NOT.sawnan2 )
THEN
383 IF(sawnan1 .OR. max2.LE.max1) indx = 2
384 IF(max2.LE.smlgrowth)
THEN
396 IF((clwdth.LT.mingap/real(128)) .AND.
397 $ (min(max1,max2).LT.fail2)
398 $ .AND.(.NOT.sawnan1).AND.(.NOT.sawnan2))
THEN
404 IF( tryrrr1 .AND. dorrr1 )
THEN
406 tmp = abs( dplus( n ) )
411 IF( prod .LE. eps )
THEN
413 $ ((dplus(i+1)*work(n+i+1))/(dplus(i)*work(n+i)))*oldp
415 prod = prod*abs(work(n+i))
418 znm2 = znm2 + prod**2
419 tmp = max( tmp, abs( dplus( i ) * prod ))
421 rrr1 = tmp/( spdiam * sqrt( znm2 ) )
422 IF (rrr1.LE.maxgrowth2)
THEN
427 ELSE IF(indx.EQ.2)
THEN
428 tmp = abs( work( n ) )
433 IF( prod .LE. eps )
THEN
434 prod = ((work(i+1)*lplus(i+1))/(work(i)*lplus(i)))*oldp
436 prod = prod*abs(lplus(i))
439 znm2 = znm2 + prod**2
440 tmp = max( tmp, abs( work( i ) * prod ))
442 rrr2 = tmp/( spdiam * sqrt( znm2 ) )
443 IF (rrr2.LE.maxgrowth2)
THEN
453 IF (ktry.LT.ktrymax)
THEN
456 lsigma = max( lsigma - ldelta,
458 rsigma = min( rsigma + rdelta,
460 ldelta = two * ldelta
461 rdelta = two * rdelta
467 IF((smlgrowth.LT.fail).OR.nofail)
THEN
479 IF (shift.EQ.sleft)
THEN
480 ELSEIF (shift.EQ.sright)
THEN
482 CALL scopy( n, work, 1, dplus, 1 )
483 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...