Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/roots/roots.for
There is 1 other file named roots.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C ROOTS.F4 (FILENAME ON LIBRARY DECTAPE)
C ROOTS, 2.1.1 (CALLING NAME, SUBLST NO.)
C ROOTS OF POLYNOMIALS
C THIS PROGRAM CAME FROM DECUS(NO. 10-101, IBM SSP), AND
C WAS IMPLEMENTED BY N. GRANT AND B. GRANET.
C LIBRARY DECTAPE PROGS. USED: USAGE. MAC
C INTERNAL SUBR. USED: POLRT
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
DIMENSION XCOEF(31),COEF(31),ROOTR(30),ROOTI(30)
C CALL USAGE('ROOTS')
30 WRITE(5,2)
2 FORMAT('1ENTER DEGREE OF POLYNOMIAL. MAXIMUM OF 30'/)
20 READ(5,1) N
1 FORMAT(I)
IF(N .GT. 0 .AND. N .LE. 30) GO TO 10
WRITE(5,11)
11 FORMAT(' INVALID DEGREE. TRY AGAIN.'/)
GO TO 20
10 WRITE(5,3)
3 FORMAT(' ENTER COEFS. OF POLY. BEGINNING WITH CONSTANT.'/
1' MAXIMUM OF 10 PER LINE.'/)
NC=N+1
READ(5,4)(XCOEF(I),I=1,NC)
4 FORMAT(10G)
CALL POLRT(XCOEF,COEF,N,ROOTR,ROOTI,IER)
WRITE(5,7)
7 FORMAT('0INPUT COEFFICIENTS IN ASCENDING POWERS OF X'/)
WRITE(5,5)(XCOEF(I),I=1,NC)
5 FORMAT(' ',G)
IF(IER .EQ. 3) GO TO 13
IF(IER .EQ. 4) GO TO 12
WRITE(5,8)
8 FORMAT('0',10X,'REAL',13X,'IMAGINARY'/)
WRITE(5,6)(ROOTR(K),ROOTI(K),K=1,N)
6 FORMAT(1X/2(5X,G15.5))
GO TO 30
12 WRITE(5,14)
14 FORMAT(' HIGH ORDER COEFFICIENT IS ZERO')
GO TO 30
13 WRITE(5,15)
15 FORMAT(' ITERATIVE PROCEDURE FAILS.')
GO TO 30
END
C
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
C---------------XCOF, M ARE INPUT. OTHER ARGS. ARE RETURNED.
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
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