001:       SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            I, INFO, N
010:       DOUBLE PRECISION   DLAM, RHO
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   D( * ), DELTA( * ), Z( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  This subroutine computes the I-th updated eigenvalue of a symmetric
020: *  rank-one modification to a diagonal matrix whose elements are
021: *  given in the array d, and that
022: *
023: *             D(i) < D(j)  for  i < j
024: *
025: *  and that RHO > 0.  This is arranged by the calling routine, and is
026: *  no loss in generality.  The rank-one modified system is thus
027: *
028: *             diag( D )  +  RHO *  Z * Z_transpose.
029: *
030: *  where we assume the Euclidean norm of Z is 1.
031: *
032: *  The method consists of approximating the rational functions in the
033: *  secular equation by simpler interpolating rational functions.
034: *
035: *  Arguments
036: *  =========
037: *
038: *  N      (input) INTEGER
039: *         The length of all arrays.
040: *
041: *  I      (input) INTEGER
042: *         The index of the eigenvalue to be computed.  1 <= I <= N.
043: *
044: *  D      (input) DOUBLE PRECISION array, dimension (N)
045: *         The original eigenvalues.  It is assumed that they are in
046: *         order, D(I) < D(J)  for I < J.
047: *
048: *  Z      (input) DOUBLE PRECISION array, dimension (N)
049: *         The components of the updating vector.
050: *
051: *  DELTA  (output) DOUBLE PRECISION array, dimension (N)
052: *         If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th
053: *         component.  If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5
054: *         for detail. The vector DELTA contains the information necessary
055: *         to construct the eigenvectors by DLAED3 and DLAED9.
056: *
057: *  RHO    (input) DOUBLE PRECISION
058: *         The scalar in the symmetric updating formula.
059: *
060: *  DLAM   (output) DOUBLE PRECISION
061: *         The computed lambda_I, the I-th updated eigenvalue.
062: *
063: *  INFO   (output) INTEGER
064: *         = 0:  successful exit
065: *         > 0:  if INFO = 1, the updating process failed.
066: *
067: *  Internal Parameters
068: *  ===================
069: *
070: *  Logical variable ORGATI (origin-at-i?) is used for distinguishing
071: *  whether D(i) or D(i+1) is treated as the origin.
072: *
073: *            ORGATI = .true.    origin at i
074: *            ORGATI = .false.   origin at i+1
075: *
076: *   Logical variable SWTCH3 (switch-for-3-poles?) is for noting
077: *   if we are working with THREE poles!
078: *
079: *   MAXIT is the maximum number of iterations allowed for each
080: *   eigenvalue.
081: *
082: *  Further Details
083: *  ===============
084: *
085: *  Based on contributions by
086: *     Ren-Cang Li, Computer Science Division, University of California
087: *     at Berkeley, USA
088: *
089: *  =====================================================================
090: *
091: *     .. Parameters ..
092:       INTEGER            MAXIT
093:       PARAMETER          ( MAXIT = 30 )
094:       DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
095:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
096:      $                   THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0,
097:      $                   TEN = 10.0D0 )
098: *     ..
099: *     .. 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
846: