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