Annotation of rpl/lapack/lapack/dlasd4.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
! 2: *
! 3: * -- LAPACK auxiliary 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 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 = 20 )
! 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack