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