Annotation of rpl/lapack/lapack/dlaed4.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
! 2: *
! 3: * -- LAPACK routine (version 3.2) --
! 4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 6: * November 2006
! 7: *
! 8: * .. Scalar Arguments ..
! 9: INTEGER I, INFO, N
! 10: DOUBLE PRECISION DLAM, RHO
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
! 14: * ..
! 15: *
! 16: * Purpose
! 17: * =======
! 18: *
! 19: * This subroutine computes the I-th updated eigenvalue of a symmetric
! 20: * rank-one modification to a diagonal matrix whose elements are
! 21: * given in the array d, and that
! 22: *
! 23: * D(i) < D(j) for i < j
! 24: *
! 25: * and that RHO > 0. This is arranged by the calling routine, and is
! 26: * no loss in generality. The rank-one modified system is thus
! 27: *
! 28: * diag( D ) + RHO * Z * Z_transpose.
! 29: *
! 30: * where we assume the Euclidean norm of Z is 1.
! 31: *
! 32: * The method consists of approximating the rational functions in the
! 33: * secular equation by simpler interpolating rational functions.
! 34: *
! 35: * Arguments
! 36: * =========
! 37: *
! 38: * N (input) INTEGER
! 39: * The length of all arrays.
! 40: *
! 41: * I (input) INTEGER
! 42: * The index of the eigenvalue to be computed. 1 <= I <= N.
! 43: *
! 44: * D (input) DOUBLE PRECISION array, dimension (N)
! 45: * The original eigenvalues. It is assumed that they are in
! 46: * order, D(I) < D(J) for I < J.
! 47: *
! 48: * Z (input) DOUBLE PRECISION array, dimension (N)
! 49: * The components of the updating vector.
! 50: *
! 51: * DELTA (output) DOUBLE PRECISION array, dimension (N)
! 52: * If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th
! 53: * component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
! 54: * for detail. The vector DELTA contains the information necessary
! 55: * to construct the eigenvectors by DLAED3 and DLAED9.
! 56: *
! 57: * RHO (input) DOUBLE PRECISION
! 58: * The scalar in the symmetric updating formula.
! 59: *
! 60: * DLAM (output) DOUBLE PRECISION
! 61: * The computed lambda_I, the I-th updated eigenvalue.
! 62: *
! 63: * INFO (output) INTEGER
! 64: * = 0: successful exit
! 65: * > 0: if INFO = 1, the updating process failed.
! 66: *
! 67: * Internal Parameters
! 68: * ===================
! 69: *
! 70: * Logical variable ORGATI (origin-at-i?) is used for distinguishing
! 71: * whether D(i) or D(i+1) is treated as the origin.
! 72: *
! 73: * ORGATI = .true. origin at i
! 74: * ORGATI = .false. origin at i+1
! 75: *
! 76: * Logical variable SWTCH3 (switch-for-3-poles?) is for noting
! 77: * if we are working with THREE poles!
! 78: *
! 79: * MAXIT is the maximum number of iterations allowed for each
! 80: * eigenvalue.
! 81: *
! 82: * Further Details
! 83: * ===============
! 84: *
! 85: * Based on contributions by
! 86: * Ren-Cang Li, Computer Science Division, University of California
! 87: * at Berkeley, USA
! 88: *
! 89: * =====================================================================
! 90: *
! 91: * .. Parameters ..
! 92: INTEGER MAXIT
! 93: PARAMETER ( MAXIT = 30 )
! 94: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
! 95: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
! 96: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0,
! 97: $ TEN = 10.0D0 )
! 98: * ..
! 99: * .. Local Scalars ..
! 100: LOGICAL ORGATI, SWTCH, SWTCH3
! 101: INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
! 102: DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
! 103: $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
! 104: $ RHOINV, TAU, TEMP, TEMP1, W
! 105: * ..
! 106: * .. Local Arrays ..
! 107: DOUBLE PRECISION ZZ( 3 )
! 108: * ..
! 109: * .. External Functions ..
! 110: DOUBLE PRECISION DLAMCH
! 111: EXTERNAL DLAMCH
! 112: * ..
! 113: * .. External Subroutines ..
! 114: EXTERNAL DLAED5, DLAED6
! 115: * ..
! 116: * .. Intrinsic Functions ..
! 117: INTRINSIC ABS, MAX, MIN, SQRT
! 118: * ..
! 119: * .. Executable Statements ..
! 120: *
! 121: * Since this routine is called in an inner loop, we do no argument
! 122: * checking.
! 123: *
! 124: * Quick return for N=1 and 2.
! 125: *
! 126: INFO = 0
! 127: IF( N.EQ.1 ) THEN
! 128: *
! 129: * Presumably, I=1 upon entry
! 130: *
! 131: DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 )
! 132: DELTA( 1 ) = ONE
! 133: RETURN
! 134: END IF
! 135: IF( N.EQ.2 ) THEN
! 136: CALL DLAED5( I, D, Z, DELTA, RHO, DLAM )
! 137: RETURN
! 138: END IF
! 139: *
! 140: * Compute machine epsilon
! 141: *
! 142: EPS = DLAMCH( 'Epsilon' )
! 143: RHOINV = ONE / RHO
! 144: *
! 145: * The case I = N
! 146: *
! 147: IF( I.EQ.N ) THEN
! 148: *
! 149: * Initialize some basic variables
! 150: *
! 151: II = N - 1
! 152: NITER = 1
! 153: *
! 154: * Calculate initial guess
! 155: *
! 156: MIDPT = RHO / TWO
! 157: *
! 158: * If ||Z||_2 is not one, then TEMP should be set to
! 159: * RHO * ||Z||_2^2 / TWO
! 160: *
! 161: DO 10 J = 1, N
! 162: DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
! 163: 10 CONTINUE
! 164: *
! 165: PSI = ZERO
! 166: DO 20 J = 1, N - 2
! 167: PSI = PSI + Z( J )*Z( J ) / DELTA( J )
! 168: 20 CONTINUE
! 169: *
! 170: C = RHOINV + PSI
! 171: W = C + Z( II )*Z( II ) / DELTA( II ) +
! 172: $ Z( N )*Z( N ) / DELTA( N )
! 173: *
! 174: IF( W.LE.ZERO ) THEN
! 175: TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) +
! 176: $ Z( N )*Z( N ) / RHO
! 177: IF( C.LE.TEMP ) THEN
! 178: TAU = RHO
! 179: ELSE
! 180: DEL = D( N ) - D( N-1 )
! 181: A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
! 182: B = Z( N )*Z( N )*DEL
! 183: IF( A.LT.ZERO ) THEN
! 184: TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
! 185: ELSE
! 186: TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
! 187: END IF
! 188: END IF
! 189: *
! 190: * It can be proved that
! 191: * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
! 192: *
! 193: DLTLB = MIDPT
! 194: DLTUB = RHO
! 195: ELSE
! 196: DEL = D( N ) - D( N-1 )
! 197: A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
! 198: B = Z( N )*Z( N )*DEL
! 199: IF( A.LT.ZERO ) THEN
! 200: TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
! 201: ELSE
! 202: TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
! 203: END IF
! 204: *
! 205: * It can be proved that
! 206: * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
! 207: *
! 208: DLTLB = ZERO
! 209: DLTUB = MIDPT
! 210: END IF
! 211: *
! 212: DO 30 J = 1, N
! 213: DELTA( J ) = ( D( J )-D( I ) ) - TAU
! 214: 30 CONTINUE
! 215: *
! 216: * Evaluate PSI and the derivative DPSI
! 217: *
! 218: DPSI = ZERO
! 219: PSI = ZERO
! 220: ERRETM = ZERO
! 221: DO 40 J = 1, II
! 222: TEMP = Z( J ) / DELTA( J )
! 223: PSI = PSI + Z( J )*TEMP
! 224: DPSI = DPSI + TEMP*TEMP
! 225: ERRETM = ERRETM + PSI
! 226: 40 CONTINUE
! 227: ERRETM = ABS( ERRETM )
! 228: *
! 229: * Evaluate PHI and the derivative DPHI
! 230: *
! 231: TEMP = Z( N ) / DELTA( N )
! 232: PHI = Z( N )*TEMP
! 233: DPHI = TEMP*TEMP
! 234: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
! 235: $ ABS( TAU )*( DPSI+DPHI )
! 236: *
! 237: W = RHOINV + PHI + PSI
! 238: *
! 239: * Test for convergence
! 240: *
! 241: IF( ABS( W ).LE.EPS*ERRETM ) THEN
! 242: DLAM = D( I ) + TAU
! 243: GO TO 250
! 244: END IF
! 245: *
! 246: IF( W.LE.ZERO ) THEN
! 247: DLTLB = MAX( DLTLB, TAU )
! 248: ELSE
! 249: DLTUB = MIN( DLTUB, TAU )
! 250: END IF
! 251: *
! 252: * Calculate the new step
! 253: *
! 254: NITER = NITER + 1
! 255: C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
! 256: A = ( DELTA( N-1 )+DELTA( N ) )*W -
! 257: $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
! 258: B = DELTA( N-1 )*DELTA( N )*W
! 259: IF( C.LT.ZERO )
! 260: $ C = ABS( C )
! 261: IF( C.EQ.ZERO ) THEN
! 262: * ETA = B/A
! 263: * ETA = RHO - TAU
! 264: ETA = DLTUB - TAU
! 265: ELSE IF( A.GE.ZERO ) THEN
! 266: ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
! 267: ELSE
! 268: ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
! 269: END IF
! 270: *
! 271: * Note, eta should be positive if w is negative, and
! 272: * eta should be negative otherwise. However,
! 273: * if for some reason caused by roundoff, eta*w > 0,
! 274: * we simply use one Newton step instead. This way
! 275: * will guarantee eta*w < 0.
! 276: *
! 277: IF( W*ETA.GT.ZERO )
! 278: $ ETA = -W / ( DPSI+DPHI )
! 279: TEMP = TAU + ETA
! 280: IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
! 281: IF( W.LT.ZERO ) THEN
! 282: ETA = ( DLTUB-TAU ) / TWO
! 283: ELSE
! 284: ETA = ( DLTLB-TAU ) / TWO
! 285: END IF
! 286: END IF
! 287: DO 50 J = 1, N
! 288: DELTA( J ) = DELTA( J ) - ETA
! 289: 50 CONTINUE
! 290: *
! 291: TAU = TAU + ETA
! 292: *
! 293: * Evaluate PSI and the derivative DPSI
! 294: *
! 295: DPSI = ZERO
! 296: PSI = ZERO
! 297: ERRETM = ZERO
! 298: DO 60 J = 1, II
! 299: TEMP = Z( J ) / DELTA( J )
! 300: PSI = PSI + Z( J )*TEMP
! 301: DPSI = DPSI + TEMP*TEMP
! 302: ERRETM = ERRETM + PSI
! 303: 60 CONTINUE
! 304: ERRETM = ABS( ERRETM )
! 305: *
! 306: * Evaluate PHI and the derivative DPHI
! 307: *
! 308: TEMP = Z( N ) / DELTA( N )
! 309: PHI = Z( N )*TEMP
! 310: DPHI = TEMP*TEMP
! 311: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
! 312: $ ABS( TAU )*( DPSI+DPHI )
! 313: *
! 314: W = RHOINV + PHI + PSI
! 315: *
! 316: * Main loop to update the values of the array DELTA
! 317: *
! 318: ITER = NITER + 1
! 319: *
! 320: DO 90 NITER = ITER, MAXIT
! 321: *
! 322: * Test for convergence
! 323: *
! 324: IF( ABS( W ).LE.EPS*ERRETM ) THEN
! 325: DLAM = D( I ) + TAU
! 326: GO TO 250
! 327: END IF
! 328: *
! 329: IF( W.LE.ZERO ) THEN
! 330: DLTLB = MAX( DLTLB, TAU )
! 331: ELSE
! 332: DLTUB = MIN( DLTUB, TAU )
! 333: END IF
! 334: *
! 335: * Calculate the new step
! 336: *
! 337: C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
! 338: A = ( DELTA( N-1 )+DELTA( N ) )*W -
! 339: $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
! 340: B = DELTA( N-1 )*DELTA( N )*W
! 341: IF( A.GE.ZERO ) THEN
! 342: ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
! 343: ELSE
! 344: ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
! 345: END IF
! 346: *
! 347: * Note, eta should be positive if w is negative, and
! 348: * eta should be negative otherwise. However,
! 349: * if for some reason caused by roundoff, eta*w > 0,
! 350: * we simply use one Newton step instead. This way
! 351: * will guarantee eta*w < 0.
! 352: *
! 353: IF( W*ETA.GT.ZERO )
! 354: $ ETA = -W / ( DPSI+DPHI )
! 355: TEMP = TAU + ETA
! 356: IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
! 357: IF( W.LT.ZERO ) THEN
! 358: ETA = ( DLTUB-TAU ) / TWO
! 359: ELSE
! 360: ETA = ( DLTLB-TAU ) / TWO
! 361: END IF
! 362: END IF
! 363: DO 70 J = 1, N
! 364: DELTA( J ) = DELTA( J ) - ETA
! 365: 70 CONTINUE
! 366: *
! 367: TAU = TAU + ETA
! 368: *
! 369: * Evaluate PSI and the derivative DPSI
! 370: *
! 371: DPSI = ZERO
! 372: PSI = ZERO
! 373: ERRETM = ZERO
! 374: DO 80 J = 1, II
! 375: TEMP = Z( J ) / DELTA( J )
! 376: PSI = PSI + Z( J )*TEMP
! 377: DPSI = DPSI + TEMP*TEMP
! 378: ERRETM = ERRETM + PSI
! 379: 80 CONTINUE
! 380: ERRETM = ABS( ERRETM )
! 381: *
! 382: * Evaluate PHI and the derivative DPHI
! 383: *
! 384: TEMP = Z( N ) / DELTA( N )
! 385: PHI = Z( N )*TEMP
! 386: DPHI = TEMP*TEMP
! 387: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
! 388: $ ABS( TAU )*( DPSI+DPHI )
! 389: *
! 390: W = RHOINV + PHI + PSI
! 391: 90 CONTINUE
! 392: *
! 393: * Return with INFO = 1, NITER = MAXIT and not converged
! 394: *
! 395: INFO = 1
! 396: DLAM = D( I ) + TAU
! 397: GO TO 250
! 398: *
! 399: * End for the case I = N
! 400: *
! 401: ELSE
! 402: *
! 403: * The case for I < N
! 404: *
! 405: NITER = 1
! 406: IP1 = I + 1
! 407: *
! 408: * Calculate initial guess
! 409: *
! 410: DEL = D( IP1 ) - D( I )
! 411: MIDPT = DEL / TWO
! 412: DO 100 J = 1, N
! 413: DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
! 414: 100 CONTINUE
! 415: *
! 416: PSI = ZERO
! 417: DO 110 J = 1, I - 1
! 418: PSI = PSI + Z( J )*Z( J ) / DELTA( J )
! 419: 110 CONTINUE
! 420: *
! 421: PHI = ZERO
! 422: DO 120 J = N, I + 2, -1
! 423: PHI = PHI + Z( J )*Z( J ) / DELTA( J )
! 424: 120 CONTINUE
! 425: C = RHOINV + PSI + PHI
! 426: W = C + Z( I )*Z( I ) / DELTA( I ) +
! 427: $ Z( IP1 )*Z( IP1 ) / DELTA( IP1 )
! 428: *
! 429: IF( W.GT.ZERO ) THEN
! 430: *
! 431: * d(i)< the ith eigenvalue < (d(i)+d(i+1))/2
! 432: *
! 433: * We choose d(i) as origin.
! 434: *
! 435: ORGATI = .TRUE.
! 436: A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
! 437: B = Z( I )*Z( I )*DEL
! 438: IF( A.GT.ZERO ) THEN
! 439: TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
! 440: ELSE
! 441: TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
! 442: END IF
! 443: DLTLB = ZERO
! 444: DLTUB = MIDPT
! 445: ELSE
! 446: *
! 447: * (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)
! 448: *
! 449: * We choose d(i+1) as origin.
! 450: *
! 451: ORGATI = .FALSE.
! 452: A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
! 453: B = Z( IP1 )*Z( IP1 )*DEL
! 454: IF( A.LT.ZERO ) THEN
! 455: TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
! 456: ELSE
! 457: TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
! 458: END IF
! 459: DLTLB = -MIDPT
! 460: DLTUB = ZERO
! 461: END IF
! 462: *
! 463: IF( ORGATI ) THEN
! 464: DO 130 J = 1, N
! 465: DELTA( J ) = ( D( J )-D( I ) ) - TAU
! 466: 130 CONTINUE
! 467: ELSE
! 468: DO 140 J = 1, N
! 469: DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU
! 470: 140 CONTINUE
! 471: END IF
! 472: IF( ORGATI ) THEN
! 473: II = I
! 474: ELSE
! 475: II = I + 1
! 476: END IF
! 477: IIM1 = II - 1
! 478: IIP1 = II + 1
! 479: *
! 480: * Evaluate PSI and the derivative DPSI
! 481: *
! 482: DPSI = ZERO
! 483: PSI = ZERO
! 484: ERRETM = ZERO
! 485: DO 150 J = 1, IIM1
! 486: TEMP = Z( J ) / DELTA( J )
! 487: PSI = PSI + Z( J )*TEMP
! 488: DPSI = DPSI + TEMP*TEMP
! 489: ERRETM = ERRETM + PSI
! 490: 150 CONTINUE
! 491: ERRETM = ABS( ERRETM )
! 492: *
! 493: * Evaluate PHI and the derivative DPHI
! 494: *
! 495: DPHI = ZERO
! 496: PHI = ZERO
! 497: DO 160 J = N, IIP1, -1
! 498: TEMP = Z( J ) / DELTA( J )
! 499: PHI = PHI + Z( J )*TEMP
! 500: DPHI = DPHI + TEMP*TEMP
! 501: ERRETM = ERRETM + PHI
! 502: 160 CONTINUE
! 503: *
! 504: W = RHOINV + PHI + PSI
! 505: *
! 506: * W is the value of the secular function with
! 507: * its ii-th element removed.
! 508: *
! 509: SWTCH3 = .FALSE.
! 510: IF( ORGATI ) THEN
! 511: IF( W.LT.ZERO )
! 512: $ SWTCH3 = .TRUE.
! 513: ELSE
! 514: IF( W.GT.ZERO )
! 515: $ SWTCH3 = .TRUE.
! 516: END IF
! 517: IF( II.EQ.1 .OR. II.EQ.N )
! 518: $ SWTCH3 = .FALSE.
! 519: *
! 520: TEMP = Z( II ) / DELTA( II )
! 521: DW = DPSI + DPHI + TEMP*TEMP
! 522: TEMP = Z( II )*TEMP
! 523: W = W + TEMP
! 524: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
! 525: $ THREE*ABS( TEMP ) + ABS( TAU )*DW
! 526: *
! 527: * Test for convergence
! 528: *
! 529: IF( ABS( W ).LE.EPS*ERRETM ) THEN
! 530: IF( ORGATI ) THEN
! 531: DLAM = D( I ) + TAU
! 532: ELSE
! 533: DLAM = D( IP1 ) + TAU
! 534: END IF
! 535: GO TO 250
! 536: END IF
! 537: *
! 538: IF( W.LE.ZERO ) THEN
! 539: DLTLB = MAX( DLTLB, TAU )
! 540: ELSE
! 541: DLTUB = MIN( DLTUB, TAU )
! 542: END IF
! 543: *
! 544: * Calculate the new step
! 545: *
! 546: NITER = NITER + 1
! 547: IF( .NOT.SWTCH3 ) THEN
! 548: IF( ORGATI ) THEN
! 549: C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )*
! 550: $ ( Z( I ) / DELTA( I ) )**2
! 551: ELSE
! 552: C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
! 553: $ ( Z( IP1 ) / DELTA( IP1 ) )**2
! 554: END IF
! 555: A = ( DELTA( I )+DELTA( IP1 ) )*W -
! 556: $ DELTA( I )*DELTA( IP1 )*DW
! 557: B = DELTA( I )*DELTA( IP1 )*W
! 558: IF( C.EQ.ZERO ) THEN
! 559: IF( A.EQ.ZERO ) THEN
! 560: IF( ORGATI ) THEN
! 561: A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )*
! 562: $ ( DPSI+DPHI )
! 563: ELSE
! 564: A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )*
! 565: $ ( DPSI+DPHI )
! 566: END IF
! 567: END IF
! 568: ETA = B / A
! 569: ELSE IF( A.LE.ZERO ) THEN
! 570: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
! 571: ELSE
! 572: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
! 573: END IF
! 574: ELSE
! 575: *
! 576: * Interpolation using THREE most relevant poles
! 577: *
! 578: TEMP = RHOINV + PSI + PHI
! 579: IF( ORGATI ) THEN
! 580: TEMP1 = Z( IIM1 ) / DELTA( IIM1 )
! 581: TEMP1 = TEMP1*TEMP1
! 582: C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) -
! 583: $ ( D( IIM1 )-D( IIP1 ) )*TEMP1
! 584: ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
! 585: ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*
! 586: $ ( ( DPSI-TEMP1 )+DPHI )
! 587: ELSE
! 588: TEMP1 = Z( IIP1 ) / DELTA( IIP1 )
! 589: TEMP1 = TEMP1*TEMP1
! 590: C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) -
! 591: $ ( D( IIP1 )-D( IIM1 ) )*TEMP1
! 592: ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*
! 593: $ ( DPSI+( DPHI-TEMP1 ) )
! 594: ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
! 595: END IF
! 596: ZZ( 2 ) = Z( II )*Z( II )
! 597: CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
! 598: $ INFO )
! 599: IF( INFO.NE.0 )
! 600: $ GO TO 250
! 601: END IF
! 602: *
! 603: * Note, eta should be positive if w is negative, and
! 604: * eta should be negative otherwise. However,
! 605: * if for some reason caused by roundoff, eta*w > 0,
! 606: * we simply use one Newton step instead. This way
! 607: * will guarantee eta*w < 0.
! 608: *
! 609: IF( W*ETA.GE.ZERO )
! 610: $ ETA = -W / DW
! 611: TEMP = TAU + ETA
! 612: IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
! 613: IF( W.LT.ZERO ) THEN
! 614: ETA = ( DLTUB-TAU ) / TWO
! 615: ELSE
! 616: ETA = ( DLTLB-TAU ) / TWO
! 617: END IF
! 618: END IF
! 619: *
! 620: PREW = W
! 621: *
! 622: DO 180 J = 1, N
! 623: DELTA( J ) = DELTA( J ) - ETA
! 624: 180 CONTINUE
! 625: *
! 626: * Evaluate PSI and the derivative DPSI
! 627: *
! 628: DPSI = ZERO
! 629: PSI = ZERO
! 630: ERRETM = ZERO
! 631: DO 190 J = 1, IIM1
! 632: TEMP = Z( J ) / DELTA( J )
! 633: PSI = PSI + Z( J )*TEMP
! 634: DPSI = DPSI + TEMP*TEMP
! 635: ERRETM = ERRETM + PSI
! 636: 190 CONTINUE
! 637: ERRETM = ABS( ERRETM )
! 638: *
! 639: * Evaluate PHI and the derivative DPHI
! 640: *
! 641: DPHI = ZERO
! 642: PHI = ZERO
! 643: DO 200 J = N, IIP1, -1
! 644: TEMP = Z( J ) / DELTA( J )
! 645: PHI = PHI + Z( J )*TEMP
! 646: DPHI = DPHI + TEMP*TEMP
! 647: ERRETM = ERRETM + PHI
! 648: 200 CONTINUE
! 649: *
! 650: TEMP = Z( II ) / DELTA( II )
! 651: DW = DPSI + DPHI + TEMP*TEMP
! 652: TEMP = Z( II )*TEMP
! 653: W = RHOINV + PHI + PSI + TEMP
! 654: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
! 655: $ THREE*ABS( TEMP ) + ABS( TAU+ETA )*DW
! 656: *
! 657: SWTCH = .FALSE.
! 658: IF( ORGATI ) THEN
! 659: IF( -W.GT.ABS( PREW ) / TEN )
! 660: $ SWTCH = .TRUE.
! 661: ELSE
! 662: IF( W.GT.ABS( PREW ) / TEN )
! 663: $ SWTCH = .TRUE.
! 664: END IF
! 665: *
! 666: TAU = TAU + ETA
! 667: *
! 668: * Main loop to update the values of the array DELTA
! 669: *
! 670: ITER = NITER + 1
! 671: *
! 672: DO 240 NITER = ITER, MAXIT
! 673: *
! 674: * Test for convergence
! 675: *
! 676: IF( ABS( W ).LE.EPS*ERRETM ) THEN
! 677: IF( ORGATI ) THEN
! 678: DLAM = D( I ) + TAU
! 679: ELSE
! 680: DLAM = D( IP1 ) + TAU
! 681: END IF
! 682: GO TO 250
! 683: END IF
! 684: *
! 685: IF( W.LE.ZERO ) THEN
! 686: DLTLB = MAX( DLTLB, TAU )
! 687: ELSE
! 688: DLTUB = MIN( DLTUB, TAU )
! 689: END IF
! 690: *
! 691: * Calculate the new step
! 692: *
! 693: IF( .NOT.SWTCH3 ) THEN
! 694: IF( .NOT.SWTCH ) THEN
! 695: IF( ORGATI ) THEN
! 696: C = W - DELTA( IP1 )*DW -
! 697: $ ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2
! 698: ELSE
! 699: C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
! 700: $ ( Z( IP1 ) / DELTA( IP1 ) )**2
! 701: END IF
! 702: ELSE
! 703: TEMP = Z( II ) / DELTA( II )
! 704: IF( ORGATI ) THEN
! 705: DPSI = DPSI + TEMP*TEMP
! 706: ELSE
! 707: DPHI = DPHI + TEMP*TEMP
! 708: END IF
! 709: C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI
! 710: END IF
! 711: A = ( DELTA( I )+DELTA( IP1 ) )*W -
! 712: $ DELTA( I )*DELTA( IP1 )*DW
! 713: B = DELTA( I )*DELTA( IP1 )*W
! 714: IF( C.EQ.ZERO ) THEN
! 715: IF( A.EQ.ZERO ) THEN
! 716: IF( .NOT.SWTCH ) THEN
! 717: IF( ORGATI ) THEN
! 718: A = Z( I )*Z( I ) + DELTA( IP1 )*
! 719: $ DELTA( IP1 )*( DPSI+DPHI )
! 720: ELSE
! 721: A = Z( IP1 )*Z( IP1 ) +
! 722: $ DELTA( I )*DELTA( I )*( DPSI+DPHI )
! 723: END IF
! 724: ELSE
! 725: A = DELTA( I )*DELTA( I )*DPSI +
! 726: $ DELTA( IP1 )*DELTA( IP1 )*DPHI
! 727: END IF
! 728: END IF
! 729: ETA = B / A
! 730: ELSE IF( A.LE.ZERO ) THEN
! 731: ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
! 732: ELSE
! 733: ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
! 734: END IF
! 735: ELSE
! 736: *
! 737: * Interpolation using THREE most relevant poles
! 738: *
! 739: TEMP = RHOINV + PSI + PHI
! 740: IF( SWTCH ) THEN
! 741: C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI
! 742: ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI
! 743: ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI
! 744: ELSE
! 745: IF( ORGATI ) THEN
! 746: TEMP1 = Z( IIM1 ) / DELTA( IIM1 )
! 747: TEMP1 = TEMP1*TEMP1
! 748: C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) -
! 749: $ ( D( IIM1 )-D( IIP1 ) )*TEMP1
! 750: ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
! 751: ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*
! 752: $ ( ( DPSI-TEMP1 )+DPHI )
! 753: ELSE
! 754: TEMP1 = Z( IIP1 ) / DELTA( IIP1 )
! 755: TEMP1 = TEMP1*TEMP1
! 756: C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) -
! 757: $ ( D( IIP1 )-D( IIM1 ) )*TEMP1
! 758: ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*
! 759: $ ( DPSI+( DPHI-TEMP1 ) )
! 760: ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
! 761: END IF
! 762: END IF
! 763: CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
! 764: $ INFO )
! 765: IF( INFO.NE.0 )
! 766: $ GO TO 250
! 767: END IF
! 768: *
! 769: * Note, eta should be positive if w is negative, and
! 770: * eta should be negative otherwise. However,
! 771: * if for some reason caused by roundoff, eta*w > 0,
! 772: * we simply use one Newton step instead. This way
! 773: * will guarantee eta*w < 0.
! 774: *
! 775: IF( W*ETA.GE.ZERO )
! 776: $ ETA = -W / DW
! 777: TEMP = TAU + ETA
! 778: IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
! 779: IF( W.LT.ZERO ) THEN
! 780: ETA = ( DLTUB-TAU ) / TWO
! 781: ELSE
! 782: ETA = ( DLTLB-TAU ) / TWO
! 783: END IF
! 784: END IF
! 785: *
! 786: DO 210 J = 1, N
! 787: DELTA( J ) = DELTA( J ) - ETA
! 788: 210 CONTINUE
! 789: *
! 790: TAU = TAU + ETA
! 791: PREW = W
! 792: *
! 793: * Evaluate PSI and the derivative DPSI
! 794: *
! 795: DPSI = ZERO
! 796: PSI = ZERO
! 797: ERRETM = ZERO
! 798: DO 220 J = 1, IIM1
! 799: TEMP = Z( J ) / DELTA( J )
! 800: PSI = PSI + Z( J )*TEMP
! 801: DPSI = DPSI + TEMP*TEMP
! 802: ERRETM = ERRETM + PSI
! 803: 220 CONTINUE
! 804: ERRETM = ABS( ERRETM )
! 805: *
! 806: * Evaluate PHI and the derivative DPHI
! 807: *
! 808: DPHI = ZERO
! 809: PHI = ZERO
! 810: DO 230 J = N, IIP1, -1
! 811: TEMP = Z( J ) / DELTA( J )
! 812: PHI = PHI + Z( J )*TEMP
! 813: DPHI = DPHI + TEMP*TEMP
! 814: ERRETM = ERRETM + PHI
! 815: 230 CONTINUE
! 816: *
! 817: TEMP = Z( II ) / DELTA( II )
! 818: DW = DPSI + DPHI + TEMP*TEMP
! 819: TEMP = Z( II )*TEMP
! 820: W = RHOINV + PHI + PSI + TEMP
! 821: ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
! 822: $ THREE*ABS( TEMP ) + ABS( TAU )*DW
! 823: IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
! 824: $ SWTCH = .NOT.SWTCH
! 825: *
! 826: 240 CONTINUE
! 827: *
! 828: * Return with INFO = 1, NITER = MAXIT and not converged
! 829: *
! 830: INFO = 1
! 831: IF( ORGATI ) THEN
! 832: DLAM = D( I ) + TAU
! 833: ELSE
! 834: DLAM = D( IP1 ) + TAU
! 835: END IF
! 836: *
! 837: END IF
! 838: *
! 839: 250 CONTINUE
! 840: *
! 841: RETURN
! 842: *
! 843: * End of DLAED4
! 844: *
! 845: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack