root / src / lapack / double / dlasd4.f @ 9
Historique | Voir | Annoter | Télécharger (27,1 ko)
| 1 |
SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO ) |
|---|---|
| 2 |
* |
| 3 |
* -- LAPACK auxiliary routine (version 3.2) -- |
| 4 |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| 5 |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| 6 |
* November 2006 |
| 7 |
* |
| 8 |
* .. Scalar Arguments .. |
| 9 |
INTEGER I, INFO, N |
| 10 |
DOUBLE PRECISION RHO, SIGMA |
| 11 |
* .. |
| 12 |
* .. Array Arguments .. |
| 13 |
DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * ) |
| 14 |
* .. |
| 15 |
* |
| 16 |
* Purpose |
| 17 |
* ======= |
| 18 |
* |
| 19 |
* This subroutine computes the square root of the I-th updated |
| 20 |
* eigenvalue of a positive symmetric rank-one modification to |
| 21 |
* a positive diagonal matrix whose entries are given as the squares |
| 22 |
* of the corresponding entries in the array d, and that |
| 23 |
* |
| 24 |
* 0 <= D(i) < D(j) for i < j |
| 25 |
* |
| 26 |
* and that RHO > 0. This is arranged by the calling routine, and is |
| 27 |
* no loss in generality. The rank-one modified system is thus |
| 28 |
* |
| 29 |
* diag( D ) * diag( D ) + RHO * Z * Z_transpose. |
| 30 |
* |
| 31 |
* where we assume the Euclidean norm of Z is 1. |
| 32 |
* |
| 33 |
* The method consists of approximating the rational functions in the |
| 34 |
* secular equation by simpler interpolating rational functions. |
| 35 |
* |
| 36 |
* Arguments |
| 37 |
* ========= |
| 38 |
* |
| 39 |
* N (input) INTEGER |
| 40 |
* The length of all arrays. |
| 41 |
* |
| 42 |
* I (input) INTEGER |
| 43 |
* The index of the eigenvalue to be computed. 1 <= I <= N. |
| 44 |
* |
| 45 |
* D (input) DOUBLE PRECISION array, dimension ( N ) |
| 46 |
* The original eigenvalues. It is assumed that they are in |
| 47 |
* order, 0 <= D(I) < D(J) for I < J. |
| 48 |
* |
| 49 |
* Z (input) DOUBLE PRECISION array, dimension ( N ) |
| 50 |
* The components of the updating vector. |
| 51 |
* |
| 52 |
* DELTA (output) DOUBLE PRECISION array, dimension ( N ) |
| 53 |
* If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th |
| 54 |
* component. If N = 1, then DELTA(1) = 1. The vector DELTA |
| 55 |
* contains the information necessary to construct the |
| 56 |
* (singular) eigenvectors. |
| 57 |
* |
| 58 |
* RHO (input) DOUBLE PRECISION |
| 59 |
* The scalar in the symmetric updating formula. |
| 60 |
* |
| 61 |
* SIGMA (output) DOUBLE PRECISION |
| 62 |
* The computed sigma_I, the I-th updated eigenvalue. |
| 63 |
* |
| 64 |
* WORK (workspace) DOUBLE PRECISION array, dimension ( N ) |
| 65 |
* If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th |
| 66 |
* component. If N = 1, then WORK( 1 ) = 1. |
| 67 |
* |
| 68 |
* INFO (output) INTEGER |
| 69 |
* = 0: successful exit |
| 70 |
* > 0: if INFO = 1, the updating process failed. |
| 71 |
* |
| 72 |
* Internal Parameters |
| 73 |
* =================== |
| 74 |
* |
| 75 |
* Logical variable ORGATI (origin-at-i?) is used for distinguishing |
| 76 |
* whether D(i) or D(i+1) is treated as the origin. |
| 77 |
* |
| 78 |
* ORGATI = .true. origin at i |
| 79 |
* ORGATI = .false. origin at i+1 |
| 80 |
* |
| 81 |
* Logical variable SWTCH3 (switch-for-3-poles?) is for noting |
| 82 |
* if we are working with THREE poles! |
| 83 |
* |
| 84 |
* MAXIT is the maximum number of iterations allowed for each |
| 85 |
* eigenvalue. |
| 86 |
* |
| 87 |
* Further Details |
| 88 |
* =============== |
| 89 |
* |
| 90 |
* Based on contributions by |
| 91 |
* Ren-Cang Li, Computer Science Division, University of California |
| 92 |
* at Berkeley, USA |
| 93 |
* |
| 94 |
* ===================================================================== |
| 95 |
* |
| 96 |
* .. Parameters .. |
| 97 |
INTEGER MAXIT |
| 98 |
PARAMETER ( MAXIT = 20 ) |
| 99 |
DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN |
| 100 |
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0, |
| 101 |
$ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0, |
| 102 |
$ TEN = 10.0D+0 ) |
| 103 |
* .. |
| 104 |
* .. Local Scalars .. |
| 105 |
LOGICAL ORGATI, SWTCH, SWTCH3 |
| 106 |
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER |
| 107 |
DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM, |
| 108 |
$ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS, |
| 109 |
$ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB, |
| 110 |
$ SG2UB, TAU, TEMP, TEMP1, TEMP2, W |
| 111 |
* .. |
| 112 |
* .. Local Arrays .. |
| 113 |
DOUBLE PRECISION DD( 3 ), ZZ( 3 ) |
| 114 |
* .. |
| 115 |
* .. External Subroutines .. |
| 116 |
EXTERNAL DLAED6, DLASD5 |
| 117 |
* .. |
| 118 |
* .. External Functions .. |
| 119 |
DOUBLE PRECISION DLAMCH |
| 120 |
EXTERNAL DLAMCH |
| 121 |
* .. |
| 122 |
* .. Intrinsic Functions .. |
| 123 |
INTRINSIC ABS, MAX, MIN, SQRT |
| 124 |
* .. |
| 125 |
* .. Executable Statements .. |
| 126 |
* |
| 127 |
* Since this routine is called in an inner loop, we do no argument |
| 128 |
* checking. |
| 129 |
* |
| 130 |
* Quick return for N=1 and 2. |
| 131 |
* |
| 132 |
INFO = 0 |
| 133 |
IF( N.EQ.1 ) THEN |
| 134 |
* |
| 135 |
* Presumably, I=1 upon entry |
| 136 |
* |
| 137 |
SIGMA = SQRT( D( 1 )*D( 1 )+RHO*Z( 1 )*Z( 1 ) ) |
| 138 |
DELTA( 1 ) = ONE |
| 139 |
WORK( 1 ) = ONE |
| 140 |
RETURN |
| 141 |
END IF |
| 142 |
IF( N.EQ.2 ) THEN |
| 143 |
CALL DLASD5( I, D, Z, DELTA, RHO, SIGMA, WORK ) |
| 144 |
RETURN |
| 145 |
END IF |
| 146 |
* |
| 147 |
* Compute machine epsilon |
| 148 |
* |
| 149 |
EPS = DLAMCH( 'Epsilon' ) |
| 150 |
RHOINV = ONE / RHO |
| 151 |
* |
| 152 |
* The case I = N |
| 153 |
* |
| 154 |
IF( I.EQ.N ) THEN |
| 155 |
* |
| 156 |
* Initialize some basic variables |
| 157 |
* |
| 158 |
II = N - 1 |
| 159 |
NITER = 1 |
| 160 |
* |
| 161 |
* Calculate initial guess |
| 162 |
* |
| 163 |
TEMP = RHO / TWO |
| 164 |
* |
| 165 |
* If ||Z||_2 is not one, then TEMP should be set to |
| 166 |
* RHO * ||Z||_2^2 / TWO |
| 167 |
* |
| 168 |
TEMP1 = TEMP / ( D( N )+SQRT( D( N )*D( N )+TEMP ) ) |
| 169 |
DO 10 J = 1, N |
| 170 |
WORK( J ) = D( J ) + D( N ) + TEMP1 |
| 171 |
DELTA( J ) = ( D( J )-D( N ) ) - TEMP1 |
| 172 |
10 CONTINUE |
| 173 |
* |
| 174 |
PSI = ZERO |
| 175 |
DO 20 J = 1, N - 2 |
| 176 |
PSI = PSI + Z( J )*Z( J ) / ( DELTA( J )*WORK( J ) ) |
| 177 |
20 CONTINUE |
| 178 |
* |
| 179 |
C = RHOINV + PSI |
| 180 |
W = C + Z( II )*Z( II ) / ( DELTA( II )*WORK( II ) ) + |
| 181 |
$ Z( N )*Z( N ) / ( DELTA( N )*WORK( N ) ) |
| 182 |
* |
| 183 |
IF( W.LE.ZERO ) THEN |
| 184 |
TEMP1 = SQRT( D( N )*D( N )+RHO ) |
| 185 |
TEMP = Z( N-1 )*Z( N-1 ) / ( ( D( N-1 )+TEMP1 )* |
| 186 |
$ ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) + |
| 187 |
$ Z( N )*Z( N ) / RHO |
| 188 |
* |
| 189 |
* The following TAU is to approximate |
| 190 |
* SIGMA_n^2 - D( N )*D( N ) |
| 191 |
* |
| 192 |
IF( C.LE.TEMP ) THEN |
| 193 |
TAU = RHO |
| 194 |
ELSE |
| 195 |
DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) ) |
| 196 |
A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) |
| 197 |
B = Z( N )*Z( N )*DELSQ |
| 198 |
IF( A.LT.ZERO ) THEN |
| 199 |
TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) |
| 200 |
ELSE |
| 201 |
TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) |
| 202 |
END IF |
| 203 |
END IF |
| 204 |
* |
| 205 |
* It can be proved that |
| 206 |
* D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO |
| 207 |
* |
| 208 |
ELSE |
| 209 |
DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) ) |
| 210 |
A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N ) |
| 211 |
B = Z( N )*Z( N )*DELSQ |
| 212 |
* |
| 213 |
* The following TAU is to approximate |
| 214 |
* SIGMA_n^2 - D( N )*D( N ) |
| 215 |
* |
| 216 |
IF( A.LT.ZERO ) THEN |
| 217 |
TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A ) |
| 218 |
ELSE |
| 219 |
TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C ) |
| 220 |
END IF |
| 221 |
* |
| 222 |
* It can be proved that |
| 223 |
* D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2 |
| 224 |
* |
| 225 |
END IF |
| 226 |
* |
| 227 |
* The following ETA is to approximate SIGMA_n - D( N ) |
| 228 |
* |
| 229 |
ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) ) |
| 230 |
* |
| 231 |
SIGMA = D( N ) + ETA |
| 232 |
DO 30 J = 1, N |
| 233 |
DELTA( J ) = ( D( J )-D( I ) ) - ETA |
| 234 |
WORK( J ) = D( J ) + D( I ) + ETA |
| 235 |
30 CONTINUE |
| 236 |
* |
| 237 |
* Evaluate PSI and the derivative DPSI |
| 238 |
* |
| 239 |
DPSI = ZERO |
| 240 |
PSI = ZERO |
| 241 |
ERRETM = ZERO |
| 242 |
DO 40 J = 1, II |
| 243 |
TEMP = Z( J ) / ( DELTA( J )*WORK( J ) ) |
| 244 |
PSI = PSI + Z( J )*TEMP |
| 245 |
DPSI = DPSI + TEMP*TEMP |
| 246 |
ERRETM = ERRETM + PSI |
| 247 |
40 CONTINUE |
| 248 |
ERRETM = ABS( ERRETM ) |
| 249 |
* |
| 250 |
* Evaluate PHI and the derivative DPHI |
| 251 |
* |
| 252 |
TEMP = Z( N ) / ( DELTA( N )*WORK( N ) ) |
| 253 |
PHI = Z( N )*TEMP |
| 254 |
DPHI = TEMP*TEMP |
| 255 |
ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + |
| 256 |
$ ABS( TAU )*( DPSI+DPHI ) |
| 257 |
* |
| 258 |
W = RHOINV + PHI + PSI |
| 259 |
* |
| 260 |
* Test for convergence |
| 261 |
* |
| 262 |
IF( ABS( W ).LE.EPS*ERRETM ) THEN |
| 263 |
GO TO 240 |
| 264 |
END IF |
| 265 |
* |
| 266 |
* Calculate the new step |
| 267 |
* |
| 268 |
NITER = NITER + 1 |
| 269 |
DTNSQ1 = WORK( N-1 )*DELTA( N-1 ) |
| 270 |
DTNSQ = WORK( N )*DELTA( N ) |
| 271 |
C = W - DTNSQ1*DPSI - DTNSQ*DPHI |
| 272 |
A = ( DTNSQ+DTNSQ1 )*W - DTNSQ*DTNSQ1*( DPSI+DPHI ) |
| 273 |
B = DTNSQ*DTNSQ1*W |
| 274 |
IF( C.LT.ZERO ) |
| 275 |
$ C = ABS( C ) |
| 276 |
IF( C.EQ.ZERO ) THEN |
| 277 |
ETA = RHO - SIGMA*SIGMA |
| 278 |
ELSE IF( A.GE.ZERO ) THEN |
| 279 |
ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) |
| 280 |
ELSE |
| 281 |
ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) |
| 282 |
END IF |
| 283 |
* |
| 284 |
* Note, eta should be positive if w is negative, and |
| 285 |
* eta should be negative otherwise. However, |
| 286 |
* if for some reason caused by roundoff, eta*w > 0, |
| 287 |
* we simply use one Newton step instead. This way |
| 288 |
* will guarantee eta*w < 0. |
| 289 |
* |
| 290 |
IF( W*ETA.GT.ZERO ) |
| 291 |
$ ETA = -W / ( DPSI+DPHI ) |
| 292 |
TEMP = ETA - DTNSQ |
| 293 |
IF( TEMP.GT.RHO ) |
| 294 |
$ ETA = RHO + DTNSQ |
| 295 |
* |
| 296 |
TAU = TAU + ETA |
| 297 |
ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) ) |
| 298 |
DO 50 J = 1, N |
| 299 |
DELTA( J ) = DELTA( J ) - ETA |
| 300 |
WORK( J ) = WORK( J ) + ETA |
| 301 |
50 CONTINUE |
| 302 |
* |
| 303 |
SIGMA = SIGMA + ETA |
| 304 |
* |
| 305 |
* Evaluate PSI and the derivative DPSI |
| 306 |
* |
| 307 |
DPSI = ZERO |
| 308 |
PSI = ZERO |
| 309 |
ERRETM = ZERO |
| 310 |
DO 60 J = 1, II |
| 311 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 312 |
PSI = PSI + Z( J )*TEMP |
| 313 |
DPSI = DPSI + TEMP*TEMP |
| 314 |
ERRETM = ERRETM + PSI |
| 315 |
60 CONTINUE |
| 316 |
ERRETM = ABS( ERRETM ) |
| 317 |
* |
| 318 |
* Evaluate PHI and the derivative DPHI |
| 319 |
* |
| 320 |
TEMP = Z( N ) / ( WORK( N )*DELTA( N ) ) |
| 321 |
PHI = Z( N )*TEMP |
| 322 |
DPHI = TEMP*TEMP |
| 323 |
ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + |
| 324 |
$ ABS( TAU )*( DPSI+DPHI ) |
| 325 |
* |
| 326 |
W = RHOINV + PHI + PSI |
| 327 |
* |
| 328 |
* Main loop to update the values of the array DELTA |
| 329 |
* |
| 330 |
ITER = NITER + 1 |
| 331 |
* |
| 332 |
DO 90 NITER = ITER, MAXIT |
| 333 |
* |
| 334 |
* Test for convergence |
| 335 |
* |
| 336 |
IF( ABS( W ).LE.EPS*ERRETM ) THEN |
| 337 |
GO TO 240 |
| 338 |
END IF |
| 339 |
* |
| 340 |
* Calculate the new step |
| 341 |
* |
| 342 |
DTNSQ1 = WORK( N-1 )*DELTA( N-1 ) |
| 343 |
DTNSQ = WORK( N )*DELTA( N ) |
| 344 |
C = W - DTNSQ1*DPSI - DTNSQ*DPHI |
| 345 |
A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI ) |
| 346 |
B = DTNSQ1*DTNSQ*W |
| 347 |
IF( A.GE.ZERO ) THEN |
| 348 |
ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) |
| 349 |
ELSE |
| 350 |
ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) |
| 351 |
END IF |
| 352 |
* |
| 353 |
* Note, eta should be positive if w is negative, and |
| 354 |
* eta should be negative otherwise. However, |
| 355 |
* if for some reason caused by roundoff, eta*w > 0, |
| 356 |
* we simply use one Newton step instead. This way |
| 357 |
* will guarantee eta*w < 0. |
| 358 |
* |
| 359 |
IF( W*ETA.GT.ZERO ) |
| 360 |
$ ETA = -W / ( DPSI+DPHI ) |
| 361 |
TEMP = ETA - DTNSQ |
| 362 |
IF( TEMP.LE.ZERO ) |
| 363 |
$ ETA = ETA / TWO |
| 364 |
* |
| 365 |
TAU = TAU + ETA |
| 366 |
ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) ) |
| 367 |
DO 70 J = 1, N |
| 368 |
DELTA( J ) = DELTA( J ) - ETA |
| 369 |
WORK( J ) = WORK( J ) + ETA |
| 370 |
70 CONTINUE |
| 371 |
* |
| 372 |
SIGMA = SIGMA + ETA |
| 373 |
* |
| 374 |
* Evaluate PSI and the derivative DPSI |
| 375 |
* |
| 376 |
DPSI = ZERO |
| 377 |
PSI = ZERO |
| 378 |
ERRETM = ZERO |
| 379 |
DO 80 J = 1, II |
| 380 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 381 |
PSI = PSI + Z( J )*TEMP |
| 382 |
DPSI = DPSI + TEMP*TEMP |
| 383 |
ERRETM = ERRETM + PSI |
| 384 |
80 CONTINUE |
| 385 |
ERRETM = ABS( ERRETM ) |
| 386 |
* |
| 387 |
* Evaluate PHI and the derivative DPHI |
| 388 |
* |
| 389 |
TEMP = Z( N ) / ( WORK( N )*DELTA( N ) ) |
| 390 |
PHI = Z( N )*TEMP |
| 391 |
DPHI = TEMP*TEMP |
| 392 |
ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV + |
| 393 |
$ ABS( TAU )*( DPSI+DPHI ) |
| 394 |
* |
| 395 |
W = RHOINV + PHI + PSI |
| 396 |
90 CONTINUE |
| 397 |
* |
| 398 |
* Return with INFO = 1, NITER = MAXIT and not converged |
| 399 |
* |
| 400 |
INFO = 1 |
| 401 |
GO TO 240 |
| 402 |
* |
| 403 |
* End for the case I = N |
| 404 |
* |
| 405 |
ELSE |
| 406 |
* |
| 407 |
* The case for I < N |
| 408 |
* |
| 409 |
NITER = 1 |
| 410 |
IP1 = I + 1 |
| 411 |
* |
| 412 |
* Calculate initial guess |
| 413 |
* |
| 414 |
DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) ) |
| 415 |
DELSQ2 = DELSQ / TWO |
| 416 |
TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) ) |
| 417 |
DO 100 J = 1, N |
| 418 |
WORK( J ) = D( J ) + D( I ) + TEMP |
| 419 |
DELTA( J ) = ( D( J )-D( I ) ) - TEMP |
| 420 |
100 CONTINUE |
| 421 |
* |
| 422 |
PSI = ZERO |
| 423 |
DO 110 J = 1, I - 1 |
| 424 |
PSI = PSI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 425 |
110 CONTINUE |
| 426 |
* |
| 427 |
PHI = ZERO |
| 428 |
DO 120 J = N, I + 2, -1 |
| 429 |
PHI = PHI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 430 |
120 CONTINUE |
| 431 |
C = RHOINV + PSI + PHI |
| 432 |
W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) + |
| 433 |
$ Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) ) |
| 434 |
* |
| 435 |
IF( W.GT.ZERO ) THEN |
| 436 |
* |
| 437 |
* d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2 |
| 438 |
* |
| 439 |
* We choose d(i) as origin. |
| 440 |
* |
| 441 |
ORGATI = .TRUE. |
| 442 |
SG2LB = ZERO |
| 443 |
SG2UB = DELSQ2 |
| 444 |
A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 ) |
| 445 |
B = Z( I )*Z( I )*DELSQ |
| 446 |
IF( A.GT.ZERO ) THEN |
| 447 |
TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) |
| 448 |
ELSE |
| 449 |
TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) |
| 450 |
END IF |
| 451 |
* |
| 452 |
* TAU now is an estimation of SIGMA^2 - D( I )^2. The |
| 453 |
* following, however, is the corresponding estimation of |
| 454 |
* SIGMA - D( I ). |
| 455 |
* |
| 456 |
ETA = TAU / ( D( I )+SQRT( D( I )*D( I )+TAU ) ) |
| 457 |
ELSE |
| 458 |
* |
| 459 |
* (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2 |
| 460 |
* |
| 461 |
* We choose d(i+1) as origin. |
| 462 |
* |
| 463 |
ORGATI = .FALSE. |
| 464 |
SG2LB = -DELSQ2 |
| 465 |
SG2UB = ZERO |
| 466 |
A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 ) |
| 467 |
B = Z( IP1 )*Z( IP1 )*DELSQ |
| 468 |
IF( A.LT.ZERO ) THEN |
| 469 |
TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) ) |
| 470 |
ELSE |
| 471 |
TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C ) |
| 472 |
END IF |
| 473 |
* |
| 474 |
* TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The |
| 475 |
* following, however, is the corresponding estimation of |
| 476 |
* SIGMA - D( IP1 ). |
| 477 |
* |
| 478 |
ETA = TAU / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+ |
| 479 |
$ TAU ) ) ) |
| 480 |
END IF |
| 481 |
* |
| 482 |
IF( ORGATI ) THEN |
| 483 |
II = I |
| 484 |
SIGMA = D( I ) + ETA |
| 485 |
DO 130 J = 1, N |
| 486 |
WORK( J ) = D( J ) + D( I ) + ETA |
| 487 |
DELTA( J ) = ( D( J )-D( I ) ) - ETA |
| 488 |
130 CONTINUE |
| 489 |
ELSE |
| 490 |
II = I + 1 |
| 491 |
SIGMA = D( IP1 ) + ETA |
| 492 |
DO 140 J = 1, N |
| 493 |
WORK( J ) = D( J ) + D( IP1 ) + ETA |
| 494 |
DELTA( J ) = ( D( J )-D( IP1 ) ) - ETA |
| 495 |
140 CONTINUE |
| 496 |
END IF |
| 497 |
IIM1 = II - 1 |
| 498 |
IIP1 = II + 1 |
| 499 |
* |
| 500 |
* Evaluate PSI and the derivative DPSI |
| 501 |
* |
| 502 |
DPSI = ZERO |
| 503 |
PSI = ZERO |
| 504 |
ERRETM = ZERO |
| 505 |
DO 150 J = 1, IIM1 |
| 506 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 507 |
PSI = PSI + Z( J )*TEMP |
| 508 |
DPSI = DPSI + TEMP*TEMP |
| 509 |
ERRETM = ERRETM + PSI |
| 510 |
150 CONTINUE |
| 511 |
ERRETM = ABS( ERRETM ) |
| 512 |
* |
| 513 |
* Evaluate PHI and the derivative DPHI |
| 514 |
* |
| 515 |
DPHI = ZERO |
| 516 |
PHI = ZERO |
| 517 |
DO 160 J = N, IIP1, -1 |
| 518 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 519 |
PHI = PHI + Z( J )*TEMP |
| 520 |
DPHI = DPHI + TEMP*TEMP |
| 521 |
ERRETM = ERRETM + PHI |
| 522 |
160 CONTINUE |
| 523 |
* |
| 524 |
W = RHOINV + PHI + PSI |
| 525 |
* |
| 526 |
* W is the value of the secular function with |
| 527 |
* its ii-th element removed. |
| 528 |
* |
| 529 |
SWTCH3 = .FALSE. |
| 530 |
IF( ORGATI ) THEN |
| 531 |
IF( W.LT.ZERO ) |
| 532 |
$ SWTCH3 = .TRUE. |
| 533 |
ELSE |
| 534 |
IF( W.GT.ZERO ) |
| 535 |
$ SWTCH3 = .TRUE. |
| 536 |
END IF |
| 537 |
IF( II.EQ.1 .OR. II.EQ.N ) |
| 538 |
$ SWTCH3 = .FALSE. |
| 539 |
* |
| 540 |
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) |
| 541 |
DW = DPSI + DPHI + TEMP*TEMP |
| 542 |
TEMP = Z( II )*TEMP |
| 543 |
W = W + TEMP |
| 544 |
ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + |
| 545 |
$ THREE*ABS( TEMP ) + ABS( TAU )*DW |
| 546 |
* |
| 547 |
* Test for convergence |
| 548 |
* |
| 549 |
IF( ABS( W ).LE.EPS*ERRETM ) THEN |
| 550 |
GO TO 240 |
| 551 |
END IF |
| 552 |
* |
| 553 |
IF( W.LE.ZERO ) THEN |
| 554 |
SG2LB = MAX( SG2LB, TAU ) |
| 555 |
ELSE |
| 556 |
SG2UB = MIN( SG2UB, TAU ) |
| 557 |
END IF |
| 558 |
* |
| 559 |
* Calculate the new step |
| 560 |
* |
| 561 |
NITER = NITER + 1 |
| 562 |
IF( .NOT.SWTCH3 ) THEN |
| 563 |
DTIPSQ = WORK( IP1 )*DELTA( IP1 ) |
| 564 |
DTISQ = WORK( I )*DELTA( I ) |
| 565 |
IF( ORGATI ) THEN |
| 566 |
C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2 |
| 567 |
ELSE |
| 568 |
C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2 |
| 569 |
END IF |
| 570 |
A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW |
| 571 |
B = DTIPSQ*DTISQ*W |
| 572 |
IF( C.EQ.ZERO ) THEN |
| 573 |
IF( A.EQ.ZERO ) THEN |
| 574 |
IF( ORGATI ) THEN |
| 575 |
A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI ) |
| 576 |
ELSE |
| 577 |
A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI ) |
| 578 |
END IF |
| 579 |
END IF |
| 580 |
ETA = B / A |
| 581 |
ELSE IF( A.LE.ZERO ) THEN |
| 582 |
ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) |
| 583 |
ELSE |
| 584 |
ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) |
| 585 |
END IF |
| 586 |
ELSE |
| 587 |
* |
| 588 |
* Interpolation using THREE most relevant poles |
| 589 |
* |
| 590 |
DTIIM = WORK( IIM1 )*DELTA( IIM1 ) |
| 591 |
DTIIP = WORK( IIP1 )*DELTA( IIP1 ) |
| 592 |
TEMP = RHOINV + PSI + PHI |
| 593 |
IF( ORGATI ) THEN |
| 594 |
TEMP1 = Z( IIM1 ) / DTIIM |
| 595 |
TEMP1 = TEMP1*TEMP1 |
| 596 |
C = ( TEMP - DTIIP*( DPSI+DPHI ) ) - |
| 597 |
$ ( D( IIM1 )-D( IIP1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1 |
| 598 |
ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) |
| 599 |
IF( DPSI.LT.TEMP1 ) THEN |
| 600 |
ZZ( 3 ) = DTIIP*DTIIP*DPHI |
| 601 |
ELSE |
| 602 |
ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI ) |
| 603 |
END IF |
| 604 |
ELSE |
| 605 |
TEMP1 = Z( IIP1 ) / DTIIP |
| 606 |
TEMP1 = TEMP1*TEMP1 |
| 607 |
C = ( TEMP - DTIIM*( DPSI+DPHI ) ) - |
| 608 |
$ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1 |
| 609 |
IF( DPHI.LT.TEMP1 ) THEN |
| 610 |
ZZ( 1 ) = DTIIM*DTIIM*DPSI |
| 611 |
ELSE |
| 612 |
ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) ) |
| 613 |
END IF |
| 614 |
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) |
| 615 |
END IF |
| 616 |
ZZ( 2 ) = Z( II )*Z( II ) |
| 617 |
DD( 1 ) = DTIIM |
| 618 |
DD( 2 ) = DELTA( II )*WORK( II ) |
| 619 |
DD( 3 ) = DTIIP |
| 620 |
CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO ) |
| 621 |
IF( INFO.NE.0 ) |
| 622 |
$ GO TO 240 |
| 623 |
END IF |
| 624 |
* |
| 625 |
* Note, eta should be positive if w is negative, and |
| 626 |
* eta should be negative otherwise. However, |
| 627 |
* if for some reason caused by roundoff, eta*w > 0, |
| 628 |
* we simply use one Newton step instead. This way |
| 629 |
* will guarantee eta*w < 0. |
| 630 |
* |
| 631 |
IF( W*ETA.GE.ZERO ) |
| 632 |
$ ETA = -W / DW |
| 633 |
IF( ORGATI ) THEN |
| 634 |
TEMP1 = WORK( I )*DELTA( I ) |
| 635 |
TEMP = ETA - TEMP1 |
| 636 |
ELSE |
| 637 |
TEMP1 = WORK( IP1 )*DELTA( IP1 ) |
| 638 |
TEMP = ETA - TEMP1 |
| 639 |
END IF |
| 640 |
IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN |
| 641 |
IF( W.LT.ZERO ) THEN |
| 642 |
ETA = ( SG2UB-TAU ) / TWO |
| 643 |
ELSE |
| 644 |
ETA = ( SG2LB-TAU ) / TWO |
| 645 |
END IF |
| 646 |
END IF |
| 647 |
* |
| 648 |
TAU = TAU + ETA |
| 649 |
ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) ) |
| 650 |
* |
| 651 |
PREW = W |
| 652 |
* |
| 653 |
SIGMA = SIGMA + ETA |
| 654 |
DO 170 J = 1, N |
| 655 |
WORK( J ) = WORK( J ) + ETA |
| 656 |
DELTA( J ) = DELTA( J ) - ETA |
| 657 |
170 CONTINUE |
| 658 |
* |
| 659 |
* Evaluate PSI and the derivative DPSI |
| 660 |
* |
| 661 |
DPSI = ZERO |
| 662 |
PSI = ZERO |
| 663 |
ERRETM = ZERO |
| 664 |
DO 180 J = 1, IIM1 |
| 665 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 666 |
PSI = PSI + Z( J )*TEMP |
| 667 |
DPSI = DPSI + TEMP*TEMP |
| 668 |
ERRETM = ERRETM + PSI |
| 669 |
180 CONTINUE |
| 670 |
ERRETM = ABS( ERRETM ) |
| 671 |
* |
| 672 |
* Evaluate PHI and the derivative DPHI |
| 673 |
* |
| 674 |
DPHI = ZERO |
| 675 |
PHI = ZERO |
| 676 |
DO 190 J = N, IIP1, -1 |
| 677 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 678 |
PHI = PHI + Z( J )*TEMP |
| 679 |
DPHI = DPHI + TEMP*TEMP |
| 680 |
ERRETM = ERRETM + PHI |
| 681 |
190 CONTINUE |
| 682 |
* |
| 683 |
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) |
| 684 |
DW = DPSI + DPHI + TEMP*TEMP |
| 685 |
TEMP = Z( II )*TEMP |
| 686 |
W = RHOINV + PHI + PSI + TEMP |
| 687 |
ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + |
| 688 |
$ THREE*ABS( TEMP ) + ABS( TAU )*DW |
| 689 |
* |
| 690 |
IF( W.LE.ZERO ) THEN |
| 691 |
SG2LB = MAX( SG2LB, TAU ) |
| 692 |
ELSE |
| 693 |
SG2UB = MIN( SG2UB, TAU ) |
| 694 |
END IF |
| 695 |
* |
| 696 |
SWTCH = .FALSE. |
| 697 |
IF( ORGATI ) THEN |
| 698 |
IF( -W.GT.ABS( PREW ) / TEN ) |
| 699 |
$ SWTCH = .TRUE. |
| 700 |
ELSE |
| 701 |
IF( W.GT.ABS( PREW ) / TEN ) |
| 702 |
$ SWTCH = .TRUE. |
| 703 |
END IF |
| 704 |
* |
| 705 |
* Main loop to update the values of the array DELTA and WORK |
| 706 |
* |
| 707 |
ITER = NITER + 1 |
| 708 |
* |
| 709 |
DO 230 NITER = ITER, MAXIT |
| 710 |
* |
| 711 |
* Test for convergence |
| 712 |
* |
| 713 |
IF( ABS( W ).LE.EPS*ERRETM ) THEN |
| 714 |
GO TO 240 |
| 715 |
END IF |
| 716 |
* |
| 717 |
* Calculate the new step |
| 718 |
* |
| 719 |
IF( .NOT.SWTCH3 ) THEN |
| 720 |
DTIPSQ = WORK( IP1 )*DELTA( IP1 ) |
| 721 |
DTISQ = WORK( I )*DELTA( I ) |
| 722 |
IF( .NOT.SWTCH ) THEN |
| 723 |
IF( ORGATI ) THEN |
| 724 |
C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2 |
| 725 |
ELSE |
| 726 |
C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2 |
| 727 |
END IF |
| 728 |
ELSE |
| 729 |
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) |
| 730 |
IF( ORGATI ) THEN |
| 731 |
DPSI = DPSI + TEMP*TEMP |
| 732 |
ELSE |
| 733 |
DPHI = DPHI + TEMP*TEMP |
| 734 |
END IF |
| 735 |
C = W - DTISQ*DPSI - DTIPSQ*DPHI |
| 736 |
END IF |
| 737 |
A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW |
| 738 |
B = DTIPSQ*DTISQ*W |
| 739 |
IF( C.EQ.ZERO ) THEN |
| 740 |
IF( A.EQ.ZERO ) THEN |
| 741 |
IF( .NOT.SWTCH ) THEN |
| 742 |
IF( ORGATI ) THEN |
| 743 |
A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ* |
| 744 |
$ ( DPSI+DPHI ) |
| 745 |
ELSE |
| 746 |
A = Z( IP1 )*Z( IP1 ) + |
| 747 |
$ DTISQ*DTISQ*( DPSI+DPHI ) |
| 748 |
END IF |
| 749 |
ELSE |
| 750 |
A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI |
| 751 |
END IF |
| 752 |
END IF |
| 753 |
ETA = B / A |
| 754 |
ELSE IF( A.LE.ZERO ) THEN |
| 755 |
ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C ) |
| 756 |
ELSE |
| 757 |
ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) |
| 758 |
END IF |
| 759 |
ELSE |
| 760 |
* |
| 761 |
* Interpolation using THREE most relevant poles |
| 762 |
* |
| 763 |
DTIIM = WORK( IIM1 )*DELTA( IIM1 ) |
| 764 |
DTIIP = WORK( IIP1 )*DELTA( IIP1 ) |
| 765 |
TEMP = RHOINV + PSI + PHI |
| 766 |
IF( SWTCH ) THEN |
| 767 |
C = TEMP - DTIIM*DPSI - DTIIP*DPHI |
| 768 |
ZZ( 1 ) = DTIIM*DTIIM*DPSI |
| 769 |
ZZ( 3 ) = DTIIP*DTIIP*DPHI |
| 770 |
ELSE |
| 771 |
IF( ORGATI ) THEN |
| 772 |
TEMP1 = Z( IIM1 ) / DTIIM |
| 773 |
TEMP1 = TEMP1*TEMP1 |
| 774 |
TEMP2 = ( D( IIM1 )-D( IIP1 ) )* |
| 775 |
$ ( D( IIM1 )+D( IIP1 ) )*TEMP1 |
| 776 |
C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2 |
| 777 |
ZZ( 1 ) = Z( IIM1 )*Z( IIM1 ) |
| 778 |
IF( DPSI.LT.TEMP1 ) THEN |
| 779 |
ZZ( 3 ) = DTIIP*DTIIP*DPHI |
| 780 |
ELSE |
| 781 |
ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI ) |
| 782 |
END IF |
| 783 |
ELSE |
| 784 |
TEMP1 = Z( IIP1 ) / DTIIP |
| 785 |
TEMP1 = TEMP1*TEMP1 |
| 786 |
TEMP2 = ( D( IIP1 )-D( IIM1 ) )* |
| 787 |
$ ( D( IIM1 )+D( IIP1 ) )*TEMP1 |
| 788 |
C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2 |
| 789 |
IF( DPHI.LT.TEMP1 ) THEN |
| 790 |
ZZ( 1 ) = DTIIM*DTIIM*DPSI |
| 791 |
ELSE |
| 792 |
ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) ) |
| 793 |
END IF |
| 794 |
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 ) |
| 795 |
END IF |
| 796 |
END IF |
| 797 |
DD( 1 ) = DTIIM |
| 798 |
DD( 2 ) = DELTA( II )*WORK( II ) |
| 799 |
DD( 3 ) = DTIIP |
| 800 |
CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO ) |
| 801 |
IF( INFO.NE.0 ) |
| 802 |
$ GO TO 240 |
| 803 |
END IF |
| 804 |
* |
| 805 |
* Note, eta should be positive if w is negative, and |
| 806 |
* eta should be negative otherwise. However, |
| 807 |
* if for some reason caused by roundoff, eta*w > 0, |
| 808 |
* we simply use one Newton step instead. This way |
| 809 |
* will guarantee eta*w < 0. |
| 810 |
* |
| 811 |
IF( W*ETA.GE.ZERO ) |
| 812 |
$ ETA = -W / DW |
| 813 |
IF( ORGATI ) THEN |
| 814 |
TEMP1 = WORK( I )*DELTA( I ) |
| 815 |
TEMP = ETA - TEMP1 |
| 816 |
ELSE |
| 817 |
TEMP1 = WORK( IP1 )*DELTA( IP1 ) |
| 818 |
TEMP = ETA - TEMP1 |
| 819 |
END IF |
| 820 |
IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN |
| 821 |
IF( W.LT.ZERO ) THEN |
| 822 |
ETA = ( SG2UB-TAU ) / TWO |
| 823 |
ELSE |
| 824 |
ETA = ( SG2LB-TAU ) / TWO |
| 825 |
END IF |
| 826 |
END IF |
| 827 |
* |
| 828 |
TAU = TAU + ETA |
| 829 |
ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) ) |
| 830 |
* |
| 831 |
SIGMA = SIGMA + ETA |
| 832 |
DO 200 J = 1, N |
| 833 |
WORK( J ) = WORK( J ) + ETA |
| 834 |
DELTA( J ) = DELTA( J ) - ETA |
| 835 |
200 CONTINUE |
| 836 |
* |
| 837 |
PREW = W |
| 838 |
* |
| 839 |
* Evaluate PSI and the derivative DPSI |
| 840 |
* |
| 841 |
DPSI = ZERO |
| 842 |
PSI = ZERO |
| 843 |
ERRETM = ZERO |
| 844 |
DO 210 J = 1, IIM1 |
| 845 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 846 |
PSI = PSI + Z( J )*TEMP |
| 847 |
DPSI = DPSI + TEMP*TEMP |
| 848 |
ERRETM = ERRETM + PSI |
| 849 |
210 CONTINUE |
| 850 |
ERRETM = ABS( ERRETM ) |
| 851 |
* |
| 852 |
* Evaluate PHI and the derivative DPHI |
| 853 |
* |
| 854 |
DPHI = ZERO |
| 855 |
PHI = ZERO |
| 856 |
DO 220 J = N, IIP1, -1 |
| 857 |
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) ) |
| 858 |
PHI = PHI + Z( J )*TEMP |
| 859 |
DPHI = DPHI + TEMP*TEMP |
| 860 |
ERRETM = ERRETM + PHI |
| 861 |
220 CONTINUE |
| 862 |
* |
| 863 |
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) ) |
| 864 |
DW = DPSI + DPHI + TEMP*TEMP |
| 865 |
TEMP = Z( II )*TEMP |
| 866 |
W = RHOINV + PHI + PSI + TEMP |
| 867 |
ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV + |
| 868 |
$ THREE*ABS( TEMP ) + ABS( TAU )*DW |
| 869 |
IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN ) |
| 870 |
$ SWTCH = .NOT.SWTCH |
| 871 |
* |
| 872 |
IF( W.LE.ZERO ) THEN |
| 873 |
SG2LB = MAX( SG2LB, TAU ) |
| 874 |
ELSE |
| 875 |
SG2UB = MIN( SG2UB, TAU ) |
| 876 |
END IF |
| 877 |
* |
| 878 |
230 CONTINUE |
| 879 |
* |
| 880 |
* Return with INFO = 1, NITER = MAXIT and not converged |
| 881 |
* |
| 882 |
INFO = 1 |
| 883 |
* |
| 884 |
END IF |
| 885 |
* |
| 886 |
240 CONTINUE |
| 887 |
RETURN |
| 888 |
* |
| 889 |
* End of DLASD4 |
| 890 |
* |
| 891 |
END |