| C
C Include file for CPZERO, Jenkins-Traub complex polynomial rootfinder
C
INTEGER*4 MAXDEG
PARAMETER (MAXDEG = 99)
INTEGER*4 NN
REAL*8 SR, SI, TR, TI, PVR, PVI, ARE, MRE, ETA, INFIN,
1 PR(MAXDEG), PI(MAXDEG), HR(MAXDEG), HI(MAXDEG),
2 QPR(MAXDEG), QPI(MAXDEG), QHR(MAXDEG), QHI(MAXDEG),
3 SHR(MAXDEG), SHI(MAXDEG)
COMMON /CPZERO/ PR, PI, HR, HI, QPR, QPI, QHR, QHI, SHR, SHI,
1 SR, SI, TR, TI, PVR, PVI, ARE, MRE, ETA, INFIN, NN
SUBROUTINE CPZERO(OPR, OPI, NDEG, ZEROR, ZEROI, FAIL)
C
C ACM Algorithm 419
C
C Find all zeros of a complex polynomial P(Z) by 3 stage variable
C shift Jenkins-Traub iteration. Zeroes are found roughly in
C order of increasing modulus, and the polynomial is deflated after
C each zero is found.
C
C OPR(K),
C OPI(K) = Real and imaginary components of coefficients in
C order of decreasing powers of Z
C
C NDEG = Degree of polynomial
C
C ZEROR(K),
C ZEROI(K) = Real and imaginary components of zeroes, roughly in
C increasing order of modulus
C
C FAIL = TRUE if major and minor passes failed to converge
C
IMPLICIT NONE
INCLUDE 'CPZERO.INC'
INTEGER*4 NDEG, CNT1, CNT2, I
REAL*8 OPR(1), OPI(1), ZEROR(1), ZEROI(1)
REAL*8 XX, YY, COSR, SINR, SMALNO, BASE, TT, ZR, ZI, BND,
1 CMOD, SCALE, CAUCHY, DSQRT
LOGICAL FAIL, CONV
C
C Initialize machine specific constants
C
CALL MCON(ETA, INFIN, SMALNO, BASE)
ARE = ETA
MRE = 2.0*DSQRT(2.0D0)*ETA
XX = 0.70710678
YY = -XX
COSR = -0.069756474
SINR = 0.99756405
FAIL = .FALSE.
NN = NDEG+1
C
C Algorithm fails if leading coefficient is zero
C
IF (OPR(1).NE.0.0.OR.OPI(1).NE.0.0) GOTO 10
FAIL = .TRUE.
RETURN
C
C Remove zeros at the origin if any
C
10 IF (OPR(NN).NE.0.0.OR.OPI(NN).NE.0.0) GOTO 20
ZEROR(NDEG-NN+2) = 0.0
ZEROI(NDEG-NN+2) = 0.0
NN = NN-1
GOTO 10
C
C Make copy of the coefficients
C
20 DO 30 I = 1,NN
PR(I) = OPR(I)
PI(I) = OPI(I)
SHR(I) = CMOD(PR(I),PI(I))
30 CONTINUE
C
C Scale the polynomial
C
BND = SCALE(NN, SHR, ETA, INFIN, SMALNO, BASE)
IF (BND.EQ.1.0) GOTO 40
DO 35 I = 1,NN
PR(I) = BND*PR(I)
PI(I) = BND*PI(I)
35 CONTINUE
C
C Start the algorithm for one zero
C
40 IF (NN.GT.2) GOTO 50
IF (NN.EQ.1) RETURN
C
C Calculate final zero and return
C
CALL CDIVID(-PR(2), -PI(2), PR(1), PI(1),
1 ZEROR(NDEG), ZEROI(NDEG))
RETURN
C
C Calculate BND, a lower bound on the modulus of the zeroes using
C Cauchy's estimate
C
50 DO 60 I = 1,NN
SHR(I) = CMOD(PR(I), PI(I))
60 CONTINUE
BND = CAUCHY(NN, SHR, SHI)
C
C Outer loop to control two major passes with different sequences of shifts
C
DO 100 CNT1 = 1,2
C
C First stage calculation with 5 no shift iterations
C
CALL NOSHFT(5)
C
C Inner loop to select a shift
C
DO 90 CNT2 = 1,9
C
C Shift is chosen with modulus BND and amplitude rotated by 94 degrees
C from previous shift
C
TT = COSR*XX-SINR*YY
YY = SINR*XX+COSR*YY
XX = TT
SR = BND*XX
SI = BND*YY
C
C Second stage calculation with fixed shift
C
CALL FXSHFT(10*CNT2, ZR, ZI, CONV)
IF (.NOT.CONV) GOTO 80
C
C The second stage jumps directly to the third stage iteration.
C If success, the zero is stored and the polynomial is deflated.
C
ZEROR(NDEG-NN+2) = ZR
ZEROI(NDEG-NN+2) = ZI
NN = NN-1
DO 70 I = 1,NN
PR(I) = QPR(I)
PI(I) = QPI(I)
70 CONTINUE
GOTO 40
80 CONTINUE
C
C If iteration was unsuccessful, another shift is tried
C
90 CONTINUE
C
C If 9 shifts fail, the outer loop is repeated with another sequence
C of shifts
C
100 CONTINUE
C
C Failed on two major passes; return empty handed
C
FAIL = .TRUE.
RETURN
END
SUBROUTINE NOSHFT(L1)
C
C Compute derivative of polynomial as initial H polynomial and
C compute L1 no-shift H polynomials
C
IMPLICIT NONE
INCLUDE 'CPZERO.INC'
INTEGER*4 L1, I, J, JJ, N
REAL*8 XNI, T1, T2, CMOD, XN
N = NN-1
XN = N
DO 10 I = 1,N
XNI = NN-I
HR(I) = XNI*PR(I)/XN
HI(I) = XNI*PI(I)/XN
10 CONTINUE
DO 50 JJ = 1,L1
IF (CMOD(HR(N), HI(N)).LE.ETA*10.0D0*CMOD(PR(N), PI(N)))
1 GOTO 30
CALL CDIVID(-PR(NN), -PI(NN), HR(N), HI(N), TR, TI)
DO 20 I = 1,N-1
J = NN-I
T1 = HR(J-1)
T2 = HI(J-1)
HR(J) = TR*T1-TI*T2+PR(J)
HI(J) = TR*T2+TI*T2+PI(J)
20 CONTINUE
HR(1) = PR(1)
HI(1) = PI(1)
GOTO 50
C
C If the constant term is essentially zero, shift H coefficients
C
30 DO 40 I = 1,N-1
J = NN-I
HR(J) = HR(J-1)
HI(J) = HI(J-1)
40 CONTINUE
HR(1) = 0.0
HI(1) = 0.0
50 CONTINUE
RETURN
END
SUBROUTINE FXSHFT(L2, ZR, ZI, CONV)
C
C Compute L2 fixed-shift H polynomials and test for convergence
C
C Initiates a variable shift iteration and return with the approximate
C zero if successful
C
C L2 = Limit of fixed shift steps
C ZR, ZI = Approximate zero if CONV is true
C CONV = Flag indicating convergence of stage 3 iteration
C
IMPLICIT NONE
INCLUDE 'CPZERO.INC'
REAL*8 ZR, ZI, OTR, OTI, SVSR, SVSI, CMOD
LOGICAL CONV, TEST, PASSED, BOOL
INTEGER*4 L2, N, I, J
N = NN-1
C
C Evaluate P at S
C
CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
TEST = .TRUE.
PASSED = .FALSE.
C
C Calculate first T = -P(S)/H(S)
C
CALL CALCT(BOOL)
C
C Main loop for one second stage step
C
DO 50 J = 1,L2
OTR = TR
OTI = TI
C
C Compute next H polynomial and new T
C
CALL NEXTH(BOOL)
CALL CALCT(BOOL)
ZR = SR+TR
ZI = SI+TI
C
C Tests for convergence unless stage 3 has failed once or this is the
C last H polynomial
C
IF (BOOL.OR..NOT.TEST.OR.J.EQ.L2) GOTO 50
IF (CMOD(TR-OTR, TI-OTI).GE.0.5*CMOD(ZR, ZI)) GOTO 40
IF (.NOT.PASSED) GOTO 30
C
C The weak convergence test has been passed twice; start the third
C stage iteration after saving the current H polynomial and shift
C
DO 10 I = 1,N
SHR(I) = HR(I)
SHI(I) = HI(I)
10 CONTINUE
SVSR = SR
SVSI = SI
CALL VRSHFT(10, ZR, ZI, CONV)
IF (CONV) RETURN
C
C The iteration failed to converge; turn off testing and restore
C H, S, PV, and T
C
TEST = .FALSE.
DO 20 I=1,N
HR(I) = SHR(I)
HI(I) = SHI(I)
20 CONTINUE
SR = SVSR
SI = SVSI
CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
CALL CALCT(BOOL)
GOTO 50
30 PASSED = .TRUE.
GOTO 50
40 PASSED = .FALSE.
50 CONTINUE
C
C Attempt an iteration with final H polynomial from second stage
C
CALL VRSHFT(10, ZR, ZI, CONV)
RETURN
END
SUBROUTINE VRSHFT(L3, ZR, ZI, CONV)
C
C Carry out third stage iteration
C
C L3 = Limit of steps in stage 3
C ZR, ZI = On entry contains initial iterate
C if iteration converges, contains final iterate
C CONV = .TRUE. if success
C
IMPLICIT NONE
INCLUDE 'CPZERO.INC'
INTEGER*4 L3, I, J
LOGICAL CONV, B, BOOL
REAL*8 ZR, ZI, MP, MS, OMP, RELSTP, R1, R2, CMOD, DSQRT,
1 ERREV, TP
CONV = .FALSE.
B = .FALSE.
SR = ZR
SI = ZI
C
C Main loop for stage 3
C
DO 60 I = 1,L3
C
C Evaluate P at S and test for convergence
C
CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
MP = CMOD(PVR, PVI)
MS = CMOD(SR, SI)
IF (MP.GT.20.0*ERREV(NN, QPR, QPI, MS, MP, ARE, MRE))
1 GOTO 10
C
C Success - polynomial value is smaller than a bound on the error
C
CONV = .TRUE.
ZR = SR
ZI = SI
RETURN
10 IF (I.EQ.1) GOTO 40
IF (B.OR.MP.LT.OMP.OR.RELSTP.GE.0.5) GOTO 30
C
C Iteration has stalled, probably due to a cluster of zeros
C Do 5 fixed shift iterations into the cluster to try and force one
C zero to dominate
C
TP = RELSTP
B = .TRUE.
IF (RELSTP.LT.ETA) TP = ETA
R1 = DSQRT(TP)
R2 = SR*(1.0+R1)-SI*R1
SI = SR*R1+SI*(1.0+R1)
SR = R2
CALL POLYEV(NN, SR, SI, PR, PI, QPR, QPI, PVR, PVI)
DO 20 J = 1,5
CALL CALCT(BOOL)
CALL NEXTH(BOOL)
20 CONTINUE
OMP = INFIN
GOTO 50
C
C Return if polynomial value increases significantly
C
30 IF (MP*0.1.GT.OMP) RETURN
40 OMP = MP
C
C Calculate next iterate
C
50 CALL CALCT(BOOL)
CALL NEXTH(BOOL)
CALL CALCT(BOOL)
IF (BOOL) GOTO 60
RELSTP = CMOD(TR, TI)/CMOD(SR, SI)
SR = SR+TR
SI = SI+TI
60 CONTINUE
RETURN
END
SUBROUTINE CALCT(BOOL)
C
C Calculate T = -P(S)/H(S)
C
C BOOL = True if H(S) is essentially zero
C
IMPLICIT NONE
INCLUDE 'CPZERO.INC'
INTEGER*4 N
REAL*8 HVR, HVI, CMOD
LOGICAL BOOL
N = NN-1
C
C Evaluate H(S)
C
CALL POLYEV(N, SR, SI, HR, HI, QHR, QHI, HVR, HVI)
BOOL = CMOD(HVR, HVI).LE.ARE*10.0*CMOD(HR(N), HI(N))
IF (BOOL) GOTO 10
CALL CDIVID(-PVR, -PVI, HVR, HVI, TR, TI)
RETURN
10 TR = 0.0
TI = 0.0
RETURN
END
SUBROUTINE NEXTH(BOOL)
C
C Calculate the next shifted H polynomial
C
IMPLICIT NONE
INCLUDE 'CPZERO.INC'
REAL*8 T1, T2
INTEGER*4 I, J
LOGICAL BOOL
IF (BOOL) GOTO 20
DO 10 J = 2,NN-1
T1 = QHR(J-1)
T2 = QHI(J-1)
HR(J) = TR*T1-TI*T2+QPR(J)
HI(J) = TR*T2+TI*T1+QPI(J)
10 CONTINUE
HR(1) = QPR(1)
HI(1) = QPI(1)
RETURN
C
C If H(S) is zero replace H with QH
C
20 DO 30 J = 2,NN-1
HR(J) = QHR(J-1)
HI(J) = QHI(J-1)
30 CONTINUE
HR(1) = 0.0
HI(1) = 0.0
RETURN
END
SUBROUTINE POLYEV(NN, SR, SI, PR, PI, QR, QI, PVR, PVI)
C
C Evaluate a polynomial via Horner recurrance placing partial
C sums in Q and computed value in PV
C
IMPLICIT NONE
INTEGER*4 I, NN
REAL*8 PR(NN), PI(NN), QR(NN), QI(NN), SR, SI, PVR, PVI, PVT
QR(1) = PR(1)
QI(1) = PI(1)
PVR = QR(1)
PVI = QI(1)
DO 10 I = 2,NN
PVT = PVR*SR-PVI*SI+PR(I)
PVI = PVR*SI+PVI*SR+PI(I)
PVR = PVT
QR(I) = PVR
QI(I) = PVI
10 CONTINUE
RETURN
END
REAL*8 FUNCTION ERREV(NN, QR, QI, MS, MP, ARE, MRE)
C
C Bounds the error in evaluating the polynomial by Horner recurrance
C
C QR, QI = Partial sums
C MS = Modulus of S
C MP = Modulus of P(S)
C ARE, MRE Relative error bounds on complex addition and multiplication
C
IMPLICIT NONE
INTEGER*4 I, NN
REAL*8 QR(NN), QI(NN), MS, MP, ARE, MRE, E, CMOD
E = CMOD(QR(1), QI(1))*MRE/(ARE+MRE)
DO 10 I = 1,NN
E = E*MS+CMOD(QR(I), QI(I))
10 CONTINUE
ERREV = E*(ARE+MRE)-MP*MRE
RETURN
END
REAL*8 FUNCTION CAUCHY(NN, PT, Q)
C
C Compute a lower bound on the modulus of the zeros of a polynomial
C using Cauchy's estimate.
C
C PT(K) = Modulus of K'th coefficient
C
IMPLICIT NONE
INTEGER*4 N, NN, I
REAL*8 Q(NN), PT(NN), X, XM, F, DX, DF, DABS, DEXP, DLOG, XN
N = NN-1
XN = N
PT(NN) = -PT(NN)
C
C Compute upper estimate of bound
C
X = DEXP((DLOG(-PT(NN))-DLOG(PT(1)))/XN)
IF (PT(N).EQ.0.0) GOTO 20
C
C If Newton step at origin is better, use it
C
XM = -PT(NN)/PT(N)
IF (XN.LT.X) X = XM
C
C Chop the interval [0,X] until P <= 0
C
20 XM = X*0.1
F = PT(1)
DO 30 I = 2,NN
F = F*XM+PT(I)
30 CONTINUE
IF (F.LE.0.0) GOTO 40
X = XM
GOTO 20
40 DX = X
C
C Do Newton iteration until X converges to two decimal places
C
50 IF (DABS(DX/X).LE.0.005) GOTO 70
Q(1) = PT(1)
DO 60 I = 2,NN
Q(I) = Q(I-1)*X+PT(I)
60 CONTINUE
F = Q(NN)
DF = Q(1)
DO 65 I=2,N
DF = DF*X+Q(I)
65 CONTINUE
DX = F/DF
X = X-DX
GOTO 50
70 CAUCHY = X
RETURN
END
REAL*8 FUNCTION SCALE(NN, PT, ETA, INFIN, SMALNO, BASE)
C
C Return a scale factor to multiply the coefficients of the polynomial
C
C The scaling is done to avoid overflow and to avoid undetected underflow
C interfering with the convergence criterion. The factor is a power of
C the base.
C
C PT(K) = Modulus of K'th coefficient
C ETA, INFIN,
C SMALNO, BASE = Constants describing machine arithmetic
C
IMPLICIT NONE
INTEGER*4 NN, I, L
REAL*8 PT(NN), ETA, INFIN, SMALNO, BASE, HI, LO,
1 MAX, MIN, X, SC, DSQRT, DLOG
C
C Find largest and smallest moduli of coefficients
C
HI = DSQRT(INFIN)
LO = SMALNO/ETA
MAX = 0.0
MIN = INFIN
DO 10 I = 1,NN
X = PT(I)
IF (X.GT.MAX) MAX = X
IF (X.NE.0.0.AND.X.LT.MIN) MIN = X
10 CONTINUE
C
C Scale only if there are very large or small components
C
SCALE = 1.0
IF (MIN.GE.LO.AND.MAX.LE.HI) RETURN
X = LO/MIN
IF (X.GT.1.0) GOTO 20
SC = 1.0/(DSQRT(MAX)*DSQRT(MIN))
GOTO 30
20 SC = X
IF (INFIN/SC.GT.MAX) SC = 1.0
30 L = DLOG(SC)/DLOG(BASE)+0.5
SCALE = BASE**L
RETURN
END
SUBROUTINE CDIVID(AR, AI, BR, BI, CR, CI)
C
C Complex division C = A/B avoiding overflow
C
IMPLICIT NONE
REAL*8 AR, AI, BR, BI, CR, CI, T, R, D, INFIN, DABS
IF (BR.NE.0.0.OR.BI.NE.0.0) goto 10
C
C Division by zero, return infinity
C
CALL MCON(T, INFIN, T, T)
CR = INFIN
CI = INFIN
RETURN
10 IF (DABS(BR).GE.DABS(BI)) GOTO 20
R = BR/BI
D = BI+R*BR
CR = (AR*R+AI)/D
CI = (AI*R-AR)/D
RETURN
20 R = BI/BR
D = BR+R*BI
CR = (AR+AI*R)/D
CI = (AI-AR*R)/D
RETURN
END
REAL*8 FUNCTION CMOD(R, I)
C
C Modulus of complex floating number avoiding overflow
C
IMPLICIT NONE
REAL*8 R, I, AR, AI, DABS, DSQRT
AR = DABS(R)
AI = DABS(I)
IF (AR.GE.AI) GOTO 10
CMOD = AI*DSQRT(1.0+(AR/AI)**2)
RETURN
10 IF (AR.LE.AI) GOTO 20
CMOD = AR*DSQRT(1.0+(AI/AR)**2)
RETURN
20 CMOD = AR*DSQRT(2.0D0)
RETURN
END
SUBROUTINE MCON(ETA, LARGENO, SMALLNO, BASE)
C
C
C MCON provides machine constants used in various parts of the program
C
C ETA = Maximum relative representation error which can
C be represented as the smallest positive floating point
C number such that 1.0 + ETA > 1.0
C LARGENO = Largest positive floating point number
C SMALLNO = Smallest positive floating point number
C BASE = Base of floating point system used
C
C Let T be the number of base digits in each floating point number
C Then ETA is either 0.5*B**(1-T) or B**(1-T) depending on whether
C rounding or truncation is used.
C
C Let M be the largest exponent and N be the smallest exponent.
C Then LARGENO = (1-BASE**(-T))*BASE**M
C and SMALLNO = BASE**N
C
C The values below correspond to VAX D floating numbers
C
IMPLICIT NONE
REAL*8 ETA, LARGENO, SMALLNO, BASE
REAL*8 XL, XS
INTEGER M, N, T
C
C Hexadecimal representations of large and small numbers:
C
C Large FFFFFF7F FFFFFFFF
C Small 00000080 00000000
C
BASE = 2.0
T = 56
M = 127
N = -128
ETA = BASE**-T
LARGENO = BASE**(126)*(1.0-BASE**-T)*BASE
SMALLNO = BASE**(-126)/BASE/BASE
RETURN
END
|