Google
 

Trailing-Edge - PDP-10 Archives - decuslib20-02 - decus/20-0026/prqd.ssp
There are 2 other files named prqd.ssp in the archive. Click here to see a list.
C                                                                       PRQD  10
C     ..................................................................PRQD  20
C                                                                       PRQD  30
C        SUBROUTINE PRQD                                                PRQD  40
C                                                                       PRQD  50
C        PURPOSE                                                        PRQD  60
C           CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN POLYNOMIAL  PRQD  70
C           WITH REAL COEFFICIENTS.                                     PRQD  80
C                                                                       PRQD  90
C        USAGE                                                          PRQD 100
C           CALL PRQD(C,IC,Q,E,POL,IR,IER)                              PRQD 110
C                                                                       PRQD 120
C        DESCRIPTION OF PARAMETERS                                      PRQD 130
C           C     - COEFFICIENT VECTOR OF GIVEN POLYNOMIAL              PRQD 140
C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH           PRQD 150
C                   THE GIVEN COEFFICIENT VECTOR GETS DIVIDED BY THE    PRQD 160
C                   LAST NONZERO TERM                                   PRQD 170
C           IC    - DIMENSION OF VECTOR C                               PRQD 180
C           Q     - WORKING STORAGE OF DIMENSION IC                     PRQD 190
C                   ON RETURN Q CONTAINS REAL PARTS OF ROOTS            PRQD 200
C           E     - WORKING STORAGE OF DIMENSION IC                     PRQD 210
C                   ON RETURN E CONTAINS COMPLEX PARTS OF ROOTS         PRQD 220
C           POL   - WORKING STORAGE OF DIMENSION IC                     PRQD 230
C                   ON RETURN POL CONTAINS THE COEFFICIENTS OF THE      PRQD 240
C                   POLYNOMIAL WITH CALCULATED ROOTS                    PRQD 250
C                   THIS RESULTING COEFFICIENT VECTOR HAS DIMENSION IR+1PRQD 260
C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH           PRQD 270
C           IR    - NUMBER OF CALCULATED ROOTS                          PRQD 280
C                   NORMALLY IR IS EQUAL TO DIMENSION IC MINUS ONE      PRQD 290
C           IER   - RESULTING ERROR PARAMETER. SEE REMARKS              PRQD 300
C                                                                       PRQD 310
C        REMARKS                                                        PRQD 320
C           THE REAL PART OF THE ROOTS IS STORED IN Q(1) UP TO Q(IR)    PRQD 330
C           CORRESPONDING COMPLEX PARTS ARE STORED IN E(1) UP TO E(IR). PRQD 340
C           IER = 0 MEANS NO ERRORS                                     PRQD 350
C           IER = 1 MEANS NO CONVERGENCE WITH FEASIBLE TOLERANCE        PRQD 360
C           IER = 2 MEANS POLYNOMIAL IS DEGENERATE (CONSTANT OR ZERO)   PRQD 370
C           IER = 3 MEANS SUBROUTINE WAS ABANDONED DUE TO ZERO DIVISOR  PRQD 380
C           IER = 4 MEANS THERE EXISTS NO S-FRACTION                    PRQD 390
C           IER =-1 MEANS CALCULATED COEFFICIENT VECTOR REVEALS POOR    PRQD 400
C                   ACCURACY OF THE CALCULATED ROOTS.                   PRQD 410
C                   THE CALCULATED COEFFICIENT VECTOR HAS LESS THAN     PRQD 420
C                   3 CORRECT DIGITS.                                   PRQD 430
C           THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED           PRQD 440
C           COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE BEEN PRQD 450
C           CALCULATED.                                                 PRQD 460
C           THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT VECTOR IS     PRQD 470
C           RECORDED IN Q(IR+1).                                        PRQD 480
C                                                                       PRQD 490
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  PRQD 500
C           NONE                                                        PRQD 510
C                                                                       PRQD 520
C        METHOD                                                         PRQD 530
C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF      PRQD 540
C           THE QUOTIENT-DIFFERENCE ALGORITHM WITH DISPLACEMENT.        PRQD 550
C           REFERENCE                                                   PRQD 560
C           H.RUTISHAUSER, DER QUOTIENTEN-DIFFERENZEN-ALGORITHMUS,      PRQD 570
C           BIRKHAEUSER, BASEL/STUTTGART, 1957.                         PRQD 580
C                                                                       PRQD 590
C     ..................................................................PRQD 600
C                                                                       PRQD 610
      SUBROUTINE PRQD(C,IC,Q,E,POL,IR,IER)                              PRQD 620
