slaed4 - is used by sstedc. Finds a single root of the secular equation
SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) INTEGER I, INFO, N REAL DLAM, RHO REAL D(*),DELTA(*), Z(*) SUBROUTINE SLAED4_64( N, I, D, Z, DELTA, RHO, DLAM, INFO ) INTEGER*8 I, INFO, N REAL DLAM, RHO REAL D(*),DELTA(*), Z(*) F95 INTERFACE SUBROUTINE LAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) INTEGER :: N, I, INFO REAL, DIMENSION(:) :: D, Z, DELTA REAL :: RHO, DLAM SUBROUTINE LAED4_64( N, I, D, Z, DELTA, RHO, DLAM, INFO ) INTEGER(8) :: N, I, INFO REAL, DIMENSION(:) :: D, Z, DELTA REAL :: RHO, DLAM C INTERFACE #include <sunperf.h> void slaed4 (int n, int i, float *d, float *z, float *delta, float rho, float *dlam, int *info); void slaed4_64 (long n, long i, float *d, float *z, float *delta, float rho, float *dlam, long *info);
Oracle Solaris Studio Performance Library slaed4(3P)
NAME
slaed4 - is used by sstedc. Finds a single root of the secular equation
SYNOPSIS
SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
INTEGER I, INFO, N
REAL DLAM, RHO
REAL D(*),DELTA(*), Z(*)
SUBROUTINE SLAED4_64( N, I, D, Z, DELTA, RHO, DLAM, INFO )
INTEGER*8 I, INFO, N
REAL DLAM, RHO
REAL D(*),DELTA(*), Z(*)
F95 INTERFACE
SUBROUTINE LAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
INTEGER :: N, I, INFO
REAL, DIMENSION(:) :: D, Z, DELTA
REAL :: RHO, DLAM
SUBROUTINE LAED4_64( N, I, D, Z, DELTA, RHO, DLAM, INFO )
INTEGER(8) :: N, I, INFO
REAL, DIMENSION(:) :: D, Z, DELTA
REAL :: RHO, DLAM
C INTERFACE
#include <sunperf.h>
void slaed4 (int n, int i, float *d, float *z, float *delta, float rho,
float *dlam, int *info);
void slaed4_64 (long n, long i, float *d, float *z, float *delta, float
rho, float *dlam, long *info);
PURPOSE
SUBROUTINE slaed4( N, I, D, Z, DELTA, RHO, DLAM, INFO ) * * --
LAPACK computational routine (version 3.4.2) -- * -- LAPACK is a soft-
ware package provided by Univ. of Tennessee, -- * -- Univ. of Cali-
fornia Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * Septem-
ber 2012 * * .. Scalar Arguments ..
INTEGER I, INFO, N
REAL DLAM, RHO * .. * .. Array Arguments
..
REAL D( * ), DELTA( * ), Z( * ) * .. * *
===================================================================== *
* .. Parameters ..
INTEGER MAXIT
PARAMETER ( MAXIT = 30 )
REAL ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
$ THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0,
$ TEN = 10.0E0 ) * .. * .. Local
Scalars ..
LOGICAL ORGATI, SWTCH, SWTCH3
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
REAL A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
$ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
$ RHOINV, TAU, TEMP, TEMP1, W * .. * ..
Local Arrays ..
REAL ZZ( 3 ) * .. * .. External Functions
..
REAL SLAMCH
EXTERNAL SLAMCH * .. * .. External Subroutines
..
EXTERNAL SLAED5, SLAED6 * .. * .. Intrinsic
Functions ..
INTRINSIC ABS, MAX, MIN, SQRT * .. * .. Exe-
cutable 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 *
DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 )
DELTA( 1 ) = ONE
RETURN
END IF
IF( N.EQ.2 ) THEN
CALL SLAED5( I, D, Z, DELTA, RHO, DLAM )
RETURN
END IF * * Compute machine epsilon *
EPS = SLAMCH( '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 *
MIDPT = RHO / TWO * * If ||Z||_2 is not one, then TEMP
should be set to * RHO * ||Z||_2^2 / TWO *
DO 10 J = 1, N
DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
10 CONTINUE *
PSI = ZERO
DO 20 J = 1, N - 2
PSI = PSI + Z( J )*Z( J ) / DELTA( J )
20 CONTINUE *
C = RHOINV + PSI
W = C + Z( II )*Z( II ) / DELTA( II ) +
$ Z( N )*Z( N ) / DELTA( N ) *
IF( W.LE.ZERO ) THEN
TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) +
$ Z( N )*Z( N ) / RHO
IF( C.LE.TEMP ) THEN
TAU = RHO
ELSE
DEL = D( N ) - D( N-1 )
A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
B = Z( N )*Z( N )*DEL
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)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO *
DLTLB = MIDPT
DLTUB = RHO
ELSE
DEL = D( N ) - D( N-1 )
A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
B = Z( N )*Z( N )*DEL
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) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2 *
DLTLB = ZERO
DLTUB = MIDPT
END IF *
DO 30 J = 1, N
DELTA( J ) = ( D( J )-D( I ) ) - TAU
30 CONTINUE * * Evaluate PSI and the derivative DPSI *
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 40 J = 1, II
TEMP = Z( J ) / DELTA( J )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
40 CONTINUE
ERRETM = ABS( ERRETM ) * * Evaluate PHI and the deriva-
tive DPHI *
TEMP = Z( 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 * * Test for convergence *
IF( ABS( W ).LE.EPS*ERRETM ) THEN
DLAM = D( I ) + TAU
GO TO 250
END IF *
IF( W.LE.ZERO ) THEN
DLTLB = MAX( DLTLB, TAU )
ELSE
DLTUB = MIN( DLTUB, TAU )
END IF * * Calculate the new step *
NITER = NITER + 1
C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
A = ( DELTA( N-1 )+DELTA( N ) )*W -
$ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
B = DELTA( N-1 )*DELTA( N )*W
IF( C.LT.ZERO )
$ C = ABS( C )
IF( C.EQ.ZERO ) THEN * ETA = B/A * ETA =
RHO - TAU
ETA = DLTUB - TAU
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 nega-
tive, 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 = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
END IF
END IF
DO 50 J = 1, N
DELTA( J ) = DELTA( J ) - ETA
50 CONTINUE *
TAU = TAU + ETA * * Evaluate PSI and the derivative
DPSI *
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 60 J = 1, II
TEMP = Z( J ) / DELTA( J )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
60 CONTINUE
ERRETM = ABS( ERRETM ) * * Evaluate PHI and the deriva-
tive DPHI *
TEMP = Z( 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 val-
ues of the array DELTA *
ITER = NITER + 1 *
DO 90 NITER = ITER, MAXIT * * Test for convergence *
IF( ABS( W ).LE.EPS*ERRETM ) THEN
DLAM = D( I ) + TAU
GO TO 250
END IF *
IF( W.LE.ZERO ) THEN
DLTLB = MAX( DLTLB, TAU )
ELSE
DLTUB = MIN( DLTUB, TAU )
END IF * * Calculate the new step *
C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
A = ( DELTA( N-1 )+DELTA( N ) )*W -
$ DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
B = DELTA( N-1 )*DELTA( N )*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 = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
END IF
END IF
DO 70 J = 1, N
DELTA( J ) = DELTA( J ) - ETA
70 CONTINUE *
TAU = TAU + ETA * * Evaluate PSI and the deriva-
tive DPSI *
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 80 J = 1, II
TEMP = Z( 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 ) / 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
DLAM = D( I ) + TAU
GO TO 250 * * End for the case I = N *
ELSE * * The case for I < N *
NITER = 1
IP1 = I + 1 * * Calculate initial guess *
DEL = D( IP1 ) - D( I )
MIDPT = DEL / TWO
DO 100 J = 1, N
DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
100 CONTINUE *
PSI = ZERO
DO 110 J = 1, I - 1
PSI = PSI + Z( J )*Z( J ) / DELTA( J )
110 CONTINUE *
PHI = ZERO
DO 120 J = N, I + 2, -1
PHI = PHI + Z( J )*Z( J ) / DELTA( J )
120 CONTINUE
C = RHOINV + PSI + PHI
W = C + Z( I )*Z( I ) / DELTA( I ) +
$ Z( IP1 )*Z( IP1 ) / DELTA( IP1 ) *
IF( W.GT.ZERO ) THEN * * d(i)< the ith eigenvalue <
(d(i)+d(i+1))/2 * * We choose d(i) as origin. *
ORGATI = .TRUE.
A = C*DEL + Z( I )*Z( I ) + Z( IP1 )*Z( IP1 )
B = Z( I )*Z( I )*DEL
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
DLTLB = ZERO
DLTUB = MIDPT
ELSE * * (d(i)+d(i+1))/2 <= the ith eigenvalue <
d(i+1) * * We choose d(i+1) as origin. *
ORGATI = .FALSE.
A = C*DEL - Z( I )*Z( I ) - Z( IP1 )*Z( IP1 )
B = Z( IP1 )*Z( IP1 )*DEL
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
DLTLB = -MIDPT
DLTUB = ZERO
END IF *
IF( ORGATI ) THEN
DO 130 J = 1, N
DELTA( J ) = ( D( J )-D( I ) ) - TAU
130 CONTINUE
ELSE
DO 140 J = 1, N
DELTA( J ) = ( D( J )-D( IP1 ) ) - TAU
140 CONTINUE
END IF
IF( ORGATI ) THEN
II = I
ELSE
II = I + 1
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 ) / DELTA( J )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
150 CONTINUE
ERRETM = ABS( ERRETM ) * * Evaluate PHI and the deriva-
tive DPHI *
DPHI = ZERO
PHI = ZERO
DO 160 J = N, IIP1, -1
TEMP = Z( 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 secu-
lar 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 ) / 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
IF( ORGATI ) THEN
DLAM = D( I ) + TAU
ELSE
DLAM = D( IP1 ) + TAU
END IF
GO TO 250
END IF *
IF( W.LE.ZERO ) THEN
DLTLB = MAX( DLTLB, TAU )
ELSE
DLTUB = MIN( DLTUB, TAU )
END IF * * Calculate the new step *
NITER = NITER + 1
IF( .NOT.SWTCH3 ) THEN
IF( ORGATI ) THEN
C = W - DELTA( IP1 )*DW - ( D( I )-D( IP1 ) )*
$ ( Z( I ) / DELTA( I ) )**2
ELSE
C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
$ ( Z( IP1 ) / DELTA( IP1 ) )**2
END IF
A = ( DELTA( I )+DELTA( IP1 ) )*W -
$ DELTA( I )*DELTA( IP1 )*DW
B = DELTA( I )*DELTA( IP1 )*W
IF( C.EQ.ZERO ) THEN
IF( A.EQ.ZERO ) THEN
IF( ORGATI ) THEN
A = Z( I )*Z( I ) + DELTA( IP1 )*DELTA( IP1 )*
$ ( DPSI+DPHI )
ELSE
A = Z( IP1 )*Z( IP1 ) + DELTA( I )*DELTA( I )*
$ ( 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 *
TEMP = RHOINV + PSI + PHI
IF( ORGATI ) THEN
TEMP1 = Z( IIM1 ) / DELTA( IIM1 )
TEMP1 = TEMP1*TEMP1
C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) -
$ ( D( IIM1 )-D( IIP1 ) )*TEMP1
ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*
$ ( ( DPSI-TEMP1 )+DPHI )
ELSE
TEMP1 = Z( IIP1 ) / DELTA( IIP1 )
TEMP1 = TEMP1*TEMP1
C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) -
$ ( D( IIP1 )-D( IIM1 ) )*TEMP1
ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*
$ ( DPSI+( DPHI-TEMP1 ) )
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
END IF
ZZ( 2 ) = Z( II )*Z( II )
CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
$ INFO )
IF( INFO.NE.0 )
$ GO TO 250
END IF * * Note, eta should be positive if w is nega-
tive, 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
TEMP = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
END IF
END IF *
PREW = W *
DO 180 J = 1, N
DELTA( J ) = DELTA( J ) - ETA
180 CONTINUE * * Evaluate PSI and the derivative DPSI *
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 190 J = 1, IIM1
TEMP = Z( J ) / DELTA( J )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
190 CONTINUE
ERRETM = ABS( ERRETM ) * * Evaluate PHI and the deriva-
tive DPHI *
DPHI = ZERO
PHI = ZERO
DO 200 J = N, IIP1, -1
TEMP = Z( J ) / DELTA( J )
PHI = PHI + Z( J )*TEMP
DPHI = DPHI + TEMP*TEMP
ERRETM = ERRETM + PHI
200 CONTINUE *
TEMP = Z( 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+ETA )*DW *
SWTCH = .FALSE.
IF( ORGATI ) THEN
IF( -W.GT.ABS( PREW ) / TEN )
$ SWTCH = .TRUE.
ELSE
IF( W.GT.ABS( PREW ) / TEN )
$ SWTCH = .TRUE.
END IF *
TAU = TAU + ETA * * Main loop to update the values of
the array DELTA *
ITER = NITER + 1 *
DO 240 NITER = ITER, MAXIT * * Test for convergence
*
IF( ABS( W ).LE.EPS*ERRETM ) THEN
IF( ORGATI ) THEN
DLAM = D( I ) + TAU
ELSE
DLAM = D( IP1 ) + TAU
END IF
GO TO 250
END IF *
IF( W.LE.ZERO ) THEN
DLTLB = MAX( DLTLB, TAU )
ELSE
DLTUB = MIN( DLTUB, TAU )
END IF * * Calculate the new step *
IF( .NOT.SWTCH3 ) THEN
IF( .NOT.SWTCH ) THEN
IF( ORGATI ) THEN
C = W - DELTA( IP1 )*DW -
$ ( D( I )-D( IP1 ) )*( Z( I ) / DELTA( I ) )**2
ELSE
C = W - DELTA( I )*DW - ( D( IP1 )-D( I ) )*
$ ( Z( IP1 ) / DELTA( IP1 ) )**2
END IF
ELSE
TEMP = Z( II ) / DELTA( II )
IF( ORGATI ) THEN
DPSI = DPSI + TEMP*TEMP
ELSE
DPHI = DPHI + TEMP*TEMP
END IF
C = W - DELTA( I )*DPSI - DELTA( IP1 )*DPHI
END IF
A = ( DELTA( I )+DELTA( IP1 ) )*W -
$ DELTA( I )*DELTA( IP1 )*DW
B = DELTA( I )*DELTA( IP1 )*W
IF( C.EQ.ZERO ) THEN
IF( A.EQ.ZERO ) THEN
IF( .NOT.SWTCH ) THEN
IF( ORGATI ) THEN
A = Z( I )*Z( I ) + DELTA( IP1 )*
$ DELTA( IP1 )*( DPSI+DPHI )
ELSE
A = Z( IP1 )*Z( IP1 ) +
$ DELTA( I )*DELTA( I )*( DPSI+DPHI )
END IF
ELSE
A = DELTA( I )*DELTA( I )*DPSI +
$ DELTA( IP1 )*DELTA( IP1 )*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 rele-
vant poles *
TEMP = RHOINV + PSI + PHI
IF( SWTCH ) THEN
C = TEMP - DELTA( IIM1 )*DPSI - DELTA( IIP1 )*DPHI
ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*DPSI
ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*DPHI
ELSE
IF( ORGATI ) THEN
TEMP1 = Z( IIM1 ) / DELTA( IIM1 )
TEMP1 = TEMP1*TEMP1
C = TEMP - DELTA( IIP1 )*( DPSI+DPHI ) -
$ ( D( IIM1 )-D( IIP1 ) )*TEMP1
ZZ( 1 ) = Z( IIM1 )*Z( IIM1 )
ZZ( 3 ) = DELTA( IIP1 )*DELTA( IIP1 )*
$ ( ( DPSI-TEMP1 )+DPHI )
ELSE
TEMP1 = Z( IIP1 ) / DELTA( IIP1 )
TEMP1 = TEMP1*TEMP1
C = TEMP - DELTA( IIM1 )*( DPSI+DPHI ) -
$ ( D( IIP1 )-D( IIM1 ) )*TEMP1
ZZ( 1 ) = DELTA( IIM1 )*DELTA( IIM1 )*
$ ( DPSI+( DPHI-TEMP1 ) )
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
END IF
END IF
CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W,
ETA,
$ INFO )
IF( INFO.NE.0 )
$ GO TO 250
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
TEMP = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
END IF
END IF *
DO 210 J = 1, N
DELTA( J ) = DELTA( J ) - ETA
210 CONTINUE *
TAU = TAU + ETA
PREW = W * * Evaluate PSI and the derivative DPSI
*
DPSI = ZERO
PSI = ZERO
ERRETM = ZERO
DO 220 J = 1, IIM1
TEMP = Z( J ) / DELTA( J )
PSI = PSI + Z( J )*TEMP
DPSI = DPSI + TEMP*TEMP
ERRETM = ERRETM + PSI
220 CONTINUE
ERRETM = ABS( ERRETM ) * * Evaluate PHI and the
derivative DPHI *
DPHI = ZERO
PHI = ZERO
DO 230 J = N, IIP1, -1
TEMP = Z( J ) / DELTA( J )
PHI = PHI + Z( J )*TEMP
DPHI = DPHI + TEMP*TEMP
ERRETM = ERRETM + PHI
230 CONTINUE *
TEMP = Z( 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 *
240 CONTINUE * * Return with INFO = 1, NITER = MAXIT and
not converged *
INFO = 1
IF( ORGATI ) THEN
DLAM = D( I ) + TAU
ELSE
DLAM = D( IP1 ) + TAU
END IF *
END IF *
250 CONTINUE *
RETURN * * End of SLAED4 *
END
ARGUMENTS
7 Nov 2015 slaed4(3P)