1: *> \brief \b DLAED4 used by DSTEDC. Finds a single root of the secular equation.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
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">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER I, INFO, N
25: * DOUBLE PRECISION DLAM, RHO
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
29: * ..
30: *
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 > 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: *
130: *> \author Univ. of Tennessee
131: *> \author Univ. of California Berkeley
132: *> \author Univ. of Colorado Denver
133: *> \author NAG Ltd.
134: *
135: *> \ingroup auxOTHERcomputational
136: *
137: *> \par Contributors:
138: * ==================
139: *>
140: *> Ren-Cang Li, Computer Science Division, University of California
141: *> at Berkeley, USA
142: *>
143: * =====================================================================
144: SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
145: *
146: * -- LAPACK computational routine --
147: * -- LAPACK is a software package provided by Univ. of Tennessee, --
148: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149: *
150: * .. Scalar Arguments ..
151: INTEGER I, INFO, N
152: DOUBLE PRECISION DLAM, RHO
153: * ..
154: * .. Array Arguments ..
155: DOUBLE PRECISION D( * ), DELTA( * ), Z( * )
156: * ..
157: *
158: * =====================================================================
159: *
160: * .. Parameters ..
161: INTEGER MAXIT
162: PARAMETER ( MAXIT = 30 )
163: DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
164: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
165: $ THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0,
166: $ TEN = 10.0D0 )
167: * ..
168: * .. Local Scalars ..
169: LOGICAL ORGATI, SWTCH, SWTCH3
170: INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
171: DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
172: $ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
173: $ RHOINV, TAU, TEMP, TEMP1, W
174: * ..
175: * .. Local Arrays ..
176: DOUBLE PRECISION ZZ( 3 )
177: * ..
178: * .. External Functions ..
179: DOUBLE PRECISION DLAMCH
180: EXTERNAL DLAMCH
181: * ..
182: * .. External Subroutines ..
183: EXTERNAL DLAED5, DLAED6
184: * ..
185: * .. Intrinsic Functions ..
186: INTRINSIC ABS, MAX, MIN, SQRT
187: * ..
188: * .. Executable Statements ..
189: *
190: * Since this routine is called in an inner loop, we do no argument
191: * checking.
192: *
193: * Quick return for N=1 and 2.
194: *
195: INFO = 0
196: IF( N.EQ.1 ) THEN
197: *
198: * Presumably, I=1 upon entry
199: *
200: DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 )
201: DELTA( 1 ) = ONE
202: RETURN
203: END IF
204: IF( N.EQ.2 ) THEN
205: CALL DLAED5( I, D, Z, DELTA, RHO, DLAM )
206: RETURN
207: END IF
208: *
209: * Compute machine epsilon
210: *
211: EPS = DLAMCH( 'Epsilon' )
212: RHOINV = ONE / RHO
213: *
214: * The case I = N
215: *
216: IF( I.EQ.N ) THEN
217: *
218: * Initialize some basic variables
219: *
220: II = N - 1
221: NITER = 1
222: *
223: * Calculate initial guess
224: *
225: MIDPT = RHO / TWO
226: *
227: * If ||Z||_2 is not one, then TEMP should be set to
228: * RHO * ||Z||_2^2 / TWO
229: *
230: DO 10 J = 1, N
231: DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
232: 10 CONTINUE
233: *
234: PSI = ZERO
235: DO 20 J = 1, N - 2
236: PSI = PSI + Z( J )*Z( J ) / DELTA( J )
237: 20 CONTINUE
238: *
239: C = RHOINV + PSI
240: W = C + Z( II )*Z( II ) / DELTA( II ) +
241: $ Z( N )*Z( N ) / DELTA( N )
242: *
243: IF( W.LE.ZERO ) THEN
244: TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) +
245: $ Z( N )*Z( N ) / RHO
246: IF( C.LE.TEMP ) THEN
247: TAU = RHO
248: ELSE
249: DEL = D( N ) - D( N-1 )
250: A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
251: B = Z( N )*Z( N )*DEL
252: IF( A.LT.ZERO ) THEN
253: TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
254: ELSE
255: TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
256: END IF
257: END IF
258: *
259: * It can be proved that
260: * D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO
261: *
262: DLTLB = MIDPT
263: DLTUB = RHO
264: ELSE
265: DEL = D( N ) - D( N-1 )
266: A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
267: B = Z( N )*Z( N )*DEL
268: IF( A.LT.ZERO ) THEN
269: TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
270: ELSE
271: TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
272: END IF
273: *
274: * It can be proved that
275: * D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2
276: *
277: DLTLB = ZERO
278: DLTUB = MIDPT
279: END IF
280: *
281: DO 30 J = 1, N
282: DELTA( J ) = ( D( J )-D( I ) ) - TAU
283: 30 CONTINUE
284: *
285: * Evaluate PSI and the derivative DPSI
286: *
287: DPSI = ZERO
288: PSI = ZERO
289: ERRETM = ZERO
290: DO 40 J = 1, II
291: TEMP = Z( J ) / DELTA( J )
292: PSI = PSI + Z( J )*TEMP
293: DPSI = DPSI + TEMP*TEMP
294: ERRETM = ERRETM + PSI
295: 40 CONTINUE
296: ERRETM = ABS( ERRETM )
297: *
298: * Evaluate PHI and the derivative DPHI
299: *
300: TEMP = Z( N ) / DELTA( N )
301: PHI = Z( N )*TEMP
302: DPHI = TEMP*TEMP
303: ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
304: $ ABS( TAU )*( DPSI+DPHI )
305: *
306: W = RHOINV + PHI + PSI
307: *
308: * Test for convergence
309: *
310: IF( ABS( W ).LE.EPS*ERRETM ) THEN
311: DLAM = D( I ) + TAU
312: GO TO 250
313: END IF
314: *
315: IF( W.LE.ZERO ) THEN
316: DLTLB = MAX( DLTLB, TAU )
317: ELSE
318: DLTUB = MIN( DLTUB, TAU )
319: END IF
320: *
321: * Calculate the new step
322: *
323: NITER = NITER + 1
324: C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
325: A = ( DELTA( N-1 )+DELTA( N ) )*W -
326: $ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
327: B = DELTA( N-1 )*DELTA( N )*W
328: IF( C.LT.ZERO )
329: $ C = ABS( C )
330: IF( C.EQ.ZERO ) THEN
331: * ETA = B/A
332: * ETA = RHO - TAU
333: * ETA = DLTUB - TAU
334: *
335: * Update proposed by Li, Ren-Cang:
336: ETA = -W / ( DPSI+DPHI )
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
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