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