C                                                                       PRQD 630
C      DIMENSIONED DUMMY VARIABLES                                      PRQD 640
      DIMENSION E(1),Q(1),C(1),POL(1)                                   PRQD 650
C                                                                       PRQD 660
C        NORMALIZATION OF GIVEN POLYNOMIAL                              PRQD 670
C           TEST OF DIMENSION                                           PRQD 680
C        IR CONTAINS INDEX OF HIGHEST COEFFICIENT                       PRQD 690
      IER=0                                                             PRQD 700
      IR=IC                                                             PRQD 710
      EPS=1.E-6                                                         PRQD 720
      TOL=1.E-3                                                         PRQD 730
      LIMIT=10*IC                                                       PRQD 740
      KOUNT=0                                                           PRQD 750
    1 IF(IR-1)79,79,2                                                   PRQD 760
C                                                                       PRQD 770
C        DROP TRAILING ZERO COEFFICIENTS                                PRQD 780
    2 IF(C(IR))4,3,4                                                    PRQD 790
    3 IR=IR-1                                                           PRQD 800
      GOTO 1                                                            PRQD 810
C                                                                       PRQD 820
C           REARRANGEMENT OF GIVEN POLYNOMIAL                           PRQD 830
C        EXTRACTION OF ZERO ROOTS                                       PRQD 840
    4 O=1./C(IR)                                                        PRQD 850
      IEND=IR-1                                                         PRQD 860
      ISTA=1                                                            PRQD 870
      NSAV=IR+1                                                         PRQD 880
      JBEG=1                                                            PRQD 890
C                                                                       PRQD 900
C        Q(J)=1.                                                        PRQD 910
C        Q(J+I)=C(IR-I)/C(IR)                                           PRQD 920
C        Q(IR)=C(J)/C(IR)                                               PRQD 930
C        WHERE J IS THE INDEX OF THE LOWEST NONZERO COEFFICIENT         PRQD 940
      DO 9 I=1,IR                                                       PRQD 950
      J=NSAV-I                                                          PRQD 960
      IF(C(I))7,5,7                                                     PRQD 970
    5 GOTO(6,8),JBEG                                                    PRQD 980
    6 NSAV=NSAV+1                                                       PRQD 990
      Q(ISTA)=0.                                                        PRQD1000
      E(ISTA)=0.                                                        PRQD1010
      ISTA=ISTA+1                                                       PRQD1020
      GOTO 9                                                            PRQD1030
    7 JBEG=2                                                            PRQD1040
    8 Q(J)=C(I)*O                                                       PRQD1050
      C(I)=Q(J)                                                         PRQD1060
    9 CONTINUE                                                          PRQD1070
C                                                                       PRQD1080
C           INITIALIZATION                                              PRQD1090
      ESAV=0.                                                           PRQD1100
      Q(ISTA)=0.                                                        PRQD1110
   10 NSAV=IR                                                           PRQD1120
C                                                                       PRQD1130
C        COMPUTATION OF DERIVATIVE                                      PRQD1140
      EXPT=IR-ISTA                                                      PRQD1150
      E(ISTA)=EXPT                                                      PRQD1160
      DO 11 I=ISTA,IEND                                                 PRQD1170
      EXPT=EXPT-1.0                                                     PRQD1180
      POL(I+1)=EPS*ABS(Q(I+1))+EPS                                      PRQD1190
   11 E(I+1)=Q(I+1)*EXPT                                                PRQD1200
C                                                                       PRQD1210
C        TEST OF REMAINING DIMENSION                                    PRQD1220
      IF(ISTA-IEND)12,20,60                                             PRQD1230
   12 JEND=IEND-1                                                       PRQD1240
C                                                                       PRQD1250
C        COMPUTATION OF S-FRACTION                                      PRQD1260
      DO 19 I=ISTA,JEND                                                 PRQD1270
      IF(I-ISTA)13,16,13                                                PRQD1280
   13 IF(ABS(E(I))-POL(I+1))14,14,16                                    PRQD1290
C                                                                       PRQD1300
C        THE GIVEN POLYNOMIAL HAS MULTIPLE ROOTS, THE COEFFICIENTS OF   PRQD1310
C        THE COMMON FACTOR ARE STORED FROM Q(NSAV) UP TO Q(IR)          PRQD1320
   14 NSAV=I                                                            PRQD1330
      DO 15 K=I,JEND                                                    PRQD1340
      IF(ABS(E(K))-POL(K+1))15,15,80                                    PRQD1350
   15 CONTINUE                                                          PRQD1360
      GOTO 21                                                           PRQD1370
C                                                                       PRQD1380
C           EUCLIDEAN ALGORITHM                                         PRQD1390
   16 DO 19 K=I,IEND                                                    PRQD1400
      E(K+1)=E(K+1)/E(I)                                                PRQD1410
      Q(K+1)=E(K+1)-Q(K+1)                                              PRQD1420
      IF(K-I)18,17,18                                                   PRQD1430
C                                                                       PRQD1440
C        TEST FOR SMALL DIVISOR                                         PRQD1450
   17 IF(ABS(Q(I+1))-POL(I+1))80,80,19                                  PRQD1460
   18 Q(K+1)=Q(K+1)/Q(I+1)                                              PRQD1470
      POL(K+1)=POL(K+1)/ABS(Q(I+1))                                     PRQD1480
      E(K)=Q(K+1)-E(K)                                                  PRQD1490
   19 CONTINUE                                                          PRQD1500
   20 Q(IR)=-Q(IR)                                                      PRQD1510
C                                                                       PRQD1520
C           THE DISPLACEMENT EXPT IS SET TO 0 AUTOMATICALLY.            PRQD1530
C           E(ISTA)=0.,Q(ISTA+1),...,E(NSAV-1),Q(NSAV),E(NSAV)=0.,      PRQD1540
C           FORM A DIAGONAL OF THE QD-ARRAY.                            PRQD1550
C        INITIALIZATION OF BOUNDARY VALUES                              PRQD1560
   21 E(ISTA)=0.                                                        PRQD1570
      NRAN=NSAV-1                                                       PRQD1580
   22 E(NRAN+1)=0.                                                      PRQD1590
C                                                                       PRQD1600
C           TEST FOR LINEAR OR CONSTANT FACTOR                          PRQD1610
C        NRAN-ISTA IS DEGREE-1                                          PRQD1620
      IF(NRAN-ISTA)24,23,31                                             PRQD1630
C                                                                       PRQD1640
C        LINEAR FACTOR                                                  PRQD1650
   23 Q(ISTA+1)=Q(ISTA+1)+EXPT                                          PRQD1660
      E(ISTA+1)=0.                                                      PRQD1670
C                                                                       PRQD1680
C        TEST FOR UNFACTORED COMMON DIVISOR                             PRQD1690
   24 E(ISTA)=ESAV                                                      PRQD1700
      IF(IR-NSAV)60,60,25                                               PRQD1710
C                                                                       PRQD1720
C        INITIALIZE QD-ALGORITHM FOR COMMON DIVISOR                     PRQD1730
   25 ISTA=NSAV                                                         PRQD1740
      ESAV=E(ISTA)                                                      PRQD1750
      GOTO 10                                                           PRQD1760
C                                                                       PRQD1770
C        COMPUTATION OF ROOT PAIR                                       PRQD1780
   26 P=P+EXPT                                                          PRQD1790
C                                                                       PRQD1800
C           TEST FOR REALITY                                            PRQD1810
      IF(O)27,28,28                                                     PRQD1820
C                                                                       PRQD1830
C           COMPLEX ROOT PAIR                                           PRQD1840
   27 Q(NRAN)=P                                                         PRQD1850
      Q(NRAN+1)=P                                                       PRQD1860
      E(NRAN)=T                                                         PRQD1870
      E(NRAN+1)=-T                                                      PRQD1880
      GOTO 29                                                           PRQD1890
C                                                                       PRQD1900
C           REAL ROOT PAIR                                              PRQD1910
   28 Q(NRAN)=P-T                                                       PRQD1920
      Q(NRAN+1)=P+T                                                     PRQD1930
      E(NRAN)=0.                                                        PRQD1940
C                                                                       PRQD1950
C           REDUCTION OF DEGREE BY 2 (DEFLATION)                        PRQD1960
   29 NRAN=NRAN-2                                                       PRQD1970
      GOTO 22                                                           PRQD1980
C                                                                       PRQD1990
C        COMPUTATION OF REAL ROOT                                       PRQD2000
   30 Q(NRAN+1)=EXPT+P                                                  PRQD2010
C                                                                       PRQD2020
C           REDUCTION OF DEGREE BY 1 (DEFLATION)                        PRQD2030
      NRAN=NRAN-1                                                       PRQD2040
      GOTO 22                                                           PRQD2050
C                                                                       PRQD2060
C        START QD-ITERATION                                             PRQD2070
   31 JBEG=ISTA+1                                                       PRQD2080
      JEND=NRAN-1                                                       PRQD2090
      TEPS=EPS                                                          PRQD2100
      TDELT=1.E-2                                                       PRQD2110
   32 KOUNT=KOUNT+1                                                     PRQD2120
      P=Q(NRAN+1)                                                       PRQD2130
      R=ABS(E(NRAN))                                                    PRQD2140
C                                                                       PRQD2150
C           TEST FOR CONVERGENCE                                        PRQD2160
      IF(R-TEPS)30,30,33                                                PRQD2170
   33 S=ABS(E(JEND))                                                    PRQD2180
C                                                                       PRQD2190
C        IS THERE A REAL ROOT NEXT                                      PRQD2200
      IF(S-R)38,38,34                                                   PRQD2210
C                                                                       PRQD2220
C        IS DISPLACEMENT SMALL ENOUGH                                   PRQD2230
   34 IF(R-TDELT)36,35,35                                               PRQD2240
   35 P=0.                                                              PRQD2250
   36 O=P                                                               PRQD2260
      DO 37 J=JBEG,NRAN                                                 PRQD2270
      Q(J)=Q(J)+E(J)-E(J-1)-O                                           PRQD2280
C                                                                       PRQD2290
C           TEST FOR SMALL DIVISOR                                      PRQD2300
      IF(ABS(Q(J))-POL(J))81,81,37                                      PRQD2310
   37 E(J)=Q(J+1)*E(J)/Q(J)                                             PRQD2320
      Q(NRAN+1)=-E(NRAN)+Q(NRAN+1)-O                                    PRQD2330
      GOTO 54                                                           PRQD2340
C                                                                       PRQD2350
C        CALCULATE DISPLACEMENT FOR DOUBLE ROOTS                        PRQD2360
C           QUADRATIC EQUATION FOR DOUBLE ROOTS                         PRQD2370
C           X**2-(Q(NRAN)+Q(NRAN+1)+E(NRAN))*X+Q(NRAN)*Q(NRAN+1)=0      PRQD2380
   38 P=0.5*(Q(NRAN)+E(NRAN)+Q(NRAN+1))                                 PRQD2390
      O=P*P-Q(NRAN)*Q(NRAN+1)                                           PRQD2400
      T=SQRT(ABS(O))                                                    PRQD2410
C                                                                       PRQD2420
C        TEST FOR CONVERGENCE                                           PRQD2430
      IF(S-TEPS)26,26,39                                                PRQD2440
C                                                                       PRQD2450
C        ARE THERE COMPLEX ROOTS                                        PRQD2460
   39 IF(O)43,40,40                                                     PRQD2470
   40 IF(P)42,41,41                                                     PRQD2480
   41 T=-T                                                              PRQD2490
   42 P=P+T                                                             PRQD2500
      R=S                                                               PRQD2510
      GOTO 34                                                           PRQD2520
C                                                                       PRQD2530
C        MODIFICATION FOR COMPLEX ROOTS                                 PRQD2540
C        IS DISPLACEMENT SMALL ENOUGH                                   PRQD2550
   43 IF(S-TDELT)44,35,35                                               PRQD2560
C                                                                       PRQD2570
C        INITIALIZATION                                                 PRQD2580
   44 O=Q(JBEG)+E(JBEG)-P                                               PRQD2590
C                                                                       PRQD2600
C        TEST FOR SMALL DIVISOR                                         PRQD2610
      IF(ABS(O)-POL(JBEG))81,81,45                                      PRQD2620
   45 T=(T/O)**2                                                        PRQD2630
      U=E(JBEG)*Q(JBEG+1)/(O*(1.+T))                                    PRQD2640
      V=O+U                                                             PRQD2650
      KOUNT=KOUNT+2                                                     PRQD2660
C                                                                       PRQD2670
C        THREEFOLD LOOP FOR COMPLEX DISPLACEMENT                        PRQD2680
      DO 53 J=JBEG,NRAN                                                 PRQD2690
      O=Q(J+1)+E(J+1)-U-P                                               PRQD2700
C                                                                       PRQD2710
C        TEST FOR SMALL DIVISOR                                         PRQD2720
      IF(ABS(V)-POL(J))46,46,49                                         PRQD2730
   46 IF(J-NRAN)81,47,81                                                PRQD2740
   47 EXPT=EXPT+P                                                       PRQD2750
      IF(ABS(E(JEND))-TOL)48,48,81                                      PRQD2760
   48 P=0.5*(V+O-E(JEND))                                               PRQD2770
      O=P*P-(V-U)*(O-U*T-O*W*(1.+T)/Q(JEND))                            PRQD2780
      T=SQRT(ABS(O))                                                    PRQD2790
      GOTO 26                                                           PRQD2800
C                                                                       PRQD2810
C           TEST FOR SMALL DIVISOR                                      PRQD2820
   49 IF(ABS(O)-POL(J+1))46,46,50                                       PRQD2830
   50 W=U*O/V                                                           PRQD2840
      T=T*(V/O)**2                                                      PRQD2850
      Q(J)=V+W-E(J-1)                                                   PRQD2860
      U=0.                                                              PRQD2870
      IF(J-NRAN)51,52,52                                                PRQD2880
   51 U=Q(J+2)*E(J+1)/(O*(1.+T))                                        PRQD2890
   52 V=O+U-W                                                           PRQD2900
C                                                                       PRQD2910
C        TEST FOR SMALL DIVISOR                                         PRQD2920
      IF(ABS(Q(J))-POL(J))81,81,53                                      PRQD2930
   53 E(J)=W*V*(1.+T)/Q(J)                                              PRQD2940
      Q(NRAN+1)=V-E(NRAN)                                               PRQD2950
   54 EXPT=EXPT+P                                                       PRQD2960
      TEPS=TEPS*1.1                                                     PRQD2970
      TDELT=TDELT*1.1                                                   PRQD2980
      IF(KOUNT-LIMIT)32,55,55                                           PRQD2990
C                                                                       PRQD3000
C        NO CONVERGENCE WITH FEASIBLE TOLERANCE                         PRQD3010
C           ERROR RETURN IN CASE OF UNSATISFACTORY CONVERGENCE          PRQD3020
   55 IER=1                                                             PRQD3030
C                                                                       PRQD3040
C        REARRANGE CALCULATED ROOTS                                     PRQD3050
   56 IEND=NSAV-NRAN-1                                                  PRQD3060
      E(ISTA)=ESAV                                                      PRQD3070
      IF(IEND)59,59,57                                                  PRQD3080
   57 DO 58 I=1,IEND                                                    PRQD3090
      J=ISTA+I                                                          PRQD3100
      K=NRAN+1+I                                                        PRQD3110
      E(J)=E(K)                                                         PRQD3120
   58 Q(J)=Q(K)                                                         PRQD3130
   59 IR=ISTA+IEND                                                      PRQD3140
