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