Annotation of rpl/lapack/lapack/dlarrf.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLARRF( N, D, L, LD, CLSTRT, CLEND,
! 2: $ W, WGAP, WERR,
! 3: $ SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
! 4: $ DPLUS, LPLUS, WORK, INFO )
! 5: *
! 6: * -- LAPACK auxiliary routine (version 3.2) --
! 7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 9: * November 2006
! 10: **
! 11: * .. Scalar Arguments ..
! 12: INTEGER CLSTRT, CLEND, INFO, N
! 13: DOUBLE PRECISION CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
! 14: * ..
! 15: * .. Array Arguments ..
! 16: DOUBLE PRECISION D( * ), DPLUS( * ), L( * ), LD( * ),
! 17: $ LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * Given the initial representation L D L^T and its cluster of close
! 24: * eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
! 25: * W( CLEND ), DLARRF finds a new relatively robust representation
! 26: * L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
! 27: * eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
! 28: *
! 29: * Arguments
! 30: * =========
! 31: *
! 32: * N (input) INTEGER
! 33: * The order of the matrix (subblock, if the matrix splitted).
! 34: *
! 35: * D (input) DOUBLE PRECISION array, dimension (N)
! 36: * The N diagonal elements of the diagonal matrix D.
! 37: *
! 38: * L (input) DOUBLE PRECISION array, dimension (N-1)
! 39: * The (N-1) subdiagonal elements of the unit bidiagonal
! 40: * matrix L.
! 41: *
! 42: * LD (input) DOUBLE PRECISION array, dimension (N-1)
! 43: * The (N-1) elements L(i)*D(i).
! 44: *
! 45: * CLSTRT (input) INTEGER
! 46: * The index of the first eigenvalue in the cluster.
! 47: *
! 48: * CLEND (input) INTEGER
! 49: * The index of the last eigenvalue in the cluster.
! 50: *
! 51: * W (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1)
! 52: * The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
! 53: * W( CLSTRT ) through W( CLEND ) form the cluster of relatively
! 54: * close eigenalues.
! 55: *
! 56: * WGAP (input/output) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1)
! 57: * The separation from the right neighbor eigenvalue in W.
! 58: *
! 59: * WERR (input) DOUBLE PRECISION array, dimension >= (CLEND-CLSTRT+1)
! 60: * WERR contain the semiwidth of the uncertainty
! 61: * interval of the corresponding eigenvalue APPROXIMATION in W
! 62: *
! 63: * SPDIAM (input) estimate of the spectral diameter obtained from the
! 64: * Gerschgorin intervals
! 65: *
! 66: * CLGAPL, CLGAPR (input) absolute gap on each end of the cluster.
! 67: * Set by the calling routine to protect against shifts too close
! 68: * to eigenvalues outside the cluster.
! 69: *
! 70: * PIVMIN (input) DOUBLE PRECISION
! 71: * The minimum pivot allowed in the Sturm sequence.
! 72: *
! 73: * SIGMA (output) DOUBLE PRECISION
! 74: * The shift used to form L(+) D(+) L(+)^T.
! 75: *
! 76: * DPLUS (output) DOUBLE PRECISION array, dimension (N)
! 77: * The N diagonal elements of the diagonal matrix D(+).
! 78: *
! 79: * LPLUS (output) DOUBLE PRECISION array, dimension (N-1)
! 80: * The first (N-1) elements of LPLUS contain the subdiagonal
! 81: * elements of the unit bidiagonal matrix L(+).
! 82: *
! 83: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 84: * Workspace.
! 85: *
! 86: * Further Details
! 87: * ===============
! 88: *
! 89: * Based on contributions by
! 90: * Beresford Parlett, University of California, Berkeley, USA
! 91: * Jim Demmel, University of California, Berkeley, USA
! 92: * Inderjit Dhillon, University of Texas, Austin, USA
! 93: * Osni Marques, LBNL/NERSC, USA
! 94: * Christof Voemel, University of California, Berkeley, USA
! 95: *
! 96: * =====================================================================
! 97: *
! 98: * .. Parameters ..
! 99: DOUBLE PRECISION FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO,
! 100: $ ZERO
! 101: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
! 102: $ FOUR = 4.0D0, QUART = 0.25D0,
! 103: $ MAXGROWTH1 = 8.D0,
! 104: $ MAXGROWTH2 = 8.D0 )
! 105: * ..
! 106: * .. Local Scalars ..
! 107: LOGICAL DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
! 108: INTEGER I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
! 109: PARAMETER ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
! 110: DOUBLE PRECISION AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
! 111: $ FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
! 112: $ MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
! 113: $ RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
! 114: * ..
! 115: * .. External Functions ..
! 116: LOGICAL DISNAN
! 117: DOUBLE PRECISION DLAMCH
! 118: EXTERNAL DISNAN, DLAMCH
! 119: * ..
! 120: * .. External Subroutines ..
! 121: EXTERNAL DCOPY
! 122: * ..
! 123: * .. Intrinsic Functions ..
! 124: INTRINSIC ABS
! 125: * ..
! 126: * .. Executable Statements ..
! 127: *
! 128: INFO = 0
! 129: FACT = DBLE(2**KTRYMAX)
! 130: EPS = DLAMCH( 'Precision' )
! 131: SHIFT = 0
! 132: FORCER = .FALSE.
! 133:
! 134:
! 135: * Note that we cannot guarantee that for any of the shifts tried,
! 136: * the factorization has a small or even moderate element growth.
! 137: * There could be Ritz values at both ends of the cluster and despite
! 138: * backing off, there are examples where all factorizations tried
! 139: * (in IEEE mode, allowing zero pivots & infinities) have INFINITE
! 140: * element growth.
! 141: * For this reason, we should use PIVMIN in this subroutine so that at
! 142: * least the L D L^T factorization exists. It can be checked afterwards
! 143: * whether the element growth caused bad residuals/orthogonality.
! 144:
! 145: * Decide whether the code should accept the best among all
! 146: * representations despite large element growth or signal INFO=1
! 147: NOFAIL = .TRUE.
! 148: *
! 149:
! 150: * Compute the average gap length of the cluster
! 151: CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
! 152: AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
! 153: MINGAP = MIN(CLGAPL, CLGAPR)
! 154: * Initial values for shifts to both ends of cluster
! 155: LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
! 156: RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )
! 157:
! 158: * Use a small fudge to make sure that we really shift to the outside
! 159: LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
! 160: RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS
! 161:
! 162: * Compute upper bounds for how much to back off the initial shifts
! 163: LDMAX = QUART * MINGAP + TWO * PIVMIN
! 164: RDMAX = QUART * MINGAP + TWO * PIVMIN
! 165:
! 166: LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
! 167: RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
! 168: *
! 169: * Initialize the record of the best representation found
! 170: *
! 171: S = DLAMCH( 'S' )
! 172: SMLGROWTH = ONE / S
! 173: FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
! 174: FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
! 175: BESTSHIFT = LSIGMA
! 176: *
! 177: * while (KTRY <= KTRYMAX)
! 178: KTRY = 0
! 179: GROWTHBOUND = MAXGROWTH1*SPDIAM
! 180:
! 181: 5 CONTINUE
! 182: SAWNAN1 = .FALSE.
! 183: SAWNAN2 = .FALSE.
! 184: * Ensure that we do not back off too much of the initial shifts
! 185: LDELTA = MIN(LDMAX,LDELTA)
! 186: RDELTA = MIN(RDMAX,RDELTA)
! 187:
! 188: * Compute the element growth when shifting to both ends of the cluster
! 189: * accept the shift if there is no element growth at one of the two ends
! 190:
! 191: * Left end
! 192: S = -LSIGMA
! 193: DPLUS( 1 ) = D( 1 ) + S
! 194: IF(ABS(DPLUS(1)).LT.PIVMIN) THEN
! 195: DPLUS(1) = -PIVMIN
! 196: * Need to set SAWNAN1 because refined RRR test should not be used
! 197: * in this case
! 198: SAWNAN1 = .TRUE.
! 199: ENDIF
! 200: MAX1 = ABS( DPLUS( 1 ) )
! 201: DO 6 I = 1, N - 1
! 202: LPLUS( I ) = LD( I ) / DPLUS( I )
! 203: S = S*LPLUS( I )*L( I ) - LSIGMA
! 204: DPLUS( I+1 ) = D( I+1 ) + S
! 205: IF(ABS(DPLUS(I+1)).LT.PIVMIN) THEN
! 206: DPLUS(I+1) = -PIVMIN
! 207: * Need to set SAWNAN1 because refined RRR test should not be used
! 208: * in this case
! 209: SAWNAN1 = .TRUE.
! 210: ENDIF
! 211: MAX1 = MAX( MAX1,ABS(DPLUS(I+1)) )
! 212: 6 CONTINUE
! 213: SAWNAN1 = SAWNAN1 .OR. DISNAN( MAX1 )
! 214:
! 215: IF( FORCER .OR.
! 216: $ (MAX1.LE.GROWTHBOUND .AND. .NOT.SAWNAN1 ) ) THEN
! 217: SIGMA = LSIGMA
! 218: SHIFT = SLEFT
! 219: GOTO 100
! 220: ENDIF
! 221:
! 222: * Right end
! 223: S = -RSIGMA
! 224: WORK( 1 ) = D( 1 ) + S
! 225: IF(ABS(WORK(1)).LT.PIVMIN) THEN
! 226: WORK(1) = -PIVMIN
! 227: * Need to set SAWNAN2 because refined RRR test should not be used
! 228: * in this case
! 229: SAWNAN2 = .TRUE.
! 230: ENDIF
! 231: MAX2 = ABS( WORK( 1 ) )
! 232: DO 7 I = 1, N - 1
! 233: WORK( N+I ) = LD( I ) / WORK( I )
! 234: S = S*WORK( N+I )*L( I ) - RSIGMA
! 235: WORK( I+1 ) = D( I+1 ) + S
! 236: IF(ABS(WORK(I+1)).LT.PIVMIN) THEN
! 237: WORK(I+1) = -PIVMIN
! 238: * Need to set SAWNAN2 because refined RRR test should not be used
! 239: * in this case
! 240: SAWNAN2 = .TRUE.
! 241: ENDIF
! 242: MAX2 = MAX( MAX2,ABS(WORK(I+1)) )
! 243: 7 CONTINUE
! 244: SAWNAN2 = SAWNAN2 .OR. DISNAN( MAX2 )
! 245:
! 246: IF( FORCER .OR.
! 247: $ (MAX2.LE.GROWTHBOUND .AND. .NOT.SAWNAN2 ) ) THEN
! 248: SIGMA = RSIGMA
! 249: SHIFT = SRIGHT
! 250: GOTO 100
! 251: ENDIF
! 252: * If we are at this point, both shifts led to too much element growth
! 253:
! 254: * Record the better of the two shifts (provided it didn't lead to NaN)
! 255: IF(SAWNAN1.AND.SAWNAN2) THEN
! 256: * both MAX1 and MAX2 are NaN
! 257: GOTO 50
! 258: ELSE
! 259: IF( .NOT.SAWNAN1 ) THEN
! 260: INDX = 1
! 261: IF(MAX1.LE.SMLGROWTH) THEN
! 262: SMLGROWTH = MAX1
! 263: BESTSHIFT = LSIGMA
! 264: ENDIF
! 265: ENDIF
! 266: IF( .NOT.SAWNAN2 ) THEN
! 267: IF(SAWNAN1 .OR. MAX2.LE.MAX1) INDX = 2
! 268: IF(MAX2.LE.SMLGROWTH) THEN
! 269: SMLGROWTH = MAX2
! 270: BESTSHIFT = RSIGMA
! 271: ENDIF
! 272: ENDIF
! 273: ENDIF
! 274:
! 275: * If we are here, both the left and the right shift led to
! 276: * element growth. If the element growth is moderate, then
! 277: * we may still accept the representation, if it passes a
! 278: * refined test for RRR. This test supposes that no NaN occurred.
! 279: * Moreover, we use the refined RRR test only for isolated clusters.
! 280: IF((CLWDTH.LT.MINGAP/DBLE(128)) .AND.
! 281: $ (MIN(MAX1,MAX2).LT.FAIL2)
! 282: $ .AND.(.NOT.SAWNAN1).AND.(.NOT.SAWNAN2)) THEN
! 283: DORRR1 = .TRUE.
! 284: ELSE
! 285: DORRR1 = .FALSE.
! 286: ENDIF
! 287: TRYRRR1 = .TRUE.
! 288: IF( TRYRRR1 .AND. DORRR1 ) THEN
! 289: IF(INDX.EQ.1) THEN
! 290: TMP = ABS( DPLUS( N ) )
! 291: ZNM2 = ONE
! 292: PROD = ONE
! 293: OLDP = ONE
! 294: DO 15 I = N-1, 1, -1
! 295: IF( PROD .LE. EPS ) THEN
! 296: PROD =
! 297: $ ((DPLUS(I+1)*WORK(N+I+1))/(DPLUS(I)*WORK(N+I)))*OLDP
! 298: ELSE
! 299: PROD = PROD*ABS(WORK(N+I))
! 300: END IF
! 301: OLDP = PROD
! 302: ZNM2 = ZNM2 + PROD**2
! 303: TMP = MAX( TMP, ABS( DPLUS( I ) * PROD ))
! 304: 15 CONTINUE
! 305: RRR1 = TMP/( SPDIAM * SQRT( ZNM2 ) )
! 306: IF (RRR1.LE.MAXGROWTH2) THEN
! 307: SIGMA = LSIGMA
! 308: SHIFT = SLEFT
! 309: GOTO 100
! 310: ENDIF
! 311: ELSE IF(INDX.EQ.2) THEN
! 312: TMP = ABS( WORK( N ) )
! 313: ZNM2 = ONE
! 314: PROD = ONE
! 315: OLDP = ONE
! 316: DO 16 I = N-1, 1, -1
! 317: IF( PROD .LE. EPS ) THEN
! 318: PROD = ((WORK(I+1)*LPLUS(I+1))/(WORK(I)*LPLUS(I)))*OLDP
! 319: ELSE
! 320: PROD = PROD*ABS(LPLUS(I))
! 321: END IF
! 322: OLDP = PROD
! 323: ZNM2 = ZNM2 + PROD**2
! 324: TMP = MAX( TMP, ABS( WORK( I ) * PROD ))
! 325: 16 CONTINUE
! 326: RRR2 = TMP/( SPDIAM * SQRT( ZNM2 ) )
! 327: IF (RRR2.LE.MAXGROWTH2) THEN
! 328: SIGMA = RSIGMA
! 329: SHIFT = SRIGHT
! 330: GOTO 100
! 331: ENDIF
! 332: END IF
! 333: ENDIF
! 334:
! 335: 50 CONTINUE
! 336:
! 337: IF (KTRY.LT.KTRYMAX) THEN
! 338: * If we are here, both shifts failed also the RRR test.
! 339: * Back off to the outside
! 340: LSIGMA = MAX( LSIGMA - LDELTA,
! 341: $ LSIGMA - LDMAX)
! 342: RSIGMA = MIN( RSIGMA + RDELTA,
! 343: $ RSIGMA + RDMAX )
! 344: LDELTA = TWO * LDELTA
! 345: RDELTA = TWO * RDELTA
! 346: KTRY = KTRY + 1
! 347: GOTO 5
! 348: ELSE
! 349: * None of the representations investigated satisfied our
! 350: * criteria. Take the best one we found.
! 351: IF((SMLGROWTH.LT.FAIL).OR.NOFAIL) THEN
! 352: LSIGMA = BESTSHIFT
! 353: RSIGMA = BESTSHIFT
! 354: FORCER = .TRUE.
! 355: GOTO 5
! 356: ELSE
! 357: INFO = 1
! 358: RETURN
! 359: ENDIF
! 360: END IF
! 361:
! 362: 100 CONTINUE
! 363: IF (SHIFT.EQ.SLEFT) THEN
! 364: ELSEIF (SHIFT.EQ.SRIGHT) THEN
! 365: * store new L and D back into DPLUS, LPLUS
! 366: CALL DCOPY( N, WORK, 1, DPLUS, 1 )
! 367: CALL DCOPY( N-1, WORK(N+1), 1, LPLUS, 1 )
! 368: ENDIF
! 369:
! 370: RETURN
! 371: *
! 372: * End of DLARRF
! 373: *
! 374: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlarrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack