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