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