 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00005 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00006 * November 2006 00007 * 00008 * .. Scalar Arguments .. 00009 INTEGER I, INFO, N 00010 DOUBLE PRECISION DLAM, RHO 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION D( * ), DELTA( * ), Z( * ) 00014 * .. 00015 * 00016 * Purpose 00017 * ======= 00018 * 00019 * This subroutine computes the I-th updated eigenvalue of a symmetric 00020 * rank-one modification to a diagonal matrix whose elements are 00021 * given in the array d, and that 00022 * 00023 * D(i) < D(j) for i < j 00024 * 00025 * and that RHO > 0. This is arranged by the calling routine, and is 00026 * no loss in generality. The rank-one modified system is thus 00027 * 00028 * diag( D ) + RHO * Z * Z_transpose. 00029 * 00030 * where we assume the Euclidean norm of Z is 1. 00031 * 00032 * The method consists of approximating the rational functions in the 00033 * secular equation by simpler interpolating rational functions. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * N (input) INTEGER 00039 * The length of all arrays. 00040 * 00041 * I (input) INTEGER 00042 * The index of the eigenvalue to be computed. 1 <= I <= N. 00043 * 00044 * D (input) DOUBLE PRECISION array, dimension (N) 00045 * The original eigenvalues. It is assumed that they are in 00046 * order, D(I) < D(J) for I < J. 00047 * 00048 * Z (input) DOUBLE PRECISION array, dimension (N) 00049 * The components of the updating vector. 00050 * 00051 * DELTA (output) DOUBLE PRECISION array, dimension (N) 00052 * If N .GT. 2, DELTA contains (D(j) - lambda_I) in its j-th 00053 * component. If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5 00054 * for detail. The vector DELTA contains the information necessary 00055 * to construct the eigenvectors by DLAED3 and DLAED9. 00056 * 00057 * RHO (input) DOUBLE PRECISION 00058 * The scalar in the symmetric updating formula. 00059 * 00060 * DLAM (output) DOUBLE PRECISION 00061 * The computed lambda_I, the I-th updated eigenvalue. 00062 * 00063 * INFO (output) INTEGER 00064 * = 0: successful exit 00065 * > 0: if INFO = 1, the updating process failed. 00066 * 00067 * Internal Parameters 00068 * =================== 00069 * 00070 * Logical variable ORGATI (origin-at-i?) is used for distinguishing 00071 * whether D(i) or D(i+1) is treated as the origin. 00072 * 00073 * ORGATI = .true. origin at i 00074 * ORGATI = .false. origin at i+1 00075 * 00076 * Logical variable SWTCH3 (switch-for-3-poles?) is for noting 00077 * if we are working with THREE poles! 00078 * 00079 * MAXIT is the maximum number of iterations allowed for each 00080 * eigenvalue. 00081 * 00082 * Further Details 00083 * =============== 00084 * 00085 * Based on contributions by 00086 * Ren-Cang Li, Computer Science Division, University of California 00087 * at Berkeley, USA 00088 * 00089 * ===================================================================== 00090 * 00091 * .. Parameters .. 00092 INTEGER MAXIT 00093 PARAMETER ( MAXIT = 30 ) 00094 DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN 00095 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, 00096 $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0, 00097 $ TEN = 10.0D0 ) 00098 * .. 00099 * .. Local Scalars .. 00100 LOGICAL ORGATI, SWTCH, SWTCH3 00101 INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER 00102 DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW, 00103 $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI, 00104 $ RHOINV, TAU, TEMP, TEMP1, W 00105 * .. 00106 * .. Local Arrays .. 00107 DOUBLE PRECISION ZZ( 3 ) 00108 * .. 00109 * .. External Functions .. 00110 DOUBLE PRECISION DLAMCH 00111 EXTERNAL DLAMCH 00112 * .. 00113 * .. External Subroutines .. 00114 EXTERNAL DLAED5, DLAED6 00115 * .. 00116 * .. Intrinsic Functions .. 00117 INTRINSIC ABS, MAX, MIN, SQRT 00118 * .. 00119 * .. Executable Statements .. 00120 * 00121 * Since this routine is called in an inner loop, we do no argument 00122 * checking. 00123 * 00124 * Quick return for N=1 and 2. 00125 * 00126 INFO = 0 00127 IF( N.EQ.1 ) THEN 00128 * 00129 * Presumably, I=1 upon entry 00130 * 00131 DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 ) 00132 DELTA( 1 ) = ONE 00133 RETURN 00134 END IF 00135 IF( N.EQ.2 ) THEN 00136 CALL DLAED5( I, D, Z, DELTA, RHO, DLAM ) 00137 RETURN 00138 END IF 00139 * 00140 * Compute machine epsilon 00141 * 00142 EPS = DLAMCH( 'Epsilon' ) 00143 RHOINV = ONE / RHO 00144 * 00145 * The case I = N 00146 * 00147 IF( I.EQ.N ) THEN 00148 * 00149 * Initialize some basic variables 00150 * 00151 II = N - 1 00152 NITER = 1 00153 * 00154 * Calculate initial guess 00155 * 00156 MIDPT = RHO / TWO 00157 * 00158 * If ||Z||_2 is not one, then TEMP should be set to 00159 * RHO * ||Z||_2^2 / TWO 00160 * 00161 DO 10 J = 1, N 00162 DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 00163 10 CONTINUE 00164 * 00165 PSI = ZERO 00166 DO 20 J = 1, N - 2 00167 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 00168 20 CONTINUE 00169 * 00170 C = RHOINV + PSI 00171 W = C + Z( II )*Z( II ) / DELTA( II ) + 00172 $ Z( N )*Z( N ) / DELTA( N ) 00173 * 00174 IF( W.LE.ZERO ) THEN 00175 TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) + 00176 $ Z( N )*Z( N ) / RHO 00177 IF( C.LE.TEMP ) THEN 00178 TAU = RHO 00179 ELSE 00180 DEL = D( N ) - D( N-1 ) 00181 A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) 00182 B = Z( N )*Z( N )*DEL 00183 IF( A.LT.ZERO ) THEN 00184 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) 00185 ELSE 00186 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) 00187 END IF 00188 END IF 00189 * 00190 * It can be proved that 00191 * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO 00192 * 00193 DLTLB = MIDPT 00194 DLTUB = RHO 00195 ELSE 00196 DEL = D( N ) - D( N-1 ) 00197 A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) 00198 B = Z( N )*Z( N )*DEL 00199 IF( A.LT.ZERO ) THEN 00200 TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) 00201 ELSE 00202 TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) 00203 END IF 00204 * 00205 * It can be proved that 00206 * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 00207 * 00208 DLTLB = ZERO 00209 DLTUB = MIDPT 00210 END IF 00211 * 00212 DO 30 J = 1, N 00213 DELTA( J ) = ( D( J )-D( I ) ) - TAU 00214 30 CONTINUE 00215 * 00216 * Evaluate PSI and the derivative DPSI 00217 * 00218 DPSI = ZERO 00219 PSI = ZERO 00220 ERRETM = ZERO 00221 DO 40 J = 1, II 00222 TEMP = Z( J ) / DELTA( J ) 00223 PSI = PSI + Z( J )*TEMP 00224 DPSI = DPSI + TEMP*TEMP 00225 ERRETM = ERRETM + PSI 00226 40 CONTINUE 00227 ERRETM = ABS( ERRETM ) 00228 * 00229 * Evaluate PHI and the derivative DPHI 00230 * 00231 TEMP = Z( N ) / DELTA( N ) 00232 PHI = Z( N )*TEMP 00233 DPHI = TEMP*TEMP 00234 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + 00235 $ ABS( TAU )*( DPSI+DPHI ) 00236 * 00237 W = RHOINV + PHI + PSI 00238 * 00239 * Test for convergence 00240 * 00241 IF( ABS( W ).LE.EPS*ERRETM ) THEN 00242 DLAM = D( I ) + TAU 00243 GO TO 250 00244 END IF 00245 * 00246 IF( W.LE.ZERO ) THEN 00247 DLTLB = MAX( DLTLB, TAU ) 00248 ELSE 00249 DLTUB = MIN( DLTUB, TAU ) 00250 END IF 00251 * 00252 * Calculate the new step 00253 * 00254 NITER = NITER + 1 00255 C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI 00256 A = ( DELTA( N-1 )+DELTA( N ) )*W - 00257 $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) 00258 B = DELTA( N-1 )*DELTA( N )*W 00259 IF( C.LT.ZERO ) 00260 $ C = ABS( C ) 00261 IF( C.EQ.ZERO ) THEN 00262 * ETA = B/A 00263 * ETA = RHO - TAU 00264 ETA = DLTUB - TAU 00265 ELSE IF( A.GE.ZERO ) THEN 00266 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00267 ELSE 00268 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) 00269 END IF 00270 * 00271 * Note, eta should be positive if w is negative, and 00272 * eta should be negative otherwise. However, 00273 * if for some reason caused by roundoff, eta*w > 0, 00274 * we simply use one Newton step instead. This way 00275 * will guarantee eta*w < 0. 00276 * 00277 IF( W*ETA.GT.ZERO ) 00278 $ ETA = -W / ( DPSI+DPHI ) 00279 TEMP = TAU + ETA 00280 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 00281 IF( W.LT.ZERO ) THEN 00282 ETA = ( DLTUB-TAU ) / TWO 00283 ELSE 00284 ETA = ( DLTLB-TAU ) / TWO 00285 END IF 00286 END IF 00287 DO 50 J = 1, N 00288 DELTA( J ) = DELTA( J ) - ETA 00289 50 CONTINUE 00290 * 00291 TAU = TAU + ETA 00292 * 00293 * Evaluate PSI and the derivative DPSI 00294 * 00295 DPSI = ZERO 00296 PSI = ZERO 00297 ERRETM = ZERO 00298 DO 60 J = 1, II 00299 TEMP = Z( J ) / DELTA( J ) 00300 PSI = PSI + Z( J )*TEMP 00301 DPSI = DPSI + TEMP*TEMP 00302 ERRETM = ERRETM + PSI 00303 60 CONTINUE 00304 ERRETM = ABS( ERRETM ) 00305 * 00306 * Evaluate PHI and the derivative DPHI 00307 * 00308 TEMP = Z( N ) / DELTA( N ) 00309 PHI = Z( N )*TEMP 00310 DPHI = TEMP*TEMP 00311 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + 00312 $ ABS( TAU )*( DPSI+DPHI ) 00313 * 00314 W = RHOINV + PHI + PSI 00315 * 00316 * Main loop to update the values of the array DELTA 00317 * 00318 ITER = NITER + 1 00319 * 00320 DO 90 NITER = ITER, MAXIT 00321 * 00322 * Test for convergence 00323 * 00324 IF( ABS( W ).LE.EPS*ERRETM ) THEN 00325 DLAM = D( I ) + TAU 00326 GO TO 250 00327 END IF 00328 * 00329 IF( W.LE.ZERO ) THEN 00330 DLTLB = MAX( DLTLB, TAU ) 00331 ELSE 00332 DLTUB = MIN( DLTUB, TAU ) 00333 END IF 00334 * 00335 * Calculate the new step 00336 * 00337 C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI 00338 A = ( DELTA( N-1 )+DELTA( N ) )*W - 00339 $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI ) 00340 B = DELTA( N-1 )*DELTA( N )*W 00341 IF( A.GE.ZERO ) THEN 00342 ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00343 ELSE 00344 ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) 00345 END IF 00346 * 00347 * Note, eta should be positive if w is negative, and 00348 * eta should be negative otherwise. However, 00349 * if for some reason caused by roundoff, eta*w > 0, 00350 * we simply use one Newton step instead. This way 00351 * will guarantee eta*w < 0. 00352 * 00353 IF( W*ETA.GT.ZERO ) 00354 $ ETA = -W / ( DPSI+DPHI ) 00355 TEMP = TAU + ETA 00356 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 00357 IF( W.LT.ZERO ) THEN 00358 ETA = ( DLTUB-TAU ) / TWO 00359 ELSE 00360 ETA = ( DLTLB-TAU ) / TWO 00361 END IF 00362 END IF 00363 DO 70 J = 1, N 00364 DELTA( J ) = DELTA( J ) - ETA 00365 70 CONTINUE 00366 * 00367 TAU = TAU + ETA 00368 * 00369 * Evaluate PSI and the derivative DPSI 00370 * 00371 DPSI = ZERO 00372 PSI = ZERO 00373 ERRETM = ZERO 00374 DO 80 J = 1, II 00375 TEMP = Z( J ) / DELTA( J ) 00376 PSI = PSI + Z( J )*TEMP 00377 DPSI = DPSI + TEMP*TEMP 00378 ERRETM = ERRETM + PSI 00379 80 CONTINUE 00380 ERRETM = ABS( ERRETM ) 00381 * 00382 * Evaluate PHI and the derivative DPHI 00383 * 00384 TEMP = Z( N ) / DELTA( N ) 00385 PHI = Z( N )*TEMP 00386 DPHI = TEMP*TEMP 00387 ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + 00388 $ ABS( TAU )*( DPSI+DPHI ) 00389 * 00390 W = RHOINV + PHI + PSI 00391 90 CONTINUE 00392 * 00393 * Return with INFO = 1, NITER = MAXIT and not converged 00394 * 00395 INFO = 1 00396 DLAM = D( I ) + TAU 00397 GO TO 250 00398 * 00399 * End for the case I = N 00400 * 00401 ELSE 00402 * 00403 * The case for I < N 00404 * 00405 NITER = 1 00406 IP1 = I + 1 00407 * 00408 * Calculate initial guess 00409 * 00410 DEL = D( IP1 ) - D( I ) 00411 MIDPT = DEL / TWO 00412 DO 100 J = 1, N 00413 DELTA( J ) = ( D( J )-D( I ) ) - MIDPT 00414 100 CONTINUE 00415 * 00416 PSI = ZERO 00417 DO 110 J = 1, I - 1 00418 PSI = PSI + Z( J )*Z( J ) / DELTA( J ) 00419 110 CONTINUE 00420 * 00421 PHI = ZERO 00422 DO 120 J = N, I + 2, -1 00423 PHI = PHI + Z( J )*Z( J ) / DELTA( J ) 00424 120 CONTINUE 00425 C = RHOINV + PSI + PHI 00426 W = C + Z( I )*Z( I ) / DELTA( I ) + 00427 $ Z( IP1 )*Z( IP1 ) / DELTA( IP1 ) 00428 * 00429 IF( W.GT.ZERO ) THEN 00430 * 00431 * d(i)< the ith eigenvalue < (d(i)+d(i+1))/2 00432 * 00433 * We choose d(i) as origin. 00434 * 00435 ORGATI = .TRUE. 00436 A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 ) 00437 B = Z( I )*Z( I )*DEL 00438 IF( A.GT.ZERO ) THEN 00439 TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 00440 ELSE 00441 TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00442 END IF 00443 DLTLB = ZERO 00444 DLTUB = MIDPT 00445 ELSE 00446 * 00447 * (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1) 00448 * 00449 * We choose d(i+1) as origin. 00450 * 00451 ORGATI = .FALSE. 00452 A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 ) 00453 B = Z( IP1 )*Z( IP1 )*DEL 00454 IF( A.LT.ZERO ) THEN 00455 TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) ) 00456 ELSE 00457 TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C ) 00458 END IF 00459 DLTLB = -MIDPT 00460 DLTUB = ZERO 00461 END IF 00462 * 00463 IF( ORGATI ) THEN 00464 DO 130 J = 1, N 00465 DELTA( J ) = ( D( J )-D( I ) ) - TAU 00466 130 CONTINUE 00467 ELSE 00468 DO 140 J = 1, N 00469 DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU 00470 140 CONTINUE 00471 END IF 00472 IF( ORGATI ) THEN 00473 II = I 00474 ELSE 00475 II = I + 1 00476 END IF 00477 IIM1 = II - 1 00478 IIP1 = II + 1 00479 * 00480 * Evaluate PSI and the derivative DPSI 00481 * 00482 DPSI = ZERO 00483 PSI = ZERO 00484 ERRETM = ZERO 00485 DO 150 J = 1, IIM1 00486 TEMP = Z( J ) / DELTA( J ) 00487 PSI = PSI + Z( J )*TEMP 00488 DPSI = DPSI + TEMP*TEMP 00489 ERRETM = ERRETM + PSI 00490 150 CONTINUE 00491 ERRETM = ABS( ERRETM ) 00492 * 00493 * Evaluate PHI and the derivative DPHI 00494 * 00495 DPHI = ZERO 00496 PHI = ZERO 00497 DO 160 J = N, IIP1, -1 00498 TEMP = Z( J ) / DELTA( J ) 00499 PHI = PHI + Z( J )*TEMP 00500 DPHI = DPHI + TEMP*TEMP 00501 ERRETM = ERRETM + PHI 00502 160 CONTINUE 00503 * 00504 W = RHOINV + PHI + PSI 00505 * 00506 * W is the value of the secular function with 00507 * its ii-th element removed. 00508 * 00509 SWTCH3 = .FALSE. 00510 IF( ORGATI ) THEN 00511 IF( W.LT.ZERO ) 00512 $ SWTCH3 = .TRUE. 00513 ELSE 00514 IF( W.GT.ZERO ) 00515 $ SWTCH3 = .TRUE. 00516 END IF 00517 IF( II.EQ.1 .OR. II.EQ.N ) 00518 $ SWTCH3 = .FALSE. 00519 * 00520 TEMP = Z( II ) / DELTA( II ) 00521 DW = DPSI + DPHI + TEMP*TEMP 00522 TEMP = Z( II )*TEMP 00523 W = W + TEMP 00524 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + 00525 $ THREE*ABS( TEMP ) + ABS( TAU )*DW 00526 * 00527 * Test for convergence 00528 * 00529 IF( ABS( W ).LE.EPS*ERRETM ) THEN 00530 IF( ORGATI ) THEN 00531 DLAM = D( I ) + TAU 00532 ELSE 00533 DLAM = D( IP1 ) + TAU 00534 END IF 00535 GO TO 250 00536 END IF 00537 * 00538 IF( W.LE.ZERO ) THEN 00539 DLTLB = MAX( DLTLB, TAU ) 00540 ELSE 00541 DLTUB = MIN( DLTUB, TAU ) 00542 END IF 00543 * 00544 * Calculate the new step 00545 * 00546 NITER = NITER + 1 00547 IF( .NOT.SWTCH3 ) THEN 00548 IF( ORGATI ) THEN 00549 C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )* 00550 $ ( Z( I ) / DELTA( I ) )**2 00551 ELSE 00552 C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* 00553 $ ( Z( IP1 ) / DELTA( IP1 ) )**2 00554 END IF 00555 A = ( DELTA( I )+DELTA( IP1 ) )*W - 00556 $ DELTA( I )*DELTA( IP1 )*DW 00557 B = DELTA( I )*DELTA( IP1 )*W 00558 IF( C.EQ.ZERO ) THEN 00559 IF( A.EQ.ZERO ) THEN 00560 IF( ORGATI ) THEN 00561 A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )* 00562 $ ( DPSI+DPHI ) 00563 ELSE 00564 A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )* 00565 $ ( DPSI+DPHI ) 00566 END IF 00567 END IF 00568 ETA = B / A 00569 ELSE IF( A.LE.ZERO ) THEN 00570 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00571 ELSE 00572 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 00573 END IF 00574 ELSE 00575 * 00576 * Interpolation using THREE most relevant poles 00577 * 00578 TEMP = RHOINV + PSI + PHI 00579 IF( ORGATI ) THEN 00580 TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) 00581 TEMP1 = TEMP1*TEMP1 00582 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - 00583 $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 00584 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) 00585 ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* 00586 $ ( ( DPSI-TEMP1 )+DPHI ) 00587 ELSE 00588 TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) 00589 TEMP1 = TEMP1*TEMP1 00590 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - 00591 $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 00592 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* 00593 $ ( DPSI+( DPHI-TEMP1 ) ) 00594 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) 00595 END IF 00596 ZZ( 2 ) = Z( II )*Z( II ) 00597 CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, 00598 $ INFO ) 00599 IF( INFO.NE.0 ) 00600 $ GO TO 250 00601 END IF 00602 * 00603 * Note, eta should be positive if w is negative, and 00604 * eta should be negative otherwise. However, 00605 * if for some reason caused by roundoff, eta*w > 0, 00606 * we simply use one Newton step instead. This way 00607 * will guarantee eta*w < 0. 00608 * 00609 IF( W*ETA.GE.ZERO ) 00610 $ ETA = -W / DW 00611 TEMP = TAU + ETA 00612 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 00613 IF( W.LT.ZERO ) THEN 00614 ETA = ( DLTUB-TAU ) / TWO 00615 ELSE 00616 ETA = ( DLTLB-TAU ) / TWO 00617 END IF 00618 END IF 00619 * 00620 PREW = W 00621 * 00622 DO 180 J = 1, N 00623 DELTA( J ) = DELTA( J ) - ETA 00624 180 CONTINUE 00625 * 00626 * Evaluate PSI and the derivative DPSI 00627 * 00628 DPSI = ZERO 00629 PSI = ZERO 00630 ERRETM = ZERO 00631 DO 190 J = 1, IIM1 00632 TEMP = Z( J ) / DELTA( J ) 00633 PSI = PSI + Z( J )*TEMP 00634 DPSI = DPSI + TEMP*TEMP 00635 ERRETM = ERRETM + PSI 00636 190 CONTINUE 00637 ERRETM = ABS( ERRETM ) 00638 * 00639 * Evaluate PHI and the derivative DPHI 00640 * 00641 DPHI = ZERO 00642 PHI = ZERO 00643 DO 200 J = N, IIP1, -1 00644 TEMP = Z( J ) / DELTA( J ) 00645 PHI = PHI + Z( J )*TEMP 00646 DPHI = DPHI + TEMP*TEMP 00647 ERRETM = ERRETM + PHI 00648 200 CONTINUE 00649 * 00650 TEMP = Z( II ) / DELTA( II ) 00651 DW = DPSI + DPHI + TEMP*TEMP 00652 TEMP = Z( II )*TEMP 00653 W = RHOINV + PHI + PSI + TEMP 00654 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + 00655 $ THREE*ABS( TEMP ) + ABS( TAU+ETA )*DW 00656 * 00657 SWTCH = .FALSE. 00658 IF( ORGATI ) THEN 00659 IF( -W.GT.ABS( PREW ) / TEN ) 00660 $ SWTCH = .TRUE. 00661 ELSE 00662 IF( W.GT.ABS( PREW ) / TEN ) 00663 $ SWTCH = .TRUE. 00664 END IF 00665 * 00666 TAU = TAU + ETA 00667 * 00668 * Main loop to update the values of the array DELTA 00669 * 00670 ITER = NITER + 1 00671 * 00672 DO 240 NITER = ITER, MAXIT 00673 * 00674 * Test for convergence 00675 * 00676 IF( ABS( W ).LE.EPS*ERRETM ) THEN 00677 IF( ORGATI ) THEN 00678 DLAM = D( I ) + TAU 00679 ELSE 00680 DLAM = D( IP1 ) + TAU 00681 END IF 00682 GO TO 250 00683 END IF 00684 * 00685 IF( W.LE.ZERO ) THEN 00686 DLTLB = MAX( DLTLB, TAU ) 00687 ELSE 00688 DLTUB = MIN( DLTUB, TAU ) 00689 END IF 00690 * 00691 * Calculate the new step 00692 * 00693 IF( .NOT.SWTCH3 ) THEN 00694 IF( .NOT.SWTCH ) THEN 00695 IF( ORGATI ) THEN 00696 C = W - DELTA( IP1 )*DW - 00697 $ ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2 00698 ELSE 00699 C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )* 00700 $ ( Z( IP1 ) / DELTA( IP1 ) )**2 00701 END IF 00702 ELSE 00703 TEMP = Z( II ) / DELTA( II ) 00704 IF( ORGATI ) THEN 00705 DPSI = DPSI + TEMP*TEMP 00706 ELSE 00707 DPHI = DPHI + TEMP*TEMP 00708 END IF 00709 C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI 00710 END IF 00711 A = ( DELTA( I )+DELTA( IP1 ) )*W - 00712 $ DELTA( I )*DELTA( IP1 )*DW 00713 B = DELTA( I )*DELTA( IP1 )*W 00714 IF( C.EQ.ZERO ) THEN 00715 IF( A.EQ.ZERO ) THEN 00716 IF( .NOT.SWTCH ) THEN 00717 IF( ORGATI ) THEN 00718 A = Z( I )*Z( I ) + DELTA( IP1 )* 00719 $ DELTA( IP1 )*( DPSI+DPHI ) 00720 ELSE 00721 A = Z( IP1 )*Z( IP1 ) + 00722 $ DELTA( I )*DELTA( I )*( DPSI+DPHI ) 00723 END IF 00724 ELSE 00725 A = DELTA( I )*DELTA( I )*DPSI + 00726 $ DELTA( IP1 )*DELTA( IP1 )*DPHI 00727 END IF 00728 END IF 00729 ETA = B / A 00730 ELSE IF( A.LE.ZERO ) THEN 00731 ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) 00732 ELSE 00733 ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) 00734 END IF 00735 ELSE 00736 * 00737 * Interpolation using THREE most relevant poles 00738 * 00739 TEMP = RHOINV + PSI + PHI 00740 IF( SWTCH ) THEN 00741 C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI 00742 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI 00743 ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI 00744 ELSE 00745 IF( ORGATI ) THEN 00746 TEMP1 = Z( IIM1 ) / DELTA( IIM1 ) 00747 TEMP1 = TEMP1*TEMP1 00748 C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) - 00749 $ ( D( IIM1 )-D( IIP1 ) )*TEMP1 00750 ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) 00751 ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )* 00752 $ ( ( DPSI-TEMP1 )+DPHI ) 00753 ELSE 00754 TEMP1 = Z( IIP1 ) / DELTA( IIP1 ) 00755 TEMP1 = TEMP1*TEMP1 00756 C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) - 00757 $ ( D( IIP1 )-D( IIM1 ) )*TEMP1 00758 ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )* 00759 $ ( DPSI+( DPHI-TEMP1 ) ) 00760 ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) 00761 END IF 00762 END IF 00763 CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA, 00764 $ INFO ) 00765 IF( INFO.NE.0 ) 00766 $ GO TO 250 00767 END IF 00768 * 00769 * Note, eta should be positive if w is negative, and 00770 * eta should be negative otherwise. However, 00771 * if for some reason caused by roundoff, eta*w > 0, 00772 * we simply use one Newton step instead. This way 00773 * will guarantee eta*w < 0. 00774 * 00775 IF( W*ETA.GE.ZERO ) 00776 $ ETA = -W / DW 00777 TEMP = TAU + ETA 00778 IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN 00779 IF( W.LT.ZERO ) THEN 00780 ETA = ( DLTUB-TAU ) / TWO 00781 ELSE 00782 ETA = ( DLTLB-TAU ) / TWO 00783 END IF 00784 END IF 00785 * 00786 DO 210 J = 1, N 00787 DELTA( J ) = DELTA( J ) - ETA 00788 210 CONTINUE 00789 * 00790 TAU = TAU + ETA 00791 PREW = W 00792 * 00793 * Evaluate PSI and the derivative DPSI 00794 * 00795 DPSI = ZERO 00796 PSI = ZERO 00797 ERRETM = ZERO 00798 DO 220 J = 1, IIM1 00799 TEMP = Z( J ) / DELTA( J ) 00800 PSI = PSI + Z( J )*TEMP 00801 DPSI = DPSI + TEMP*TEMP 00802 ERRETM = ERRETM + PSI 00803 220 CONTINUE 00804 ERRETM = ABS( ERRETM ) 00805 * 00806 * Evaluate PHI and the derivative DPHI 00807 * 00808 DPHI = ZERO 00809 PHI = ZERO 00810 DO 230 J = N, IIP1, -1 00811 TEMP = Z( J ) / DELTA( J ) 00812 PHI = PHI + Z( J )*TEMP 00813 DPHI = DPHI + TEMP*TEMP 00814 ERRETM = ERRETM + PHI 00815 230 CONTINUE 00816 * 00817 TEMP = Z( II ) / DELTA( II ) 00818 DW = DPSI + DPHI + TEMP*TEMP 00819 TEMP = Z( II )*TEMP 00820 W = RHOINV + PHI + PSI + TEMP 00821 ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + 00822 $ THREE*ABS( TEMP ) + ABS( TAU )*DW 00823 IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN ) 00824 $ SWTCH = .NOT.SWTCH 00825 * 00826 240 CONTINUE 00827 * 00828 * Return with INFO = 1, NITER = MAXIT and not converged 00829 * 00830 INFO = 1 00831 IF( ORGATI ) THEN 00832 DLAM = D( I ) + TAU 00833 ELSE 00834 DLAM = D( IP1 ) + TAU 00835 END IF 00836 * 00837 END IF 00838 * 00839 250 CONTINUE 00840 * 00841 RETURN 00842 * 00843 * End of DLAED4 00844 * 00845 END
1.7.3 