C                                                                       PRQD3150
C        NORMAL RETURN                                                  PRQD3160
   60 IR=IR-1                                                           PRQD3170
      IF(IR)78,78,61                                                    PRQD3180
C                                                                       PRQD3190
C        REARRANGE CALCULATED ROOTS                                     PRQD3200
   61 DO 62 I=1,IR                                                      PRQD3210
      Q(I)=Q(I+1)                                                       PRQD3220
   62 E(I)=E(I+1)                                                       PRQD3230
C                                                                       PRQD3240
C        CALCULATE COEFFICIENT VECTOR FROM ROOTS                        PRQD3250
      POL(IR+1)=1.                                                      PRQD3260
      IEND=IR-1                                                         PRQD3270
      JBEG=1                                                            PRQD3280
      DO 69 J=1,IR                                                      PRQD3290
      ISTA=IR+1-J                                                       PRQD3300
      O=0.                                                              PRQD3310
      P=Q(ISTA)                                                         PRQD3320
      T=E(ISTA)                                                         PRQD3330
      IF(T)65,63,65                                                     PRQD3340
C                                                                       PRQD3350
C        MULTIPLY WITH LINEAR FACTOR                                    PRQD3360
   63 DO 64 I=ISTA,IR                                                   PRQD3370
      POL(I)=O-P*POL(I+1)                                               PRQD3380
   64 O=POL(I+1)                                                        PRQD3390
      GOTO 69                                                           PRQD3400
   65 GOTO(66,67),JBEG                                                  PRQD3410
   66 JBEG=2                                                            PRQD3420
      POL(ISTA)=0.                                                      PRQD3430
      GOTO 69                                                           PRQD3440
C                                                                       PRQD3450
C        MULTIPLY WITH QUADRATIC FACTOR                                 PRQD3460
   67 JBEG=1                                                            PRQD3470
      U=P*P+T*T                                                         PRQD3480
      P=P+P                                                             PRQD3490
      DO 68 I=ISTA,IEND                                                 PRQD3500
      POL(I)=O-P*POL(I+1)+U*POL(I+2)                                    PRQD3510
   68 O=POL(I+1)                                                        PRQD3520
      POL(IR)=O-P                                                       PRQD3530
   69 CONTINUE                                                          PRQD3540
      IF(IER)78,70,78                                                   PRQD3550
C                                                                       PRQD3560
C        COMPARISON OF COEFFICIENT VECTORS, IE. TEST OF ACCURACY        PRQD3570
   70 P=0.                                                              PRQD3580
      DO 75 I=1,IR                                                      PRQD3590
      IF(C(I))72,71,72                                                  PRQD3600
   71 O=ABS(POL(I))                                                     PRQD3610
      GOTO 73                                                           PRQD3620
   72 O=ABS((POL(I)-C(I))/C(I))                                         PRQD3630
   73 IF(P-O)74,75,75                                                   PRQD3640
   74 P=O                                                               PRQD3650
   75 CONTINUE                                                          PRQD3660
      IF(P-TOL)77,76,76                                                 PRQD3670
   76 IER=-1                                                            PRQD3680
   77 Q(IR+1)=P                                                         PRQD3690
      E(IR+1)=0.                                                        PRQD3700
   78 RETURN                                                            PRQD3710
C                                                                       PRQD3720
C        ERROR RETURNS                                                  PRQD3730
C           ERROR RETURN FOR POLYNOMIALS OF DEGREE LESS THAN 1          PRQD3740
   79 IER=2                                                             PRQD3750
      IR=0                                                              PRQD3760
      RETURN                                                            PRQD3770
C                                                                       PRQD3780
C           ERROR RETURN IF THERE EXISTS NO S-FRACTION                  PRQD3790
   80 IER=4                                                             PRQD3800
      IR=ISTA                                                           PRQD3810
      GOTO 60                                                           PRQD3820
C                                                                       PRQD3830
C           ERROR RETURN IN CASE OF INSTABLE QD-ALGORITHM               PRQD3840
   81 IER=3                                                             PRQD3850
      GOTO 56                                                           PRQD3860
      END                                                               PRQD3870