Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap2_198111 - decus/20-0026/polrt.ssp
There are 2 other files named polrt.ssp in the archive. Click here to see a list.
C                                                                       PLRT  10
C     ..................................................................PLRT  20
C                                                                       PLRT  30
C        SUBROUTINE POLRT                                               PLRT  40
C                                                                       PLRT  50
C        PURPOSE                                                        PLRT  60
C           COMPUTES THE REAL AND COMPLEX ROOTS OF A REAL POLYNOMIAL    PLRT  70
C                                                                       PLRT  80
C        USAGE                                                          PLRT  90
C           CALL POLRT(XCOF,COF,M,ROOTR,ROOTI,IER)                      PLRT 100
C                                                                       PLRT 110
C        DESCRIPTION OF PARAMETERS                                      PLRT 120
C           XCOF -VECTOR OF M+1 COEFFICIENTS OF THE POLYNOMIAL          PLRT 130
C                 ORDERED FROM SMALLEST TO LARGEST POWER                PLRT 140
C           COF  -WORKING VECTOR OF LENGTH M+1                          PLRT 150
C           M    -ORDER OF POLYNOMIAL                                   PLRT 160
C           ROOTR-RESULTANT VECTOR OF LENGTH M CONTAINING REAL ROOTS    PLRT 170
C                 OF THE POLYNOMIAL                                     PLRT 180
C           ROOTI-RESULTANT VECTOR OF LENGTH M CONTAINING THE           PLRT 190
C                 CORRESPONDING IMAGINARY ROOTS OF THE POLYNOMIAL       PLRT 200
C           IER  -ERROR CODE WHERE                                      PLRT 210
C                 IER=0  NO ERROR                                       PLRT 220
C                 IER=1  M LESS THAN ONE                                PLRT 230
C                 IER=2  M GREATER THAN 36                              PLRT 240
C                 IER=3  UNABLE TO DETERMINE ROOT WITH 500 INTERATIONS  PLRT 250
C                        ON 5 STARTING VALUES                           PLRT 260
C                 IER=4  HIGH ORDER COEFFICIENT IS ZERO                 PLRT 270
C                                                                       PLRT 280
C        REMARKS                                                        PLRT 290
C           LIMITED TO 36TH ORDER POLYNOMIAL OR LESS.                   PLRT 300
C           FLOATING POINT OVERFLOW MAY OCCUR FOR HIGH ORDER            PLRT 310
C           POLYNOMIALS BUT WILL NOT AFFECT THE ACCURACY OF THE RESULTS.PLRT 320
C                                                                       PLRT 330
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  PLRT 340
C           NONE                                                        PLRT 350
C                                                                       PLRT 360
C        METHOD                                                         PLRT 370
C           NEWTON-RAPHSON ITERATIVE TECHNIQUE.  THE FINAL ITERATIONS   PLRT 380
C           ON EACH ROOT ARE PERFORMED USING THE ORIGINAL POLYNOMIAL    PLRT 390
C           RATHER THAN THE REDUCED POLYNOMIAL TO AVOID ACCUMULATED     PLRT 400
C           ERRORS IN THE REDUCED POLYNOMIAL.                           PLRT 410
C                                                                       PLRT 420
C     ..................................................................PLRT 430
C                                                                       PLRT 440
      SUBROUTINE POLRT(XCOF,COF,M,ROOTR,ROOTI,IER)                      PLRT 450
      DIMENSION XCOF(1),COF(1),ROOTR(1),ROOTI(1)                        PLRT 460
      DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2,SUMSQ, PLRT 470
     1 DX,DY,TEMP,ALPHA                                                 PLRT 480
