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