+++ /dev/null
-*
-* $Id$
-*
-* $Log$
-* Revision 1.1.1.1 1996/04/01 15:01:53 mclareni
-* Mathlib gen
-*
-*
-#include "gen/pilot.h"
-#if defined(CERNLIB_DOUBLE)
- SUBROUTINE WPOLYZ(C,N,MAXIT,Z0,H)
-C
-#include "gen/imp64.inc"
-C
- CHARACTER NAME*(*)
- PARAMETER (NAME = 'WPOLYZ')
-#endif
-#if !defined(CERNLIB_DOUBLE)
- SUBROUTINE CPOLYZ(C,N,MAXIT,Z0,H)
-C
- CHARACTER NAME*(*)
- PARAMETER (NAME = 'CPOLYZ')
-#endif
-C
-C Based on
-C T. Pomentale, Homotopy methods for polynomial equations,
-C J. Inst. Maths Applics 13 (1974) 201-213
-
-C Revised version of a routine originally written by T. Pomentale
-
-#include "gen/defc64.inc"
- + C(0:*),Z0(*),F(100),B(100),Z,CN,CD,FZ,FF,F1,F2,V,ZZ,Z1,ZV,IU
- LOGICAL LMD
- DIMENSION H(*)
- CHARACTER*80 ERRTXT
-
- PARAMETER (EPS = 1D-14)
-#if defined(CERNLIB_CRAY)
- PARAMETER (IU = (0E0,1E0))
-#endif
-#if !defined(CERNLIB_CRAY)
- PARAMETER (IU = (0D0,1D0))
-#endif
- PARAMETER (R0 = 0, R1 = 1, BIG1 = 1D7, BIG2 = BIG1/10, HF = R1/2)
- PARAMETER (FAC1 = 9*R1/10, FAC2 = 2+100*FAC1*EPS)
- PARAMETER (PI = 3.14159 26535 89793 24D0, PI2 = 2*PI)
-
-#if defined(CERNLIB_QF2C)
-#include "defdr.inc"
-#endif
-#if defined(CERNLIB_DOUBLE)
- TST(Z1,ZV)=DREAL(Z1)*DREAL(ZV)+DIMAG(Z1)*DIMAG(ZV)
-#endif
-#if !defined(CERNLIB_DOUBLE)
- TST(Z1,ZV)=REAL(Z1)*REAL(ZV)+AIMAG(Z1)*AIMAG(ZV)
-#endif
-
- IF(N .LE. 0) RETURN
- IF(C(0) .EQ. 0) THEN
- CALL MTLPRT(NAME,'C209.1','A(0) = 0')
- RETURN
- ENDIF
- IF(N .EQ. 1) THEN
- Z0(1)=-C(1)/C(0)
- H(1)=0
- RETURN
- ENDIF
- IF(N .EQ. 2) THEN
- H(1)=0
- H(2)=0
- Z0(2)=(-C(1)+SQRT(C(1)**2-4*C(0)*C(2)))/(2*C(0))
- Z0(1)=-Z0(2)-C(1)/C(0)
- RETURN
- ENDIF
- LMD=.FALSE.
- DO 1 K = 1,N
- 1 LMD=LMD .OR. Z0(K) .NE. 0
- M=0
- WN=SQRT(N+R0)
- Z=1+IU
- N1=N
-
- DO 3 K = 1,N+1
- H(K)=1
- 3 B(K)=C(K-1)
- 6 DO 10 K = 1,N-2
- IF(B(N1+1) .EQ. 0) GO TO 18
- IF(LMD) Z=Z0(K)
- IF(M .GE. MAXIT) THEN
- CALL MTLPRT(NAME,'C209.2','TOO MANY ITERATIONS')
- RETURN
- ENDIF
-
- S=BIG1
- DO 61 I = 1,N1+1
- 61 F(I)=B(I)
- ZA=ABS(Z)
- HH=HF*ABS(B(1))
- DO 62 I = 2,N1+1
- F(I)=Z*F(I-1)+F(I)
- D=ABS(F(I))
- 62 HH=ZA*HH+D
- DO 63 L = 1,2
- DO 63 I = 2,N1-L+1
- 63 F(I)=Z*F(I-1)+F(I)
- M=M+1
- IF(D .LE. EPS*(FAC2*HH-D)) GO TO 13
-
- 7 D1=S
- 8 NI=0
- V=Z
- FZ=F(N1+1)
- FF=FZ
- F1=F(N1)
- F2=F(N1-1)
-
- 9 NI=NI+1
- IH=0
- Z=V
- F(N1+1)=FF
- F(N1)=F1
- F(N1-1)=F2
- 11 IH=IH+1
- IF(IH .GT. NI) GO TO 9
- CN=F(N1+1)-(NI-IH)*FZ/(NI+R0)
- CD=F(N1)**2-2*F(N1-1)*CN
- CN=CN*F(N1)
- IF(ABS(CN) .GT. 10*ABS(CD)) GO TO 9
- Z=Z-CN/CD
- IF(M .GE. MAXIT) THEN
- CALL MTLPRT(NAME,'C209.2','TOO MANY ITERATIONS')
- RETURN
- ENDIF
-
- DO 71 I = 1,N1+1
- 71 F(I)=B(I)
- ZA=ABS(Z)
- HH=HF*ABS(B(1))
- DO 72 I = 2,N1+1
- F(I)=Z*F(I-1)+F(I)
- D=ABS(F(I))
- 72 HH=ZA*HH+D
- DO 73 L = 1,2
- DO 73 I = 2,N1-L+1
- 73 F(I)=Z*F(I-1)+F(I)
- M=M+1
- ERR=EPS*(FAC2*HH-D)
- IF(D .LE. ERR) GO TO 13
- S=ABS(F(N1+1))
- IF(S .LE. FAC1*D1) GO TO 7
- IF(S .LT. BIG2*ERR) GO TO 13
- S1=ABS(F(N1))
- S2=ABS(F(N1-1))
- S3=1+ABS(Z)
- IF(2*S1 .GT. S2*S3 .OR. 10*S .LT. S1*S3) GO TO 11
-
- ZZ=Z
- KK=2
- M1=1
- S=S1/S2
- 40 DT=PI2/KK
- DO 42 J=1,KK,M1
- ZV=S*EXP(IU*(J*DT))
- Z=ZZ+ZV
- IF(M .GE. MAXIT) GO TO 44
- DO 41 I = 1,N1+1
- 41 F(I)=B(I)
- DO 43 L = 0,1
- DO 43 I = 2,N1-L+1
- 43 F(I)=Z*F(I-1)+F(I)
- M=M+1
- IF(F(N1) .NE. 0) THEN
- Z1=(ABS(F(N1))*F(N1+1))/F(N1)
- IF(TST(Z1,ZV) .LT. 0) GO TO 44
- ENDIF
- 42 CONTINUE
- M1=2
- KK=2*KK
- IF(KK .LE. 8) GO TO 40
-
- 44 S2=1
- S=ABS(F(N1+1))
- S1=ABS(F(N1))
- IF(S .GT. S1*(1+ABS(Z))) S2=S1/(2*S)
- Z=Z-S2*F(N1+1)/F(N1)
- IF(M .GE. MAXIT) THEN
- CALL MTLPRT(NAME,'C209.2','TOO MANY ITERATIONS')
- RETURN
- ENDIF
- DO 81 I = 1,N1+1
- 81 F(I)=B(I)
- DO 83 L = 0,2
- DO 83 I = 2,N1-L+1
- 83 F(I)=Z*F(I-1)+F(I)
- M=M+1
- GO TO 8
-
- 18 Z=0
- H(K)=0
-
- 13 IF(Z .EQ. 0) GO TO 53
- X=ABS(B(N1+1))
- Z1=1
- JR=0
- DO 51 J = 2,N1
- Z1=Z*Z1
- Y=ABS(Z1*B(N1-J+2))
- IF(Y .GT. X) THEN
- X=Y
- JR=J-1
- ENDIF
- 51 CONTINUE
- IF(JR .GT. 0) THEN
- B(N1+1)=-B(N1+1)/Z
- DO 52 J = N1,N1-JR+2,-1
- 52 B(J)=(B(J+1)-B(J))/Z
- DO 55 J = N1-JR+1,N1
- 55 B(J)=B(J+1)
- ENDIF
- DO 54 J = 2,N1-JR
- 54 B(J)=B(J)+Z*B(J-1)
-
- 53 N1=N1-1
- Z0(K)=Z
- 10 CONTINUE
- Z0(N)=(-B(2)+SQRT(B(2)**2-4*B(1)*B(3)))/(2*B(1))
- Z0(N-1)=-Z0(N)-B(2)/B(1)
- IF(B(3) .EQ. 0) H(N)=0
- IF(B(2) .EQ. 0 .AND. B(3) .EQ. 0) H(N-1)=0
- 12 CONTINUE
-
- DO 14 K=1,N
- IF(H(K) .EQ. 0) GO TO 14
- H(K)=0
-
- DO 91 I = 1,N+1
- 91 F(I)=C(I-1)
- ZA=ABS(Z0(K))
- HH=HF*ABS(C(0))
- DO 92 I = 2,N+1
- F(I)=Z0(K)*F(I-1)+F(I)
- D=ABS(F(I))
- 92 HH=ZA*HH+D
- DO 93 L = 1,2
- DO 93 I = 2,N-L+1
- 93 F(I)=Z0(K)*F(I-1)+F(I)
- IF(F(N+1) .EQ. 0 .AND. F(N) .EQ. 0) GO TO 14
- ERR=EPS*(FAC2*HH-D)
- IF(D .LE. ERR) THEN
- S=WN*ERR
- ELSE
- S=WN*ABS(F(N+1))
- ENDIF
- S1=SQRT(ABS(F(N)**2-2*F(N-1)*F(N+1)))
- IF(S .LE. 10*S1*ABS(Z0(K))) THEN
- H(K)=S/S1
- ELSE
- WRITE(ERRTXT,103) M
- CALL MTLPRT(NAME,'C209.3',ERRTXT)
- ENDIF
- 14 CONTINUE
- RETURN
-103 FORMAT('BOUND FOR Z0',I3,' CANNOT BE FOUND')
- END
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff