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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed4.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack