SUBROUTINE DLASD4( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
*
* -- LAPACK auxiliary routine (version 3.1) --
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
* November 2006
*
* .. Scalar Arguments ..
INTEGER I, INFO, N
DOUBLE PRECISION RHO, SIGMA
* ..
* .. Array Arguments ..
DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), Z( * )
* ..
*
* Purpose
* =======
*
* This subroutine computes the square root of the I-th updated
* eigenvalue of a positive symmetric rank-one modification to
* a positive diagonal matrix whose entries are given as the squares
* of the corresponding entries in the array d, and that
*
* 0 <= D(i) < D(j) for i < j
*
* and that RHO > 0. This is arranged by the calling routine, and is
* no loss in generality. The rank-one modified system is thus
*
* diag( D ) * diag( D ) + RHO * Z * Z_transpose.
*
* where we assume the Euclidean norm of Z is 1.
*
* The method consists of approximating the rational functions in the
* secular equation by simpler interpolating rational functions.
*
* Arguments
* =========
*
* N (input) INTEGER
* The length of all arrays.
*
* I (input) INTEGER
* The index of the eigenvalue to be computed. 1 <= I <= N.
*
* D (input) DOUBLE PRECISION array, dimension ( N )
* The original eigenvalues. It is assumed that they are in
* order, 0 <= D(I) < D(J) for I < J.
*
* Z (input) DOUBLE PRECISION array, dimension ( N )
* The components of the updating vector.
*
* DELTA (output) DOUBLE PRECISION array, dimension ( N )
* If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th
* component. If N = 1, then DELTA(1) = 1. The vector DELTA
* contains the information necessary to construct the
* (singular) eigenvectors.
*
* RHO (input) DOUBLE PRECISION
* The scalar in the symmetric updating formula.
*
* SIGMA (output) DOUBLE PRECISION
* The computed sigma_I, the I-th updated eigenvalue.
*
* WORK (workspace) DOUBLE PRECISION array, dimension ( N )
* If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th
* component. If N = 1, then WORK( 1 ) = 1.
*
* INFO (output) INTEGER
* = 0: successful exit
* > 0: if INFO = 1, the updating process failed.
*
* Internal Parameters
* ===================
*
* Logical variable ORGATI (origin-at-i?) is used for distinguishing
* whether D(i) or D(i+1) is treated as the origin.
*
* ORGATI = .true. origin at i
* ORGATI = .false. origin at i+1
*
* Logical variable SWTCH3 (switch-for-3-poles?) is for noting
* if we are working with THREE poles!
*
* MAXIT is the maximum number of iterations allowed for each
* eigenvalue.
*
* Further Details
* ===============
*
* Based on contributions by
* Ren-Cang Li, Computer Science Division, University of California
* at Berkeley, USA
*
* =====================================================================
*
* .. Parameters ..
INTEGER MAXIT
PARAMETER ( MAXIT = 20 )
DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
$ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
$ TEN = 10.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL ORGATI, SWTCH, SWTCH3
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM,
$ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
$ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB,
$ SG2UB, TAU, TEMP, TEMP1, TEMP2, W
* ..
* .. Local Arrays ..
DOUBLE PRECISION DD( 3 ), ZZ( 3 )
* ..
* .. External Subroutines ..
EXTERNAL DLAED6, DLASD5
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
* ..
* .. Executable Statements ..
*
* Since this routine is called in an inner loop, we do no argument
* checking.
*
* Quick return for N=1 and 2.
*
INFO = 0
IF( N.EQ.1 ) THEN
*
* Presumably, I=1 upon entry
*
SIGMA = SQRT( D( 1 )*D( 1 )+RHO*Z( 1 )*Z( 1 ) )
DELTA( 1 ) = ONE
WORK( 1 ) = ONE
RETURN
END IF
IF( N.EQ.2 ) THEN
CALL DLASD5( I, D, Z, DELTA, RHO, SIGMA, WORK )
RETURN
END IF
*
* Compute machine epsilon
*
EPS = DLAMCH( 'Epsilon' )
RHOINV = ONE / RHO
*
* The case I = N
*
IF( I.EQ.N ) THEN
*
* Initialize some basic variables
*
II = N - 1
NITER = 1
*
* Calculate initial guess
*
TEMP = RHO / TWO
*
* If ||Z||_2 is not one, then TEMP should be set to
* RHO * ||Z||_2^2 / TWO
*
TEMP1 = TEMP / ( D( N )+SQRT( D( N )*D( N )+TEMP ) )
DO 10 J = 1, N
WORK( J ) = D( J ) + D( N ) + TEMP1
DELTA( J ) = ( D( J )-D( N ) ) - TEMP1
10 CONTINUE
*
PSI = ZERO
DO 20 J = 1, N - 2
PSI = PSI + Z( J )*Z( J ) / ( DELTA( J )*WORK( J ) )
20 CONTINUE
*
C = RHOINV + PSI
W = C + Z( II )*Z( II ) / ( DELTA( II )*WORK( II ) ) +
$ Z( N )*Z( N ) / ( DELTA( N )*WORK( N ) )
*
IF( W.LE.ZERO ) THEN
TEMP1 = SQRT( D( N )*D( N )+RHO )
TEMP = Z( N-1 )*Z( N-1 ) / ( ( D( N-1 )+TEMP1 )*
$ ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +
$ Z( N )*Z( N ) / RHO
*
* The following TAU is to approximate
* SIGMA_n^2 - D( N )*D( N )
*
IF( C.LE.TEMP ) THEN
TAU = RHO
ELSE
DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
B = Z( N )*Z( N )*DELSQ
IF( A.LT.ZERO ) THEN
TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
ELSE
TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
END IF
END IF
*
* It can be proved that
* D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO
*
ELSE
DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
B = Z( N )*Z( N )*DELSQ
*
* The following TAU is to approximate
* SIGMA_n^2 - D( N )*D( N )
*
IF( A.LT.ZERO ) THEN
TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
ELSE
TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
END IF
*
* It can be proved that
* D(N)^2 < D(N)^2+TAU < SIGMA(N)^2 < D(N)^2+RHO/2
*
END IF
*
* The following ETA is to approximate SIGMA_n - D( N )
*
ETA = TAU / ( D( N )+SQRT( D( N )*D( N )+TAU ) )
*
SIGMA = D( N ) + ETA
DO 30 J = 1, N
DELTA( J ) = ( D( J )-D( I ) ) - ETA
WORK( J ) = D( J ) + D( I ) + ETA
30 CONTINUE
*
* Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 40 J = 1, II
TEMP = Z( J ) / ( DELTA( J )*WORK( J ) )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
40 CONTINUE
ERRETM = ABS( ERRETM )
*
* Evaluate PHI and the derivative DPHI
*
TEMP = Z( N ) / ( DELTA( N )*WORK( N ) )
PHI = Z( N )*TEMP
DPHI = TEMP*TEMP
ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
$ ABS( TAU )*( DPSI+DPHI )
*
W = RHOINV + PHI + PSI
*
* Test for convergence
*
IF( ABS( W ).LE.EPS*ERRETM ) THEN
GO TO 240
END IF
*
* Calculate the new step
*
NITER = NITER + 1
DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
DTNSQ = WORK( N )*DELTA( N )
C = W - DTNSQ1*DPSI - DTNSQ*DPHI
A = ( DTNSQ+DTNSQ1 )*W - DTNSQ*DTNSQ1*( DPSI+DPHI )
B = DTNSQ*DTNSQ1*W
IF( C.LT.ZERO )
$ C = ABS( C )
IF( C.EQ.ZERO ) THEN
ETA = RHO - SIGMA*SIGMA
ELSE IF( A.GE.ZERO ) THEN
ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE
ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
END IF
*
* Note, eta should be positive if w is negative, and
* eta should be negative otherwise. However,
* if for some reason caused by roundoff, eta*w > 0,
* we simply use one Newton step instead. This way
* will guarantee eta*w < 0.
*
IF( W*ETA.GT.ZERO )
$ ETA = -W / ( DPSI+DPHI )
TEMP = ETA - DTNSQ
IF( TEMP.GT.RHO )
$ ETA = RHO + DTNSQ
*
TAU = TAU + ETA
ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
DO 50 J = 1, N
DELTA( J ) = DELTA( J ) - ETA
WORK( J ) = WORK( J ) + ETA
50 CONTINUE
*
SIGMA = SIGMA + ETA
*
* Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 60 J = 1, II
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
60 CONTINUE
ERRETM = ABS( ERRETM )
*
* Evaluate PHI and the derivative DPHI
*
TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
PHI = Z( N )*TEMP
DPHI = TEMP*TEMP
ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
$ ABS( TAU )*( DPSI+DPHI )
*
W = RHOINV + PHI + PSI
*
* Main loop to update the values of the array DELTA
*
ITER = NITER + 1
*
DO 90 NITER = ITER, MAXIT
*
* Test for convergence
*
IF( ABS( W ).LE.EPS*ERRETM ) THEN
GO TO 240
END IF
*
* Calculate the new step
*
DTNSQ1 = WORK( N-1 )*DELTA( N-1 )
DTNSQ = WORK( N )*DELTA( N )
C = W - DTNSQ1*DPSI - DTNSQ*DPHI
A = ( DTNSQ+DTNSQ1 )*W - DTNSQ1*DTNSQ*( DPSI+DPHI )
B = DTNSQ1*DTNSQ*W
IF( A.GE.ZERO ) THEN
ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE
ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
END IF
*
* Note, eta should be positive if w is negative, and
* eta should be negative otherwise. However,
* if for some reason caused by roundoff, eta*w > 0,
* we simply use one Newton step instead. This way
* will guarantee eta*w < 0.
*
IF( W*ETA.GT.ZERO )
$ ETA = -W / ( DPSI+DPHI )
TEMP = ETA - DTNSQ
IF( TEMP.LE.ZERO )
$ ETA = ETA / TWO
*
TAU = TAU + ETA
ETA = ETA / ( SIGMA+SQRT( ETA+SIGMA*SIGMA ) )
DO 70 J = 1, N
DELTA( J ) = DELTA( J ) - ETA
WORK( J ) = WORK( J ) + ETA
70 CONTINUE
*
SIGMA = SIGMA + ETA
*
* Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 80 J = 1, II
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
80 CONTINUE
ERRETM = ABS( ERRETM )
*
* Evaluate PHI and the derivative DPHI
*
TEMP = Z( N ) / ( WORK( N )*DELTA( N ) )
PHI = Z( N )*TEMP
DPHI = TEMP*TEMP
ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
$ ABS( TAU )*( DPSI+DPHI )
*
W = RHOINV + PHI + PSI
90 CONTINUE
*
* Return with INFO = 1, NITER = MAXIT and not converged
*
INFO = 1
GO TO 240
*
* End for the case I = N
*
ELSE
*
* The case for I < N
*
NITER = 1
IP1 = I + 1
*
* Calculate initial guess
*
DELSQ = ( D( IP1 )-D( I ) )*( D( IP1 )+D( I ) )
DELSQ2 = DELSQ / TWO
TEMP = DELSQ2 / ( D( I )+SQRT( D( I )*D( I )+DELSQ2 ) )
DO 100 J = 1, N
WORK( J ) = D( J ) + D( I ) + TEMP
DELTA( J ) = ( D( J )-D( I ) ) - TEMP
100 CONTINUE
*
PSI = ZERO
DO 110 J = 1, I - 1
PSI = PSI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
110 CONTINUE
*
PHI = ZERO
DO 120 J = N, I + 2, -1
PHI = PHI + Z( J )*Z( J ) / ( WORK( J )*DELTA( J ) )
120 CONTINUE
C = RHOINV + PSI + PHI
W = C + Z( I )*Z( I ) / ( WORK( I )*DELTA( I ) ) +
$ Z( IP1 )*Z( IP1 ) / ( WORK( IP1 )*DELTA( IP1 ) )
*
IF( W.GT.ZERO ) THEN
*
* d(i)^2 < the ith sigma^2 < (d(i)^2+d(i+1)^2)/2
*
* We choose d(i) as origin.
*
ORGATI = .TRUE.
SG2LB = ZERO
SG2UB = DELSQ2
A = C*DELSQ + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
B = Z( I )*Z( I )*DELSQ
IF( A.GT.ZERO ) THEN
TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
ELSE
TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
END IF
*
* TAU now is an estimation of SIGMA^2 - D( I )^2. The
* following, however, is the corresponding estimation of
* SIGMA - D( I ).
*
ETA = TAU / ( D( I )+SQRT( D( I )*D( I )+TAU ) )
ELSE
*
* (d(i)^2+d(i+1)^2)/2 <= the ith sigma^2 < d(i+1)^2/2
*
* We choose d(i+1) as origin.
*
ORGATI = .FALSE.
SG2LB = -DELSQ2
SG2UB = ZERO
A = C*DELSQ - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
B = Z( IP1 )*Z( IP1 )*DELSQ
IF( A.LT.ZERO ) THEN
TAU = TWO*B / ( A-SQRT( ABS( A*A+FOUR*B*C ) ) )
ELSE
TAU = -( A+SQRT( ABS( A*A+FOUR*B*C ) ) ) / ( TWO*C )
END IF
*
* TAU now is an estimation of SIGMA^2 - D( IP1 )^2. The
* following, however, is the corresponding estimation of
* SIGMA - D( IP1 ).
*
ETA = TAU / ( D( IP1 )+SQRT( ABS( D( IP1 )*D( IP1 )+
$ TAU ) ) )
END IF
*
IF( ORGATI ) THEN
II = I
SIGMA = D( I ) + ETA
DO 130 J = 1, N
WORK( J ) = D( J ) + D( I ) + ETA
DELTA( J ) = ( D( J )-D( I ) ) - ETA
130 CONTINUE
ELSE
II = I + 1
SIGMA = D( IP1 ) + ETA
DO 140 J = 1, N
WORK( J ) = D( J ) + D( IP1 ) + ETA
DELTA( J ) = ( D( J )-D( IP1 ) ) - ETA
140 CONTINUE
END IF
IIM1 = II - 1
IIP1 = II + 1
*
* Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 150 J = 1, IIM1
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
150 CONTINUE
ERRETM = ABS( ERRETM )
*
* Evaluate PHI and the derivative DPHI
*
DPHI = ZERO
PHI = ZERO
DO 160 J = N, IIP1, -1
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PHI = PHI + Z( J )*TEMP
DPHI = DPHI + TEMP*TEMP
ERRETM = ERRETM + PHI
160 CONTINUE
*
W = RHOINV + PHI + PSI
*
* W is the value of the secular function with
* its ii-th element removed.
*
SWTCH3 = .FALSE.
IF( ORGATI ) THEN
IF( W.LT.ZERO )
$ SWTCH3 = .TRUE.
ELSE
IF( W.GT.ZERO )
$ SWTCH3 = .TRUE.
END IF
IF( II.EQ.1 .OR. II.EQ.N )
$ SWTCH3 = .FALSE.
*
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
DW = DPSI + DPHI + TEMP*TEMP
TEMP = Z( II )*TEMP
W = W + TEMP
ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
$ THREE*ABS( TEMP ) + ABS( TAU )*DW
*
* Test for convergence
*
IF( ABS( W ).LE.EPS*ERRETM ) THEN
GO TO 240
END IF
*
IF( W.LE.ZERO ) THEN
SG2LB = MAX( SG2LB, TAU )
ELSE
SG2UB = MIN( SG2UB, TAU )
END IF
*
* Calculate the new step
*
NITER = NITER + 1
IF( .NOT.SWTCH3 ) THEN
DTIPSQ = WORK( IP1 )*DELTA( IP1 )
DTISQ = WORK( I )*DELTA( I )
IF( ORGATI ) THEN
C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
ELSE
C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
END IF
A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
B = DTIPSQ*DTISQ*W
IF( C.EQ.ZERO ) THEN
IF( A.EQ.ZERO ) THEN
IF( ORGATI ) THEN
A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*( DPSI+DPHI )
ELSE
A = Z( IP1 )*Z( IP1 ) + DTISQ*DTISQ*( DPSI+DPHI )
END IF
END IF
ETA = B / A
ELSE IF( A.LE.ZERO ) THEN
ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE
ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
END IF
ELSE
*
* Interpolation using THREE most relevant poles
*
DTIIM = WORK( IIM1 )*DELTA( IIM1 )
DTIIP = WORK( IIP1 )*DELTA( IIP1 )
TEMP = RHOINV + PSI + PHI
IF( ORGATI ) THEN
TEMP1 = Z( IIM1 ) / DTIIM
TEMP1 = TEMP1*TEMP1
C = ( TEMP - DTIIP*( DPSI+DPHI ) ) -
$ ( D( IIM1 )-D( IIP1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
IF( DPSI.LT.TEMP1 ) THEN
ZZ( 3 ) = DTIIP*DTIIP*DPHI
ELSE
ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
END IF
ELSE
TEMP1 = Z( IIP1 ) / DTIIP
TEMP1 = TEMP1*TEMP1
C = ( TEMP - DTIIM*( DPSI+DPHI ) ) -
$ ( D( IIP1 )-D( IIM1 ) )*( D( IIM1 )+D( IIP1 ) )*TEMP1
IF( DPHI.LT.TEMP1 ) THEN
ZZ( 1 ) = DTIIM*DTIIM*DPSI
ELSE
ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
END IF
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
END IF
ZZ( 2 ) = Z( II )*Z( II )
DD( 1 ) = DTIIM
DD( 2 ) = DELTA( II )*WORK( II )
DD( 3 ) = DTIIP
CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
IF( INFO.NE.0 )
$ GO TO 240
END IF
*
* Note, eta should be positive if w is negative, and
* eta should be negative otherwise. However,
* if for some reason caused by roundoff, eta*w > 0,
* we simply use one Newton step instead. This way
* will guarantee eta*w < 0.
*
IF( W*ETA.GE.ZERO )
$ ETA = -W / DW
IF( ORGATI ) THEN
TEMP1 = WORK( I )*DELTA( I )
TEMP = ETA - TEMP1
ELSE
TEMP1 = WORK( IP1 )*DELTA( IP1 )
TEMP = ETA - TEMP1
END IF
IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
IF( W.LT.ZERO ) THEN
ETA = ( SG2UB-TAU ) / TWO
ELSE
ETA = ( SG2LB-TAU ) / TWO
END IF
END IF
*
TAU = TAU + ETA
ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
*
PREW = W
*
SIGMA = SIGMA + ETA
DO 170 J = 1, N
WORK( J ) = WORK( J ) + ETA
DELTA( J ) = DELTA( J ) - ETA
170 CONTINUE
*
* Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 180 J = 1, IIM1
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
180 CONTINUE
ERRETM = ABS( ERRETM )
*
* Evaluate PHI and the derivative DPHI
*
DPHI = ZERO
PHI = ZERO
DO 190 J = N, IIP1, -1
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PHI = PHI + Z( J )*TEMP
DPHI = DPHI + TEMP*TEMP
ERRETM = ERRETM + PHI
190 CONTINUE
*
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
DW = DPSI + DPHI + TEMP*TEMP
TEMP = Z( II )*TEMP
W = RHOINV + PHI + PSI + TEMP
ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
$ THREE*ABS( TEMP ) + ABS( TAU )*DW
*
IF( W.LE.ZERO ) THEN
SG2LB = MAX( SG2LB, TAU )
ELSE
SG2UB = MIN( SG2UB, TAU )
END IF
*
SWTCH = .FALSE.
IF( ORGATI ) THEN
IF( -W.GT.ABS( PREW ) / TEN )
$ SWTCH = .TRUE.
ELSE
IF( W.GT.ABS( PREW ) / TEN )
$ SWTCH = .TRUE.
END IF
*
* Main loop to update the values of the array DELTA and WORK
*
ITER = NITER + 1
*
DO 230 NITER = ITER, MAXIT
*
* Test for convergence
*
IF( ABS( W ).LE.EPS*ERRETM ) THEN
GO TO 240
END IF
*
* Calculate the new step
*
IF( .NOT.SWTCH3 ) THEN
DTIPSQ = WORK( IP1 )*DELTA( IP1 )
DTISQ = WORK( I )*DELTA( I )
IF( .NOT.SWTCH ) THEN
IF( ORGATI ) THEN
C = W - DTIPSQ*DW + DELSQ*( Z( I ) / DTISQ )**2
ELSE
C = W - DTISQ*DW - DELSQ*( Z( IP1 ) / DTIPSQ )**2
END IF
ELSE
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
IF( ORGATI ) THEN
DPSI = DPSI + TEMP*TEMP
ELSE
DPHI = DPHI + TEMP*TEMP
END IF
C = W - DTISQ*DPSI - DTIPSQ*DPHI
END IF
A = ( DTIPSQ+DTISQ )*W - DTIPSQ*DTISQ*DW
B = DTIPSQ*DTISQ*W
IF( C.EQ.ZERO ) THEN
IF( A.EQ.ZERO ) THEN
IF( .NOT.SWTCH ) THEN
IF( ORGATI ) THEN
A = Z( I )*Z( I ) + DTIPSQ*DTIPSQ*
$ ( DPSI+DPHI )
ELSE
A = Z( IP1 )*Z( IP1 ) +
$ DTISQ*DTISQ*( DPSI+DPHI )
END IF
ELSE
A = DTISQ*DTISQ*DPSI + DTIPSQ*DTIPSQ*DPHI
END IF
END IF
ETA = B / A
ELSE IF( A.LE.ZERO ) THEN
ETA = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
ELSE
ETA = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
END IF
ELSE
*
* Interpolation using THREE most relevant poles
*
DTIIM = WORK( IIM1 )*DELTA( IIM1 )
DTIIP = WORK( IIP1 )*DELTA( IIP1 )
TEMP = RHOINV + PSI + PHI
IF( SWTCH ) THEN
C = TEMP - DTIIM*DPSI - DTIIP*DPHI
ZZ( 1 ) = DTIIM*DTIIM*DPSI
ZZ( 3 ) = DTIIP*DTIIP*DPHI
ELSE
IF( ORGATI ) THEN
TEMP1 = Z( IIM1 ) / DTIIM
TEMP1 = TEMP1*TEMP1
TEMP2 = ( D( IIM1 )-D( IIP1 ) )*
$ ( D( IIM1 )+D( IIP1 ) )*TEMP1
C = TEMP - DTIIP*( DPSI+DPHI ) - TEMP2
ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
IF( DPSI.LT.TEMP1 ) THEN
ZZ( 3 ) = DTIIP*DTIIP*DPHI
ELSE
ZZ( 3 ) = DTIIP*DTIIP*( ( DPSI-TEMP1 )+DPHI )
END IF
ELSE
TEMP1 = Z( IIP1 ) / DTIIP
TEMP1 = TEMP1*TEMP1
TEMP2 = ( D( IIP1 )-D( IIM1 ) )*
$ ( D( IIM1 )+D( IIP1 ) )*TEMP1
C = TEMP - DTIIM*( DPSI+DPHI ) - TEMP2
IF( DPHI.LT.TEMP1 ) THEN
ZZ( 1 ) = DTIIM*DTIIM*DPSI
ELSE
ZZ( 1 ) = DTIIM*DTIIM*( DPSI+( DPHI-TEMP1 ) )
END IF
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
END IF
END IF
DD( 1 ) = DTIIM
DD( 2 ) = DELTA( II )*WORK( II )
DD( 3 ) = DTIIP
CALL DLAED6( NITER, ORGATI, C, DD, ZZ, W, ETA, INFO )
IF( INFO.NE.0 )
$ GO TO 240
END IF
*
* Note, eta should be positive if w is negative, and
* eta should be negative otherwise. However,
* if for some reason caused by roundoff, eta*w > 0,
* we simply use one Newton step instead. This way
* will guarantee eta*w < 0.
*
IF( W*ETA.GE.ZERO )
$ ETA = -W / DW
IF( ORGATI ) THEN
TEMP1 = WORK( I )*DELTA( I )
TEMP = ETA - TEMP1
ELSE
TEMP1 = WORK( IP1 )*DELTA( IP1 )
TEMP = ETA - TEMP1
END IF
IF( TEMP.GT.SG2UB .OR. TEMP.LT.SG2LB ) THEN
IF( W.LT.ZERO ) THEN
ETA = ( SG2UB-TAU ) / TWO
ELSE
ETA = ( SG2LB-TAU ) / TWO
END IF
END IF
*
TAU = TAU + ETA
ETA = ETA / ( SIGMA+SQRT( SIGMA*SIGMA+ETA ) )
*
SIGMA = SIGMA + ETA
DO 200 J = 1, N
WORK( J ) = WORK( J ) + ETA
DELTA( J ) = DELTA( J ) - ETA
200 CONTINUE
*
PREW = W
*
* Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 210 J = 1, IIM1
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
210 CONTINUE
ERRETM = ABS( ERRETM )
*
* Evaluate PHI and the derivative DPHI
*
DPHI = ZERO
PHI = ZERO
DO 220 J = N, IIP1, -1
TEMP = Z( J ) / ( WORK( J )*DELTA( J ) )
PHI = PHI + Z( J )*TEMP
DPHI = DPHI + TEMP*TEMP
ERRETM = ERRETM + PHI
220 CONTINUE
*
TEMP = Z( II ) / ( WORK( II )*DELTA( II ) )
DW = DPSI + DPHI + TEMP*TEMP
TEMP = Z( II )*TEMP
W = RHOINV + PHI + PSI + TEMP
ERRETM = EIGHT*( PHI-PSI ) + ERRETM + TWO*RHOINV +
$ THREE*ABS( TEMP ) + ABS( TAU )*DW
IF( W*PREW.GT.ZERO .AND. ABS( W ).GT.ABS( PREW ) / TEN )
$ SWTCH = .NOT.SWTCH
*
IF( W.LE.ZERO ) THEN
SG2LB = MAX( SG2LB, TAU )
ELSE
SG2UB = MIN( SG2UB, TAU )
END IF
*
230 CONTINUE
*
* Return with INFO = 1, NITER = MAXIT and not converged
*
INFO = 1
*
END IF
*
240 CONTINUE
RETURN
*
* End of DLASD4
*
END