This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / c / cpolyz64.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:01:53  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10 #if defined(CERNLIB_DOUBLE)
11       SUBROUTINE WPOLYZ(C,N,MAXIT,Z0,H)
12 C
13 #include "gen/imp64.inc"
14 C
15       CHARACTER NAME*(*)
16       PARAMETER (NAME = 'WPOLYZ')
17 #endif
18 #if !defined(CERNLIB_DOUBLE)
19       SUBROUTINE CPOLYZ(C,N,MAXIT,Z0,H)
20 C
21       CHARACTER NAME*(*)
22       PARAMETER (NAME = 'CPOLYZ')
23 #endif
24 C
25 C     Based on
26 C      T. Pomentale, Homotopy methods for polynomial equations,
27 C      J. Inst. Maths Applics 13 (1974) 201-213
28
29 C     Revised version of a routine originally written by T. Pomentale
30
31 #include "gen/defc64.inc"
32      + C(0:*),Z0(*),F(100),B(100),Z,CN,CD,FZ,FF,F1,F2,V,ZZ,Z1,ZV,IU
33       LOGICAL LMD
34       DIMENSION H(*)
35       CHARACTER*80 ERRTXT
36
37       PARAMETER (EPS = 1D-14)
38 #if defined(CERNLIB_CRAY)
39       PARAMETER (IU = (0E0,1E0))
40 #endif
41 #if !defined(CERNLIB_CRAY)
42       PARAMETER (IU = (0D0,1D0))
43 #endif
44       PARAMETER (R0 = 0, R1 = 1, BIG1 = 1D7, BIG2 = BIG1/10, HF = R1/2)
45       PARAMETER (FAC1 = 9*R1/10, FAC2 = 2+100*FAC1*EPS)
46       PARAMETER (PI = 3.14159 26535 89793 24D0, PI2 = 2*PI)
47
48 #if defined(CERNLIB_QF2C)
49 #include "defdr.inc"
50 #endif
51 #if defined(CERNLIB_DOUBLE)
52       TST(Z1,ZV)=DREAL(Z1)*DREAL(ZV)+DIMAG(Z1)*DIMAG(ZV)
53 #endif
54 #if !defined(CERNLIB_DOUBLE)
55       TST(Z1,ZV)=REAL(Z1)*REAL(ZV)+AIMAG(Z1)*AIMAG(ZV)
56 #endif
57
58       IF(N .LE. 0) RETURN
59       IF(C(0) .EQ. 0) THEN
60        CALL MTLPRT(NAME,'C209.1','A(0) = 0')
61        RETURN
62       ENDIF
63       IF(N .EQ. 1) THEN
64        Z0(1)=-C(1)/C(0)
65        H(1)=0
66        RETURN
67       ENDIF
68       IF(N .EQ. 2) THEN
69        H(1)=0
70        H(2)=0
71        Z0(2)=(-C(1)+SQRT(C(1)**2-4*C(0)*C(2)))/(2*C(0))
72        Z0(1)=-Z0(2)-C(1)/C(0)
73        RETURN
74       ENDIF
75       LMD=.FALSE.
76       DO 1 K = 1,N
77     1 LMD=LMD .OR. Z0(K) .NE. 0
78       M=0
79       WN=SQRT(N+R0)
80       Z=1+IU
81       N1=N
82
83       DO 3 K = 1,N+1
84       H(K)=1
85     3 B(K)=C(K-1)
86     6 DO 10 K = 1,N-2
87       IF(B(N1+1) .EQ. 0) GO TO 18
88       IF(LMD) Z=Z0(K)
89       IF(M .GE. MAXIT) THEN
90        CALL MTLPRT(NAME,'C209.2','TOO MANY ITERATIONS')
91        RETURN
92       ENDIF
93
94       S=BIG1
95       DO 61 I = 1,N1+1
96    61 F(I)=B(I)
97       ZA=ABS(Z)
98       HH=HF*ABS(B(1))
99       DO 62 I = 2,N1+1
100       F(I)=Z*F(I-1)+F(I)
101       D=ABS(F(I))
102    62 HH=ZA*HH+D
103       DO 63 L = 1,2
104       DO 63 I = 2,N1-L+1
105    63 F(I)=Z*F(I-1)+F(I)
106       M=M+1
107       IF(D .LE. EPS*(FAC2*HH-D)) GO TO 13
108
109     7 D1=S
110     8 NI=0
111       V=Z
112       FZ=F(N1+1)
113       FF=FZ
114       F1=F(N1)
115       F2=F(N1-1)
116
117     9 NI=NI+1
118       IH=0
119       Z=V
120       F(N1+1)=FF
121       F(N1)=F1
122       F(N1-1)=F2
123    11 IH=IH+1
124       IF(IH .GT. NI) GO TO 9
125       CN=F(N1+1)-(NI-IH)*FZ/(NI+R0)
126       CD=F(N1)**2-2*F(N1-1)*CN
127       CN=CN*F(N1)
128       IF(ABS(CN) .GT. 10*ABS(CD)) GO TO 9
129       Z=Z-CN/CD
130       IF(M .GE. MAXIT) THEN
131        CALL MTLPRT(NAME,'C209.2','TOO MANY ITERATIONS')
132        RETURN
133       ENDIF
134
135       DO 71 I = 1,N1+1
136    71 F(I)=B(I)
137       ZA=ABS(Z)
138       HH=HF*ABS(B(1))
139       DO 72 I = 2,N1+1
140       F(I)=Z*F(I-1)+F(I)
141       D=ABS(F(I))
142    72 HH=ZA*HH+D
143       DO 73 L = 1,2
144       DO 73 I = 2,N1-L+1
145    73 F(I)=Z*F(I-1)+F(I)
146       M=M+1
147       ERR=EPS*(FAC2*HH-D)
148       IF(D .LE. ERR) GO TO 13
149       S=ABS(F(N1+1))
150       IF(S .LE. FAC1*D1) GO TO 7
151       IF(S .LT. BIG2*ERR) GO TO 13
152       S1=ABS(F(N1))
153       S2=ABS(F(N1-1))
154       S3=1+ABS(Z)
155       IF(2*S1 .GT. S2*S3 .OR. 10*S .LT. S1*S3) GO TO 11
156
157       ZZ=Z
158       KK=2
159       M1=1
160       S=S1/S2
161    40 DT=PI2/KK
162       DO 42 J=1,KK,M1
163       ZV=S*EXP(IU*(J*DT))
164       Z=ZZ+ZV
165       IF(M .GE. MAXIT) GO TO 44
166       DO 41 I = 1,N1+1
167    41 F(I)=B(I)
168       DO 43 L = 0,1
169       DO 43 I = 2,N1-L+1
170    43 F(I)=Z*F(I-1)+F(I)
171       M=M+1
172       IF(F(N1) .NE. 0) THEN
173        Z1=(ABS(F(N1))*F(N1+1))/F(N1)
174        IF(TST(Z1,ZV) .LT. 0) GO TO 44
175       ENDIF
176    42 CONTINUE
177       M1=2
178       KK=2*KK
179       IF(KK .LE. 8) GO TO 40
180
181    44 S2=1
182       S=ABS(F(N1+1))
183       S1=ABS(F(N1))
184       IF(S .GT. S1*(1+ABS(Z))) S2=S1/(2*S)
185       Z=Z-S2*F(N1+1)/F(N1)
186       IF(M .GE. MAXIT) THEN
187        CALL MTLPRT(NAME,'C209.2','TOO MANY ITERATIONS')
188        RETURN
189       ENDIF
190       DO 81 I = 1,N1+1
191    81 F(I)=B(I)
192       DO 83 L = 0,2
193       DO 83 I = 2,N1-L+1
194    83 F(I)=Z*F(I-1)+F(I)
195       M=M+1
196       GO TO 8
197
198    18 Z=0
199       H(K)=0
200
201    13 IF(Z .EQ. 0) GO TO 53
202       X=ABS(B(N1+1))
203       Z1=1
204       JR=0
205       DO 51 J = 2,N1
206       Z1=Z*Z1
207       Y=ABS(Z1*B(N1-J+2))
208       IF(Y .GT. X) THEN
209        X=Y
210        JR=J-1
211       ENDIF
212    51 CONTINUE
213       IF(JR .GT. 0) THEN
214        B(N1+1)=-B(N1+1)/Z
215        DO 52 J = N1,N1-JR+2,-1
216    52  B(J)=(B(J+1)-B(J))/Z
217        DO 55 J = N1-JR+1,N1
218    55  B(J)=B(J+1)
219       ENDIF
220       DO 54 J = 2,N1-JR
221    54 B(J)=B(J)+Z*B(J-1)
222
223    53 N1=N1-1
224       Z0(K)=Z
225    10 CONTINUE
226       Z0(N)=(-B(2)+SQRT(B(2)**2-4*B(1)*B(3)))/(2*B(1))
227       Z0(N-1)=-Z0(N)-B(2)/B(1)
228       IF(B(3) .EQ. 0) H(N)=0
229       IF(B(2) .EQ. 0 .AND. B(3) .EQ. 0) H(N-1)=0
230    12 CONTINUE
231
232       DO 14 K=1,N
233       IF(H(K) .EQ. 0) GO TO 14
234       H(K)=0
235
236       DO 91 I = 1,N+1
237    91 F(I)=C(I-1)
238       ZA=ABS(Z0(K))
239       HH=HF*ABS(C(0))
240       DO 92 I = 2,N+1
241       F(I)=Z0(K)*F(I-1)+F(I)
242       D=ABS(F(I))
243    92 HH=ZA*HH+D
244       DO 93 L = 1,2
245       DO 93 I = 2,N-L+1
246    93 F(I)=Z0(K)*F(I-1)+F(I)
247       IF(F(N+1) .EQ. 0 .AND. F(N) .EQ. 0) GO TO 14
248       ERR=EPS*(FAC2*HH-D)
249       IF(D .LE. ERR) THEN
250        S=WN*ERR
251       ELSE
252        S=WN*ABS(F(N+1))
253       ENDIF
254       S1=SQRT(ABS(F(N)**2-2*F(N-1)*F(N+1)))
255       IF(S .LE. 10*S1*ABS(Z0(K))) THEN
256        H(K)=S/S1
257       ELSE
258        WRITE(ERRTXT,103) M
259        CALL MTLPRT(NAME,'C209.3',ERRTXT)
260       ENDIF
261    14 CONTINUE
262       RETURN
263 103   FORMAT('BOUND FOR Z0',I3,' CANNOT BE FOUND')
264       END