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