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