C                                                                       PLRT 490
C        ...............................................................PLRT 500
C                                                                       PLRT 510
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE  PLRT 520
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION      PLRT 530
C        STATEMENT WHICH FOLLOWS.                                       PLRT 540
C                                                                       PLRT 550
C     DOUBLE PRECISION XCOF,COF,ROOTR,ROOTI                             PLRT 560
C                                                                       PLRT 570
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS    PLRT 580
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS      PLRT 590
C        ROUTINE.                                                       PLRT 600
C        THE DOUBLE PRECISION VERSION MAY BE MODIFIED BY CHANGING THE   PLRT 610
C        CONSTANT IN STATEMENT 78 TO 1.0D-12 AND IN STATEMENT 122 TO    PLRT 620
C        1.0D-10.  THIS WILL PROVIDE HIGHER PRECISION RESULTS AT THE    PLRT 630
C        COST OF EXECUTION TIME                                         PLRT 640
C                                                                       PLRT 650
C        ...............................................................PLRT 660
C                                                                       PLRT 670
      IFIT=0                                                            PLRT 680
      N=M                                                               PLRT 690
      IER=0                                                             PLRT 700
      IF(XCOF(N+1))10,25,10                                             PLRT 710
   10 IF(N) 15,15,32                                                    PLRT 720
C                                                                       PLRT 730
C        SET ERROR CODE TO 1                                            PLRT 740
C                                                                       PLRT 750
   15 IER=1                                                             PLRT 760
   20 RETURN                                                            PLRT 770
C                                                                       PLRT 780
C        SET ERROR CODE TO 4                                            PLRT 790
C                                                                       PLRT 800
   25 IER=4                                                             PLRT 810
      GO TO 20                                                          PLRT 820
C                                                                       PLRT 830
C        SET ERROR CODE TO 2                                            PLRT 840
C                                                                       PLRT 850
   30 IER=2                                                             PLRT 860
      GO TO 20                                                          PLRT 870
   32 IF(N-36) 35,35,30                                                 PLRT 880
   35 NX=N                                                              PLRT 890
      NXX=N+1                                                           PLRT 900
      N2=1                                                              PLRT 910
      KJ1 = N+1                                                         PLRT 920
      DO 40 L=1,KJ1                                                     PLRT 930
      MT=KJ1-L+1                                                        PLRT 940
   40 COF(MT)=XCOF(L)                                                   PLRT 950
C                                                                       PLRT 960
C        SET INITIAL VALUES                                             PLRT 970
C                                                                       PLRT 980
   45 XO=.00500101                                                      PLRT 990
      YO=0.01000101                                                     PLRT1000
C                                                                       PLRT1010
C        ZERO INITIAL VALUE COUNTER                                     PLRT1020
C                                                                       PLRT1030
      IN=0                                                              PLRT1040
   50 X=XO                                                              PLRT1050
C                                                                       PLRT1060
C        INCREMENT INITIAL VALUES AND COUNTER                           PLRT1070
C                                                                       PLRT1080
      XO=-10.0*YO                                                       PLRT1090
      YO=-10.0*X                                                        PLRT1100
C                                                                       PLRT1110
C        SET X AND Y TO CURRENT VALUE                                   PLRT1120
C                                                                       PLRT1130
      X=XO                                                              PLRT1140
      Y=YO                                                              PLRT1150
      IN=IN+1                                                           PLRT1160
      GO TO 59                                                          PLRT1170
   55 IFIT=1                                                            PLRT1180
      XPR=X                                                             PLRT1190
      YPR=Y                                                             PLRT1200
