| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140 |
- *DECK CPZERO
- SUBROUTINE CPZERO (IN, A, R, T, IFLG, S)
- C***BEGIN PROLOGUE CPZERO
- C***PURPOSE Find the zeros of a polynomial with complex coefficients.
- C***LIBRARY SLATEC
- C***CATEGORY F1A1B
- C***TYPE COMPLEX (RPZERO-S, CPZERO-C)
- C***KEYWORDS POLYNOMIAL ROOTS, POLYNOMIAL ZEROS, REAL ROOTS
- C***AUTHOR Kahaner, D. K., (NBS)
- C***DESCRIPTION
- C
- C Find the zeros of the complex polynomial
- C P(Z)= A(1)*Z**N + A(2)*Z**(N-1) +...+ A(N+1)
- C
- C Input...
- C IN = degree of P(Z)
- C A = complex vector containing coefficients of P(Z),
- C A(I) = coefficient of Z**(N+1-i)
- C R = N word complex vector containing initial estimates for zeros
- C if these are known.
- C T = 4(N+1) word array used for temporary storage
- C IFLG = flag to indicate if initial estimates of
- C zeros are input.
- C If IFLG .EQ. 0, no estimates are input.
- C If IFLG .NE. 0, the vector R contains estimates of
- C the zeros
- C ** WARNING ****** If estimates are input, they must
- C be separated, that is, distinct or
- C not repeated.
- C S = an N word array
- C
- C Output...
- C R(I) = Ith zero,
- C S(I) = bound for R(I) .
- C IFLG = error diagnostic
- C Error Diagnostics...
- C If IFLG .EQ. 0 on return, all is well
- C If IFLG .EQ. 1 on return, A(1)=0.0 or N=0 on input
- C If IFLG .EQ. 2 on return, the program failed to converge
- C after 25*N iterations. Best current estimates of the
- C zeros are in R(I). Error bounds are not calculated.
- C
- C***REFERENCES (NONE)
- C***ROUTINES CALLED CPEVL
- C***REVISION HISTORY (YYMMDD)
- C 810223 DATE WRITTEN
- C 890531 Changed all specific intrinsics to generic. (WRB)
- C 890531 REVISION DATE from Version 3.2
- C 891214 Prologue converted to Version 4.0 format. (BAB)
- C***END PROLOGUE CPZERO
- C
- REAL S(*)
- COMPLEX R(*),T(*),A(*),PN,TEMP
- C***FIRST EXECUTABLE STATEMENT CPZERO
- IF( IN .LE. 0 .OR. ABS(A(1)) .EQ. 0.0 ) GO TO 30
- C
- C CHECK FOR EASILY OBTAINED ZEROS
- C
- N=IN
- N1=N+1
- IF(IFLG .NE. 0) GO TO 14
- 1 N1=N+1
- IF(N .GT. 1) GO TO 2
- R(1)=-A(2)/A(1)
- S(1)=0.0
- RETURN
- 2 IF( ABS(A(N1)) .NE. 0.0 ) GO TO 3
- R(N)=0.0
- S(N)=0.0
- N=N-1
- GO TO 1
- C
- C IF INITIAL ESTIMATES FOR ZEROS NOT GIVEN, FIND SOME
- C
- 3 TEMP=-A(2)/(A(1)*N)
- CALL CPEVL(N,N,A,TEMP,T,T,.FALSE.)
- IMAX=N+2
- T(N1)=ABS(T(N1))
- DO 6 I=2,N1
- T(N+I)=-ABS(T(N+2-I))
- IF(REAL(T(N+I)) .LT. REAL(T(IMAX))) IMAX=N+I
- 6 CONTINUE
- X=(-REAL(T(IMAX))/REAL(T(N1)))**(1./(IMAX-N1))
- 7 X=2.*X
- CALL CPEVL(N,0,T(N1),CMPLX(X,0.0),PN,PN,.FALSE.)
- IF (REAL(PN).LT.0.) GO TO 7
- U=.5*X
- V=X
- 10 X=.5*(U+V)
- CALL CPEVL(N,0,T(N1),CMPLX(X,0.0),PN,PN,.FALSE.)
- IF (REAL(PN).GT.0.) V=X
- IF (REAL(PN).LE.0.) U=X
- IF((V-U) .GT. .001*(1.+V)) GO TO 10
- DO 13 I=1,N
- U=(3.14159265/N)*(2*I-1.5)
- 13 R(I)=MAX(X,.001*ABS(TEMP))*CMPLX(COS(U),SIN(U))+TEMP
- C
- C MAIN ITERATION LOOP STARTS HERE
- C
- 14 NR=0
- NMAX=25*N
- DO 19 NIT=1,NMAX
- DO 18 I=1,N
- IF(NIT .NE. 1 .AND. ABS(T(I)) .EQ. 0.) GO TO 18
- CALL CPEVL(N,0,A,R(I),PN,TEMP,.TRUE.)
- IF(ABS(REAL(PN))+ABS(AIMAG(PN)) .GT. REAL(TEMP)+
- 1 AIMAG(TEMP)) GO TO 16
- T(I)=0.0
- NR=NR+1
- GO TO 18
- 16 TEMP=A(1)
- DO 17 J=1,N
- 17 IF(J .NE. I) TEMP=TEMP*(R(I)-R(J))
- T(I)=PN/TEMP
- 18 CONTINUE
- DO 15 I=1,N
- 15 R(I)=R(I)-T(I)
- IF(NR .EQ. N) GO TO 21
- 19 CONTINUE
- GO TO 26
- C
- C CALCULATE ERROR BOUNDS FOR ZEROS
- C
- 21 DO 25 NR=1,N
- CALL CPEVL(N,N,A,R(NR),T,T(N+2),.TRUE.)
- X=ABS(CMPLX(ABS(REAL(T(1))),ABS(AIMAG(T(1))))+T(N+2))
- S(NR)=0.0
- DO 23 I=1,N
- X=X*REAL(N1-I)/I
- TEMP=CMPLX(MAX(ABS(REAL(T(I+1)))-REAL(T(N1+I)),0.0),
- 1 MAX(ABS(AIMAG(T(I+1)))-AIMAG(T(N1+I)),0.0))
- 23 S(NR)=MAX(S(NR),(ABS(TEMP)/X)**(1./I))
- 25 S(NR)=1./S(NR)
- RETURN
- C ERROR EXITS
- 26 IFLG=2
- RETURN
- 30 IFLG=1
- RETURN
- END
|
