File: POLRT.FT of Tape: Various/ETH/eth11-1
(Source file text)
C .................................................................. C C SUBROUTINE POLRT C C PURPOSE C COMPUTES THE REAL AND COMPLEX ROOTS OF A REAL POLYNOMIAL C C USAGE C CALL POLRT(XCOF,COF,M,ROOTR,ROOTI,IER) C C DESCRIPTION OF PARAMETERS C XCOF -VECTOR OF M+1 COEFFICIENTS OF THE POLYNOMIAL C ORDERED FROM SMALLEST TO LARGEST POWER C COF -WORKING VECTOR OF LENGTH M+1 C M -ORDER OF POLYNOMIAL C ROOTR-RESULTANT VECTOR OF LENGTH M CONTAINING REAL ROOTS C OF THE POLYNOMIAL C ROOTI-RESULTANT VECTOR OF LENGTH M CONTAINING THE C CORRESPONDING IMAGINARY ROOTS OF THE POLYNOMIAL C IER -ERROR CODE WHERE C IER=0 NO ERROR C IER=1 M LESS THAN ONE C IER=2 M GREATER THAN 36 C IER=3 UNABLE TO DETERMINE ROOT WITH 500 INTERATIONS C ON 5 STARTING VALUES C IER=4 HIGH ORDER COEFFICIENT IS ZERO C C REMARKS C LIMITED TO 36TH ORDER POLYNOMIAL OR LESS. C FLOATING POINT OVERFLOW MAY OCCUR FOR HIGH ORDER C POLYNOMIALS BUT WILL NOT AFFECT THE ACCURACY OF THE RESULTS. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C NONE C C METHOD C NEWTON-RAPHSON ITERATIVE TECHNIQUE. THE FINAL ITERATIONS C ON EACH ROOT ARE PERFORMED USING THE ORIGINAL POLYNOMIAL C RATHER THAN THE REDUCED POLYNOMIAL TO AVOID ACCUMULATED C ERRORS IN THE REDUCED POLYNOMIAL. C C .................................................................. C SUBROUTINE POLRT(XCOF,COF,M,ROOTR,ROOTI,IER) DIMENSION XCOF(1),COF(1),ROOTR(1),ROOTI(1) DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2,SUMSQ, 1 DX,DY,TEMP,ALPHA,DABS C C ............................................................... C C IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE C C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION C STATEMENT WHICH FOLLOWS. C C DOUBLE PRECISION XCOF,COF,ROOTR,ROOTI C C THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS C APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS C ROUTINE. C THE DOUBLE PRECISION VERSION MAY BE MODIFIED BY CHANGING THE C CONSTANT IN STATEMENT 78 TO 1.0D-12 AND IN STATEMENT 122 TO C 1.0D-10. THIS WILL PROVIDE HIGHER PRECISION RESULTS AT THE C COST OF EXECUTION TIME C C ............................................................... C IFIT=0 N=M IER=0 IF(XCOF(N+1))10,25,10 10 IF(N) 15,15,32 C C SET ERROR CODE TO 1 C 15 IER=1 20 RETURN C C SET ERROR CODE TO 4 C 25 IER=4 GO TO 20 C C SET ERROR CODE TO 2 C 30 IER=2 GO TO 20 32 IF(N-36) 35,35,30 35 NX=N NXX=N+1 N2=1 KJ1 = N+1 DO 40 L=1,KJ1 MT=KJ1-L+1 40 COF(MT)=XCOF(L) C C SET INITIAL VALUES C 45 XO=.00500101 YO=0.01000101 C C ZERO INITIAL VALUE COUNTER C IN=0 50 X=XO C C INCREMENT INITIAL VALUES AND COUNTER C XO=-10.0*YO YO=-10.0*X C C SET X AND Y TO CURRENT VALUE C X=XO Y=YO IN=IN+1 GO TO 59 55 IFIT=1 XPR=X YPR=Y C C EVALUATE POLYNOMIAL AND DERIVATIVES C 59 ICT=0 60 UX=0.0 UY=0.0 V =0.0 YT=0.0 XT=1.0 U=COF(N+1) IF(U) 65,130,65 65 DO 70 I=1,N L =N-I+1 TEMP=COF(L) XT2=X*XT-Y*YT YT2=X*YT+Y*XT U=U+TEMP*XT2 V=V+TEMP*YT2 FI=I UX=UX+FI*XT*TEMP UY=UY-FI*YT*TEMP XT=XT2 70 YT=YT2 SUMSQ=UX*UX+UY*UY IF(SUMSQ) 75,110,75 75 DX=(V*UY-U*UX)/SUMSQ X=X+DX DY=-(U*UY+V*UX)/SUMSQ Y=Y+DY 78 IF(DABS(DY)+DABS(DX)-1.0D-05) 100,80,80 C C STEP ITERATION COUNTER C 80 ICT=ICT+1 IF(ICT-500) 60,85,85 85 IF(IFIT)100,90,100 90 IF(IN-5) 50,95,95 C C SET ERROR CODE TO 3 C 95 IER=3 GO TO 20 100 DO 105 L=1,NXX MT=KJ1-L+1 TEMP=XCOF(MT) XCOF(MT)=COF(L) 105 COF(L)=TEMP ITEMP=N N=NX NX=ITEMP IF(IFIT) 120,55,120 110 IF(IFIT) 115,50,115 115 X=XPR Y=YPR 120 IFIT=0 122 IF(DABS(Y)-1.0D-4*DABS(X)) 135,125,125 125 ALPHA=X+X SUMSQ=X*X+Y*Y N=N-2 GO TO 140 130 X=0.0 NX=NX-1 NXX=NXX-1 135 Y=0.0 SUMSQ=0.0 ALPHA=X N=N-1 140 COF(2)=COF(2)+ALPHA*COF(1) 145 DO 150 L=2,N 150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1) 155 ROOTI(N2)=Y ROOTR(N2)=X N2=N2+1 IF(SUMSQ) 160,165,160 160 Y=-Y SUMSQ=0.0 GO TO 155 165 IF(N) 20,20,45 END C