C                                                                       PLRT1210
C        EVALUATE POLYNOMIAL AND DERIVATIVES                            PLRT1220
C                                                                       PLRT1230
   59 ICT=0                                                             PLRT1240
   60 UX=0.0                                                            PLRT1250
      UY=0.0                                                            PLRT1260
      V =0.0                                                            PLRT1270
      YT=0.0                                                            PLRT1280
      XT=1.0                                                            PLRT1290
      U=COF(N+1)                                                        PLRT1300
      IF(U) 65,130,65                                                   PLRT1310
   65 DO 70 I=1,N                                                       PLRT1320
      L =N-I+1                                                          PLRT1330
      TEMP=COF(L)                                                       PLRT1340
      XT2=X*XT-Y*YT                                                     PLRT1350
      YT2=X*YT+Y*XT                                                     PLRT1360
      U=U+TEMP*XT2                                                      PLRT1370
      V=V+TEMP*YT2                                                      PLRT1380
      FI=I                                                              PLRT1390
      UX=UX+FI*XT*TEMP                                                  PLRT1400
      UY=UY-FI*YT*TEMP                                                  PLRT1410
      XT=XT2                                                            PLRT1420
   70 YT=YT2                                                            PLRT1430
      SUMSQ=UX*UX+UY*UY                                                 PLRT1440
      IF(SUMSQ) 75,110,75                                               PLRT1450
   75 DX=(V*UY-U*UX)/SUMSQ                                              PLRT1460
      X=X+DX                                                            PLRT1470
      DY=-(U*UY+V*UX)/SUMSQ                                             PLRT1480
      Y=Y+DY                                                            PLRT1490
   78 IF(DABS(DY)+DABS(DX)-1.0D-05) 100,80,80                           PLRT1500
C                                                                       PLRT1510
C        STEP ITERATION COUNTER                                         PLRT1520
C                                                                       PLRT1530
   80 ICT=ICT+1                                                         PLRT1540
      IF(ICT-500) 60,85,85                                              PLRT1550
   85 IF(IFIT)100,90,100                                                PLRT1560
   90 IF(IN-5) 50,95,95                                                 PLRT1570
C                                                                       PLRT1580
C        SET ERROR CODE TO 3                                            PLRT1590
C                                                                       PLRT1600
   95 IER=3                                                             PLRT1610
      GO TO 20                                                          PLRT1620
  100 DO 105 L=1,NXX                                                    PLRT1630
      MT=KJ1-L+1                                                        PLRT1640
      TEMP=XCOF(MT)                                                     PLRT1650
      XCOF(MT)=COF(L)                                                   PLRT1660
  105 COF(L)=TEMP                                                       PLRT1670
      ITEMP=N                                                           PLRT1680
      N=NX                                                              PLRT1690
      NX=ITEMP                                                          PLRT1700
      IF(IFIT) 120,55,120                                               PLRT1710
  110 IF(IFIT) 115,50,115                                               PLRT1720
  115 X=XPR                                                             PLRT1730
      Y=YPR                                                             PLRT1740
  120 IFIT=0                                                            PLRT1750
  122 IF(DABS(Y)-1.0D-4*DABS(X)) 135,125,125                            PLRT1760
  125 ALPHA=X+X                                                         PLRT1770
      SUMSQ=X*X+Y*Y                                                     PLRT1780
      N=N-2                                                             PLRT1790
      GO TO 140                                                         PLRT1800
  130 X=0.0                                                             PLRT1810
      NX=NX-1                                                           PLRT1820
      NXX=NXX-1                                                         PLRT1830
  135 Y=0.0                                                             PLRT1840
      SUMSQ=0.0                                                         PLRT1850
      ALPHA=X                                                           PLRT1860
      N=N-1                                                             PLRT1870
  140 COF(2)=COF(2)+ALPHA*COF(1)                                        PLRT1880
  145 DO 150 L=2,N                                                      PLRT1890
  150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1)                     PLRT1900
  155 ROOTI(N2)=Y                                                       PLRT1910
      ROOTR(N2)=X                                                       PLRT1920
      N2=N2+1                                                           PLRT1930
      IF(SUMSQ) 160,165,160                                             PLRT1940
  160 Y=-Y                                                              PLRT1950
      SUMSQ=0.0                                                         PLRT1960
      GO TO 155                                                         PLRT1970
  165 IF(N) 20,20,45                                                    PLRT1980
      END                                                               PLRT1990