C DPRQ 10
C ..................................................................DPRQ 20
C DPRQ 30
C SUBROUTINE DPRQD DPRQ 40
C DPRQ 50
C PURPOSE DPRQ 60
C CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN POLYNOMIAL DPRQ 70
C WITH REAL COEFFICIENTS. DPRQ 80
C DPRQ 90
C USAGE DPRQ 100
C CALL DPRQD(C,IC,Q,E,POL,IR,IER) DPRQ 110
C DPRQ 120
C DESCRIPTION OF PARAMETERS DPRQ 130
C C - COEFFICIENT VECTOR OF GIVEN POLYNOMIAL DPRQ 140
C COEFFICIENTS ARE ORDERED FROM LOW TO HIGH DPRQ 150
C THE GIVEN COEFFICIENT VECTOR GETS DIVIDED BY THE DPRQ 160
C LAST NONZERO TERM DPRQ 170
C DOUBLE PRECISION ARRAY DPRQ 180
C IC - DIMENSION OF VECTOR C DPRQ 190
C Q - WORKING STORAGE OF DIMENSION IC DPRQ 200
C ON RETURN Q CONTAINS REAL PARTS OF ROOTS DPRQ 210
C DOUBLE PRECISION ARRAY DPRQ 220
C E - WORKING STORAGE OF DIMENSION IC DPRQ 230
C ON RETURN E CONTAINS COMPLEX PARTS OF ROOTS DPRQ 240
C DOUBLE PRECISION ARRAY DPRQ 250
C POL - WORKING STORAGE OF DIMENSION IC DPRQ 260
C ON RETURN POL CONTAINS THE COEFFICIENTS OF THE DPRQ 270
C POLYNOMIAL WITH CALCULATED ROOTS DPRQ 280
C THIS RESULTING COEFFICIENT VECTOR HAS DIMENSION IR+1DPRQ 290
C COEFFICIENTS ARE ORDERED FROM LOW TO HIGH DPRQ 300
C DOUBLE PRECISION ARRAY DPRQ 310
C IR - NUMBER OF CALCULATED ROOTS DPRQ 320
C NORMALLY IR IS EQUAL TO DIMENSION IC MINUS ONE DPRQ 330
C IER - RESULTING ERROR PARAMETER. SEE REMARKS DPRQ 340
C DPRQ 350
C REMARKS DPRQ 360
C THE REAL PART OF THE ROOTS IS STORED IN Q(1) UP TO Q(IR) DPRQ 370
C CORRESPONDING COMPLEX PARTS ARE STORED IN E(1) UP TO E(IR). DPRQ 380
C IER = 0 MEANS NO ERRORS DPRQ 390
C IER = 1 MEANS NO CONVERGENCE WITH FEASIBLE TOLERANCE DPRQ 400
C IER = 2 MEANS POLYNOMIAL IS DEGENERATE (CONSTANT OR ZERO) DPRQ 410
C IER = 3 MEANS SUBROUTINE WAS ABANDONED DUE TO ZERO DIVISOR DPRQ 420
C IER = 4 MEANS THERE EXISTS NO S-FRACTION DPRQ 430
C IER =-1 MEANS CALCULATED COEFFICIENT VECTOR REVEALS POOR DPRQ 440
C ACCURACY OF THE CALCULATED ROOTS. DPRQ 450
C THE CALCULATED COEFFICIENT VECTOR HAS LESS THAN DPRQ 460
C 6 CORRECT DIGITS. DPRQ 470
C THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED DPRQ 480
C COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE BEEN DPRQ 490
C CALCULATED. DPRQ 500
C THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT VECTOR IS DPRQ 510
C RECORDED IN Q(IR+1). DPRQ 520
C DPRQ 530
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED DPRQ 540
C NONE DPRQ 550
C DPRQ 560
C METHOD DPRQ 570
C THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF DPRQ 580
C THE QUOTIENT-DIFFERENCE ALGORITHM WITH DISPLACEMENT. DPRQ 590
C REFERENCE DPRQ 600
C H.RUTISHAUSER, DER QUOTIENTEN-DIFFERENZEN-ALGORITHMUS, DPRQ 610
C BIRKHAEUSER, BASEL/STUTTGART, 1957. DPRQ 620
C DPRQ 630
C ..................................................................DPRQ 640
C DPRQ 650
SUBROUTINE DPRQD(C,IC,Q,E,POL,IR,IER) DPRQ 660
C DPRQ 670
C DIMENSIONED DUMMY VARIABLES DPRQ 680
DIMENSION E(1),Q(1),C(1),POL(1) DPRQ 690
DOUBLE PRECISION Q,E,O,P,T,EXPT,ESAV,U,V,W,C,POL,EPS DPRQ 700
C DPRQ 710
C NORMALIZATION OF GIVEN POLYNOMIAL DPRQ 720
C TEST OF DIMENSION DPRQ 730
C IR CONTAINS INDEX OF HIGHEST COEFFICIENT DPRQ 740
IR=IC DPRQ 750
IER=0 DPRQ 760
EPS=1.D-16 DPRQ 770
TOL=1.E-6 DPRQ 780
LIMIT=10*IC DPRQ 790
KOUNT=0 DPRQ 800
1 IF(IR-1)79,79,2 DPRQ 810
C DPRQ 820
C DROP TRAILING ZERO COEFFICIENTS DPRQ 830
2 IF(C(IR))4,3,4 DPRQ 840
3 IR=IR-1 DPRQ 850
GOTO 1 DPRQ 860
C DPRQ 870
C REARRANGEMENT OF GIVEN POLYNOMIAL DPRQ 880
C EXTRACTION OF ZERO ROOTS DPRQ 890
4 O=1.0D0/C(IR) DPRQ 900
IEND=IR-1 DPRQ 910
ISTA=1 DPRQ 920
NSAV=IR+1 DPRQ 930
JBEG=1 DPRQ 940
C DPRQ 950
C Q(J)=1. DPRQ 960
C Q(J+I)=C(IR-I)/C(IR) DPRQ 970
C Q(IR)=C(J)/C(IR) DPRQ 980
C WHERE J IS THE INDEX OF THE LOWEST NONZERO COEFFICIENT DPRQ 990
DO 9 I=1,IR DPRQ1000
J=NSAV-I DPRQ1010
IF(C(I))7,5,7 DPRQ1020
5 GOTO(6,8),JBEG DPRQ1030
6 NSAV=NSAV+1 DPRQ1040
Q(ISTA)=0.D0 DPRQ1050
E(ISTA)=0.D0 DPRQ1060
ISTA=ISTA+1 DPRQ1070
GOTO 9 DPRQ1080
7 JBEG=2 DPRQ1090
8 Q(J)=C(I)*O DPRQ1100
C(I)=Q(J) DPRQ1110
9 CONTINUE DPRQ1120
C DPRQ1130
C INITIALIZATION DPRQ1140
ESAV=0.D0 DPRQ1150
Q(ISTA)=0.D0 DPRQ1160
10 NSAV=IR DPRQ1170
C DPRQ1180
C COMPUTATION OF DERIVATIVE DPRQ1190
EXPT=IR-ISTA DPRQ1200
E(ISTA)=EXPT DPRQ1210
DO 11 I=ISTA,IEND DPRQ1220
EXPT=EXPT-1.0D0 DPRQ1230
POL(I+1)=EPS*DABS(Q(I+1))+EPS DPRQ1240
11 E(I+1)=Q(I+1)*EXPT DPRQ1250
C DPRQ1260
C TEST OF REMAINING DIMENSION DPRQ1270
IF(ISTA-IEND)12,20,60 DPRQ1280
12 JEND=IEND-1 DPRQ1290
C DPRQ1300
C COMPUTATION OF S-FRACTION DPRQ1310
DO 19 I=ISTA,JEND DPRQ1320
IF(I-ISTA)13,16,13 DPRQ1330
13 IF(DABS(E(I))-POL(I+1))14,14,16 DPRQ1340
C DPRQ1350
C THE GIVEN POLYNOMIAL HAS MULTIPLE ROOTS, THE COEFFICIENTS OF DPRQ1360
C THE COMMON FACTOR ARE STORED FROM Q(NSAV) UP TO Q(IR) DPRQ1370
14 NSAV=I DPRQ1380
DO 15 K=I,JEND DPRQ1390
IF(DABS(E(K))-POL(K+1))15,15,80 DPRQ1400
15 CONTINUE DPRQ1410
GOTO 21 DPRQ1420
C DPRQ1430
C EUCLIDEAN ALGORITHM DPRQ1440
16 DO 19 K=I,IEND DPRQ1450
E(K+1)=E(K+1)/E(I) DPRQ1460
Q(K+1)=E(K+1)-Q(K+1) DPRQ1470
IF(K-I)18,17,18 DPRQ1480
C DPRQ1490
C TEST FOR SMALL DIVISOR DPRQ1500
17 IF(DABS(Q(I+1))-POL(I+1))80,80,19 DPRQ1510
18 Q(K+1)=Q(K+1)/Q(I+1) DPRQ1520
POL(K+1)=POL(K+1)/DABS(Q(I+1)) DPRQ1530
E(K)=Q(K+1)-E(K) DPRQ1540
19 CONTINUE DPRQ1550
20 Q(IR)=-Q(IR) DPRQ1560
C DPRQ1570
C THE DISPLACEMENT EXPT IS SET TO 0 AUTOMATICALLY. DPRQ1580
C E(ISTA)=0.,Q(ISTA+1),...,E(NSAV-1),Q(NSAV),E(NSAV)=0., DPRQ1590
C FORM A DIAGONAL OF THE QD-ARRAY. DPRQ1600
C INITIALIZATION OF BOUNDARY VALUES DPRQ1610
21 E(ISTA)=0.D0 DPRQ1620
NRAN=NSAV-1 DPRQ1630
22 E(NRAN+1)=0.D0 DPRQ1640
C DPRQ1650
C TEST FOR LINEAR OR CONSTANT FACTOR DPRQ1660
C NRAN-ISTA IS DEGREE-1 DPRQ1670
IF(NRAN-ISTA)24,23,31 DPRQ1680
C DPRQ1690
C LINEAR FACTOR DPRQ1700
23 Q(ISTA+1)=Q(ISTA+1)+EXPT DPRQ1710
E(ISTA+1)=0.D0 DPRQ1720
C DPRQ1730
C TEST FOR UNFACTORED COMMON DIVISOR DPRQ1740
24 E(ISTA)=ESAV DPRQ1750
IF(IR-NSAV)60,60,25 DPRQ1760
C DPRQ1770
C INITIALIZE QD-ALGORITHM FOR COMMON DIVISOR DPRQ1780
25 ISTA=NSAV DPRQ1790
ESAV=E(ISTA) DPRQ1800
GOTO 10 DPRQ1810
C DPRQ1820
C COMPUTATION OF ROOT PAIR DPRQ1830
26 P=P+EXPT DPRQ1840
C DPRQ1850
C TEST FOR REALITY DPRQ1860
IF(O)27,28,28 DPRQ1870
C DPRQ1880
C COMPLEX ROOT PAIR DPRQ1890
27 Q(NRAN)=P DPRQ1900
Q(NRAN+1)=P DPRQ1910
E(NRAN)=T DPRQ1920
E(NRAN+1)=-T DPRQ1930
GOTO 29 DPRQ1940
C DPRQ1950
C REAL ROOT PAIR DPRQ1960
28 Q(NRAN)=P-T DPRQ1970
Q(NRAN+1)=P+T DPRQ1980
E(NRAN)=0.D0 DPRQ1990
C DPRQ2000
C REDUCTION OF DEGREE BY 2 (DEFLATION) DPRQ2010
29 NRAN=NRAN-2 DPRQ2020
GOTO 22 DPRQ2030
C DPRQ2040
C COMPUTATION OF REAL ROOT DPRQ2050
30 Q(NRAN+1)=EXPT+P DPRQ2060
C DPRQ2070
C REDUCTION OF DEGREE BY 1 (DEFLATION) DPRQ2080
NRAN=NRAN-1 DPRQ2090
GOTO 22 DPRQ2100
C DPRQ2110
C START QD-ITERATION DPRQ2120
31 JBEG=ISTA+1 DPRQ2130
JEND=NRAN-1 DPRQ2140
TEPS=EPS DPRQ2150
TDELT=1.E-2 DPRQ2160
32 KOUNT=KOUNT+1 DPRQ2170
P=Q(NRAN+1) DPRQ2180
R=ABS(SNGL(E(NRAN))) DPRQ2190
C DPRQ2200
C TEST FOR CONVERGENCE DPRQ2210
IF(R-TEPS)30,30,33 DPRQ2220
33 S=ABS(SNGL(E(JEND))) DPRQ2230
C DPRQ2240
C IS THERE A REAL ROOT NEXT DPRQ2250
IF(S-R)38,38,34 DPRQ2260
C DPRQ2270
C IS DISPLACEMENT SMALL ENOUGH DPRQ2280
34 IF(R-TDELT)36,35,35 DPRQ2290
35 P=0.D0 DPRQ2300
36 O=P DPRQ2310
DO 37 J=JBEG,NRAN DPRQ2320
Q(J)=Q(J)+E(J)-E(J-1)-O DPRQ2330
C DPRQ2340
C TEST FOR SMALL DIVISOR DPRQ2350
IF(DABS(Q(J))-POL(J))81,81,37 DPRQ2360
37 E(J)=Q(J+1)*E(J)/Q(J) DPRQ2370
Q(NRAN+1)=-E(NRAN)+Q(NRAN+1)-O DPRQ2380
GOTO 54 DPRQ2390
C DPRQ2400
C CALCULATE DISPLACEMENT FOR DOUBLE ROOTS DPRQ2410
C QUADRATIC EQUATION FOR DOUBLE ROOTS DPRQ2420
C X**2-(Q(NRAN)+Q(NRAN+1)+E(NRAN))*X+Q(NRAN)*Q(NRAN+1)=0 DPRQ2430
38 P=0.5D0*(Q(NRAN)+E(NRAN)+Q(NRAN+1)) DPRQ2440
O=P*P-Q(NRAN)*Q(NRAN+1) DPRQ2450
T=DSQRT(DABS(O)) DPRQ2460
C DPRQ2470
C TEST FOR CONVERGENCE DPRQ2480
IF(S-TEPS)26,26,39 DPRQ2490
C DPRQ2500
C ARE THERE COMPLEX ROOTS DPRQ2510
39 IF(O)43,40,40 DPRQ2520
40 IF(P)42,41,41 DPRQ2530
41 T=-T DPRQ2540
42 P=P+T DPRQ2550
R=S DPRQ2560
GOTO 34 DPRQ2570
C DPRQ2580
C MODIFICATION FOR COMPLEX ROOTS DPRQ2590
C IS DISPLACEMENT SMALL ENOUGH DPRQ2600
43 IF(S-TDELT)44,35,35 DPRQ2610
C DPRQ2620
C INITIALIZATION DPRQ2630
44 O=Q(JBEG)+E(JBEG)-P DPRQ2640
C DPRQ2650
C TEST FOR SMALL DIVISOR DPRQ2660
IF(DABS(O)-POL(JBEG))81,81,45 DPRQ2670
45 T=(T/O)**2 DPRQ2680
U=E(JBEG)*Q(JBEG+1)/(O*(1.0D0+T)) DPRQ2690
V=O+U DPRQ2700
C DPRQ2710
C THREEFOLD LOOP FOR COMPLEX DISPLACEMENT DPRQ2720
KOUNT=KOUNT+2 DPRQ2730
DO 53 J=JBEG,NRAN DPRQ2740
O=Q(J+1)+E(J+1)-U-P DPRQ2750
C DPRQ2760
C TEST FOR SMALL DIVISOR DPRQ2770
IF(DABS(V)-POL(J))46,46,49 DPRQ2780
46 IF(J-NRAN)81,47,81 DPRQ2790
47 EXPT=EXPT+P DPRQ2800
IF(ABS(SNGL(E(JEND)))-TOL)48,48,81 DPRQ2810
48 P=0.5D0*(V+O-E(JEND)) DPRQ2820
O=P*P-(V-U)*(O-U*T-O*W*(1.D0+T)/Q(JEND)) DPRQ2830
T=DSQRT(DABS(O)) DPRQ2840
GOTO 26 DPRQ2850
C DPRQ2860
C TEST FOR SMALL DIVISOR DPRQ2870
49 IF(DABS(O)-POL(J+1))46,46,50 DPRQ2880
50 W=U*O/V DPRQ2890
T=T*(V/O)**2 DPRQ2900
Q(J)=V+W-E(J-1) DPRQ2910
U=0.D0 DPRQ2920
IF(J-NRAN)51,52,52 DPRQ2930
51 U=Q(J+2)*E(J+1)/(O*(1.D0+T)) DPRQ2940
52 V=O+U-W DPRQ2950
C DPRQ2960
C TEST FOR SMALL DIVISOR DPRQ2970
IF(DABS(Q(J))-POL(J))81,81,53 DPRQ2980
53 E(J)=W*V*(1.0D0+T)/Q(J) DPRQ2990
Q(NRAN+1)=V-E(NRAN) DPRQ3000
54 EXPT=EXPT+P DPRQ3010
TEPS=TEPS*1.1 DPRQ3020
TDELT=TDELT*1.1 DPRQ3030
IF(KOUNT-LIMIT)32,55,55 DPRQ3040
C DPRQ3050
C NO CONVERGENCE WITH FEASIBLE TOLERANCE DPRQ3060
C ERROR RETURN IN CASE OF UNSATISFACTORY CONVERGENCE DPRQ3070
55 IER=1 DPRQ3080
C DPRQ3090
C REARRANGE CALCULATED ROOTS DPRQ3100
56 IEND=NSAV-NRAN-1 DPRQ3110
E(ISTA)=ESAV DPRQ3120
IF(IEND)59,59,57 DPRQ3130
57 DO 58 I=1,IEND DPRQ3140
J=ISTA+I DPRQ3150
K=NRAN+1+I DPRQ3160
E(J)=E(K) DPRQ3170
58 Q(J)=Q(K) DPRQ3180
59 IR=ISTA+IEND DPRQ3190
C DPRQ3200
C NORMAL RETURN DPRQ3210
60 IR=IR-1 DPRQ3220
IF(IR)78,78,61 DPRQ3230
C DPRQ3240
C REARRANGE CALCULATED ROOTS DPRQ3250
61 DO 62 I=1,IR DPRQ3260
Q(I)=Q(I+1) DPRQ3270
62 E(I)=E(I+1) DPRQ3280
C DPRQ3290
C CALCULATE COEFFICIENT VECTOR FROM ROOTS DPRQ3300
POL(IR+1)=1.D0 DPRQ3310
IEND=IR-1 DPRQ3320
JBEG=1 DPRQ3330
DO 69 J=1,IR DPRQ3340
ISTA=IR+1-J DPRQ3350
O=0.D0 DPRQ3360
P=Q(ISTA) DPRQ3370
T=E(ISTA) DPRQ3380
IF(T)65,63,65 DPRQ3390
C DPRQ3400
C MULTIPLY WITH LINEAR FACTOR DPRQ3410
63 DO 64 I=ISTA,IR DPRQ3420
POL(I)=O-P*POL(I+1) DPRQ3430
64 O=POL(I+1) DPRQ3440
GOTO 69 DPRQ3450
65 GOTO(66,67),JBEG DPRQ3460
66 JBEG=2 DPRQ3470
POL(ISTA)=0.D0 DPRQ3480
GOTO 69 DPRQ3490
C DPRQ3500
C MULTIPLY WITH QUADRATIC FACTOR DPRQ3510
67 JBEG=1 DPRQ3520
U=P*P+T*T DPRQ3530
P=P+P DPRQ3540
DO 68 I=ISTA,IEND DPRQ3550
POL(I)=O-P*POL(I+1)+U*POL(I+2) DPRQ3560
68 O=POL(I+1) DPRQ3570
POL(IR)=O-P DPRQ3580
69 CONTINUE DPRQ3590
IF(IER)78,70,78 DPRQ3600
C DPRQ3610
C COMPARISON OF COEFFICIENT VECTORS, IE. TEST OF ACCURACY DPRQ3620
70 P=0.D0 DPRQ3630
DO 75 I=1,IR DPRQ3640
IF(C(I))72,71,72 DPRQ3650
71 O=DABS(POL(I)) DPRQ3660
GOTO 73 DPRQ3670
72 O=DABS((POL(I)-C(I))/C(I)) DPRQ3680
73 IF(P-O)74,75,75 DPRQ3690
74 P=O DPRQ3700
75 CONTINUE DPRQ3710
IF(SNGL(P)-TOL)77,76,76 DPRQ3720
76 IER=-1 DPRQ3730
77 Q(IR+1)=P DPRQ3740
E(IR+1)=0.D0 DPRQ3750
78 RETURN DPRQ3760
C DPRQ3770
C ERROR RETURNS DPRQ3780
C ERROR RETURN FOR POLYNOMIALS OF DEGREE LESS THAN 1 DPRQ3790
79 IER=2 DPRQ3800
IR=0 DPRQ3810
RETURN DPRQ3820
C DPRQ3830
C ERROR RETURN IF THERE EXISTS NO S-FRACTION DPRQ3840
80 IER=4 DPRQ3850
IR=ISTA DPRQ3860
GOTO 60 DPRQ3870
C DPRQ3880
C ERROR RETURN IN CASE OF INSTABLE QD-ALGORITHM DPRQ3890
81 IER=3 DPRQ3900
GOTO 56 DPRQ3910
END DPRQ3920