Google
 

Trailing-Edge - PDP-10 Archives - decuslib10-02 - 43,50145/dprbm.ssp
There are 2 other files named dprbm.ssp in the archive. Click here to see a list.
C                                                                       DPRB  10
C     ..................................................................DPRB  20
C                                                                       DPRB  30
C        SUBROUTINE DPRBM                                               DPRB  40
C                                                                       DPRB  50
C        PURPOSE                                                        DPRB  60
C           TO CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN          DPRB  70
C           POLYNOMIAL WITH REAL COEFFICIENTS.                          DPRB  80
C                                                                       DPRB  90
C        USAGE                                                          DPRB 100
C           CALL DPRBM (C,IC,RR,RC,POL,IR,IER)                          DPRB 110
C                                                                       DPRB 120
C        DESCRIPTION OF PARAMETERS                                      DPRB 130
C           C      - DOUBLE PRECISION INPUT VECTOR CONTAINING THE       DPRB 140
C                    COEFFICIENTS OF THE GIVEN POLYNOMIAL. COEFFICIENTS DPRB 150
C                    ARE ORDERED FROM LOW TO HIGH. ON RETURN COEFFI-    DPRB 160
C                    CIENTS ARE DIVIDED BY THE LAST NONZERO TERM.       DPRB 170
C           IC     - DIMENSION OF VECTORS C, RR, RC, AND POL.           DPRB 180
C           RR     - RESULTANT DOUBLE PRECISION VECTOR OF REAL PARTS    DPRB 190
C                    OF THE ROOTS.                                      DPRB 200
C           RC     - RESULTANT DOUBLE PRECISION VECTOR OF COMPLEX PARTS DPRB 210
C                    OF THE ROOTS.                                      DPRB 220
C           POL    - RESULTANT DOUBLE PRECISION VECTOR OF COEFFICIENTS  DPRB 230
C                    OF THE POLYNOMIAL WITH CALCULATED ROOTS.           DPRB 240
C                    COEFFICIENTS ARE ORDERED FROM LOW TO HIGH (SEE     DPRB 250
C                    REMARK 4).                                         DPRB 260
C           IR     - OUTPUT VALUE SPECIFYING THE NUMBER OF CALCULATED   DPRB 270
C                    ROOTS. NORMALLY IR IS EQUAL TO IC-1.               DPRB 280
C           IER    - RESULTANT ERROR PARAMETER CODED AS FOLLOWS         DPRB 290
C                     IER=0  - NO ERROR,                                DPRB 300
C                     IER=1  - SUBROUTINE DPQFB RECORDS POOR CONVERGENCEDPRB 310
C                              AT SOME QUADRATIC FACTORIZATION WITHIN   DPRB 320
C                              100 ITERATION STEPS,                     DPRB 330
C                     IER=2  - POLYNOMIAL IS DEGENERATE, I.E. ZERO OR   DPRB 340
C                              CONSTANT,                                DPRB 350
C                              OR OVERFLOW IN NORMALIZATION OF GIVEN    DPRB 360
C                              POLYNOMIAL,                              DPRB 370
C                     IER=3  - THE SUBROUTINE IS BYPASSED DUE TO        DPRB 380
C                              SUCCESSIVE ZERO DIVISORS OR OVERFLOWS    DPRB 390
C                              IN QUADRATIC FACTORIZATION OR DUE TO     DPRB 400
C                              COMPLETELY UNSATISFACTORY ACCURACY,      DPRB 410
C                     IER=-1 - CALCULATED COEFFICIENT VECTOR HAS LESS   DPRB 420
C                              THAN SIX CORRECT SIGNIFICANT DIGITS.     DPRB 430
C                              THIS REVEALS POOR ACCURACY OF CALCULATED DPRB 440
C                              ROOTS.                                   DPRB 450
C                                                                       DPRB 460
C        REMARKS                                                        DPRB 470
C           (1) REAL PARTS OF THE ROOTS ARE STORED IN RR(1) UP TO RR(IR)DPRB 480
C               AND CORRESPONDING COMPLEX PARTS IN RC(1) UP TO RC(IR).  DPRB 490
C           (2) ERROR MESSAGE IER=1 INDICATES POOR CONVERGENCE WITHIN   DPRB 500
C               100 ITERATION STEPS AT SOME QUADRATIC FACTORIZATION     DPRB 510
C               PERFORMED BY SUBROUTINE DPQFB.                          DPRB 520
C           (3) NO ACTION BESIDES ERROR MESSAGE IER=2 IN CASE OF A ZERO DPRB 530
C               OR CONSTANT POLYNOMIAL. THE SAME ERROR MESSAGE IS GIVEN DPRB 540
C               IN CASE OF AN OVERFLOW IN NORMALIZATION OF GIVEN        DPRB 550
C               POLYNOMIAL.                                             DPRB 560
C           (4) ERROR MESSAGE IER=3 INDICATES SUCCESSIVE ZERO DIVISORS  DPRB 570
C               OR OVERFLOWS OR COMPLETELY UNSATISFACTORY ACCURACY AT   DPRB 580
C               ANY QUADRATIC FACTORIZATION PERFORMED BY                DPRB 590
C               SUBROUTINE DPQFB. IN THIS CASE CALCULATION IS BYPASSED. DPRB 600
C               IR RECORDS THE NUMBER OF CALCULATED ROOTS.              DPRB 610
C               POL(1),...,POL(J-IR) ARE THE COEFFICIENTS OF THE        DPRB 620
C               REMAINING POLYNOMIAL, WHERE J IS THE ACTUAL NUMBER OF   DPRB 630
C               COEFFICIENTS IN VECTOR C (NORMALLY J=IC).               DPRB 640
C           (5) IF CALCULATED COEFFICIENT VECTOR HAS LESS THAN SIX      DPRB 650
C               CORRECT SIGNIFICANT DIGITS THOUGH ALL QUADRATIC         DPRB 660
C               FACTORIZATIONS SHOWED SATISFACTORY ACCURACY, THE ERROR  DPRB 670
C               MESSAGE IER=-1 IS GIVEN.                                DPRB 680
C           (6) THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED       DPRB 690
C               COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE  DPRB 700
C               BEEN CALCULATED. IN THIS CASE THE NUMBER OF ROOTS IR IS DPRB 710
C               EQUAL TO THE ACTUAL DEGREE OF THE POLYNOMIAL (NORMALLY  DPRB 720
C               IR=IC-1). THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT DPRB 730
C               VECTOR IS RECORDED IN RR(IR+1).                         DPRB 740
C                                                                       DPRB 750
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DPRB 760
C           SUBROUTINE DPQFB    QUADRATIC FACTORIZATION OF A POLYNOMIAL DPRB 770
C                               BY BAIRSTOW ITERATION.                  DPRB 780
C                                                                       DPRB 790
C        METHOD                                                         DPRB 800
C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF      DPRB 810
C           SUCCESSIVE QUADRATIC FACTORIZATION PERFORMED BY BAIRSTOW    DPRB 820
C           ITERATION. X**2 IS USED AS INITIAL GUESS FOR THE FIRST      DPRB 830
C           QUADRATIC FACTOR, AND FURTHER EACH CALCULATED QUADRATIC     DPRB 840
C           FACTOR IS USED AS INITIAL GUESS FOR THE NEXT ONE. AFTER     DPRB 850
C           COMPUTATION OF ALL ROOTS THE COEFFICIENT VECTOR IS          DPRB 860
C           CALCULATED AND COMPARED WITH THE GIVEN ONE.                 DPRB 870
C           FOR REFERENCE, SEE J. H. WILKINSON, THE EVALUATION OF THE   DPRB 880
C           ZEROS OF ILL-CONDITIONED POLYNOMIALS (PART ONE AND TWO),    DPRB 890
C           NUMERISCHE MATHEMATIK, VOL.1 (1959), PP.150-180.            DPRB 900
C                                                                       DPRB 910
C     ..................................................................DPRB 920
C                                                                       DPRB 930
      SUBROUTINE DPRBM(C,IC,RR,RC,POL,IR,IER)                           DPRB 940
C                                                                       DPRB 950
C                                                                       DPRB 960
      DIMENSION C(1),RR(1),RC(1),POL(1),Q(4)                            DPRB 970
      DOUBLE PRECISION C,RR,RC,POL,Q,EPS,A,B,H,Q1,Q2                    DPRB 980
C                                                                       DPRB 990
C        TEST ON LEADING ZERO COEFFICIENTS                              DPRB1000
      EPS=1.D-6                                                         DPRB1010
      LIM=100                                                           DPRB1020
      IR=IC+1                                                           DPRB1030
    1 IR=IR-1                                                           DPRB1040
      IF(IR-1)42,42,2                                                   DPRB1050
    2 IF(C(IR))3,1,3                                                    DPRB1060
C                                                                       DPRB1070
C        WORK UP ZERO ROOTS AND NORMALIZE REMAINING POLYNOMIAL          DPRB1080
    3 IER=0                                                             DPRB1090
      J=IR                                                              DPRB1100
      L=0                                                               DPRB1110
      A=C(IR)                                                           DPRB1120
      DO 8 I=1,IR                                                       DPRB1130
      IF(L)4,4,7                                                        DPRB1140
    4 IF(C(I))6,5,6                                                     DPRB1150
    5 RR(I)=0.D0                                                        DPRB1160
      RC(I)=0.D0                                                        DPRB1170
      POL(J)=0.D0                                                       DPRB1180
      J=J-1                                                             DPRB1190
      GO TO 8                                                           DPRB1200
    6 L=1                                                               DPRB1210
      IST=I                                                             DPRB1220
      J=0                                                               DPRB1230
    7 J=J+1                                                             DPRB1240
      C(I)=C(I)/A                                                       DPRB1250
      POL(J)=C(I)                                                       DPRB1260
      CALL OVERFL(N)                                                    DPRB1270
      IF(N-2)42,8,8                                                     DPRB1280
    8 CONTINUE                                                          DPRB1290
C                                                                       DPRB1300
C        START BAIRSTOW ITERATION                                       DPRB1310
      Q1=0.D0                                                           DPRB1320
      Q2=0.D0                                                           DPRB1330
    9 IF(J-2)33,10,14                                                   DPRB1340
C                                                                       DPRB1350
C        DEGREE OF RESTPOLYNOMIAL IS EQUAL TO ONE                       DPRB1360
   10 A=POL(1)                                                          DPRB1370
      RR(IST)=-A                                                        DPRB1380
      RC(IST)=0.D0                                                      DPRB1390
      IR=IR-1                                                           DPRB1400
      Q2=0.D0                                                           DPRB1410
      IF(IR-1)13,13,11                                                  DPRB1420
   11 DO 12 I=2,IR                                                      DPRB1430
      Q1=Q2                                                             DPRB1440
      Q2=POL(I+1)                                                       DPRB1450
   12 POL(I)=A*Q2+Q1                                                    DPRB1460
   13 POL(IR+1)=A+Q2                                                    DPRB1470
      GO TO 34                                                          DPRB1480
C        THIS IS BRANCH TO COMPARISON OF COEFFICIENT VECTORS C AND POL  DPRB1490
C                                                                       DPRB1500
C        DEGREE OF RESTPOLYNOMIAL IS GREATER THAN ONE                   DPRB1510
   14 DO 22 L=1,10                                                      DPRB1520
      N=1                                                               DPRB1530
   15 Q(1)=Q1                                                           DPRB1540
      Q(2)=Q2                                                           DPRB1550
      CALL DPQFB(POL,J,Q,LIM,I)                                         DPRB1560
      IF(I)16,24,23                                                     DPRB1570
   16 IF(Q1)18,17,18                                                    DPRB1580
   17 IF(Q2)18,21,18                                                    DPRB1590
   18 GO TO (19,20,19,21),N                                             DPRB1600
   19 Q1=-Q1                                                            DPRB1610
      N=N+1                                                             DPRB1620
      GO TO 15                                                          DPRB1630
   20 Q2=-Q2                                                            DPRB1640
      N=N+1                                                             DPRB1650
      GO TO 15                                                          DPRB1660
   21 Q1=1.D0+Q1                                                        DPRB1670
   22 Q2=1.D0-Q2                                                        DPRB1680
C                                                                       DPRB1690
C        ERROR EXIT DUE TO UNSATISFACTORY RESULTS OF FACTORIZATION      DPRB1700
      IER=3                                                             DPRB1710
      IR=IR-J                                                           DPRB1720
      RETURN                                                            DPRB1730
C                                                                       DPRB1740
C        WORK UP RESULTS OF QUADRATIC FACTORIZATION                     DPRB1750
   23 IER=1                                                             DPRB1760
   24 Q1=Q(1)                                                           DPRB1770
      Q2=Q(2)                                                           DPRB1780
C                                                                       DPRB1790
C        PERFORM DIVISION OF FACTORIZED POLYNOMIAL BY QUADRATIC FACTOR  DPRB1800
      B=0.D0                                                            DPRB1810
      A=0.D0                                                            DPRB1820
      I=J                                                               DPRB1830
   25 H=-Q1*B-Q2*A+POL(I)                                               DPRB1840
      POL(I)=B                                                          DPRB1850
      B=A                                                               DPRB1860
      A=H                                                               DPRB1870
      I=I-1                                                             DPRB1880
      IF(I-2)26,26,25                                                   DPRB1890
   26 POL(2)=B                                                          DPRB1900
      POL(1)=A                                                          DPRB1910
C                                                                       DPRB1920
C        MULTIPLY POLYNOMIAL WITH CALCULATED ROOTS BY QUADRATIC FACTOR  DPRB1930
      L=IR-1                                                            DPRB1940
      IF(J-L)27,27,29                                                   DPRB1950
   27 DO 28 I=J,L                                                       DPRB1960
   28 POL(I-1)=POL(I-1)+POL(I)*Q2+POL(I+1)*Q1                           DPRB1970
   29 POL(L)=POL(L)+POL(L+1)*Q2+Q1                                      DPRB1980
      POL(IR)=POL(IR)+Q2                                                DPRB1990
C                                                                       DPRB2000
C        CALCULATE ROOT-PAIR FROM QUADRATIC FACTOR X*X+Q2*X+Q1          DPRB2010
      H=-.5D0*Q2                                                        DPRB2020
      A=H*H-Q1                                                          DPRB2030
      B=DSQRT(DABS(A))                                                  DPRB2040
      IF(A)30,30,31                                                     DPRB2050
   30 RR(IST)=H                                                         DPRB2060
      RC(IST)=B                                                         DPRB2070
      IST=IST+1                                                         DPRB2080
      RR(IST)=H                                                         DPRB2090
      RC(IST)=-B                                                        DPRB2100
      GO TO 32                                                          DPRB2110
   31 B=H+DSIGN(B,H)                                                    DPRB2120
      RR(IST)=Q1/B                                                      DPRB2130
      RC(IST)=0.D0                                                      DPRB2140
      IST=IST+1                                                         DPRB2150
      RR(IST)=B                                                         DPRB2160
      RC(IST)=0.D0                                                      DPRB2170
   32 IST=IST+1                                                         DPRB2180
      J=J-2                                                             DPRB2190
      GO TO 9                                                           DPRB2200
C                                                                       DPRB2210
C        SHIFT BACK ELEMENTS OF POL BY 1 AND COMPARE VECTORS POL AND C  DPRB2220
   33 IR=IR-1                                                           DPRB2230
   34 A=0.D0                                                            DPRB2240
      DO 38 I=1,IR                                                      DPRB2250
      Q1=C(I)                                                           DPRB2260
      Q2=POL(I+1)                                                       DPRB2270
      POL(I)=Q2                                                         DPRB2280
      IF(Q1)35,36,35                                                    DPRB2290
   35 Q2=(Q1-Q2)/Q1                                                     DPRB2300
   36 Q2=DABS(Q2)                                                       DPRB2310
      IF(Q2-A)38,38,37                                                  DPRB2320
   37 A=Q2                                                              DPRB2330
   38 CONTINUE                                                          DPRB2340
      I=IR+1                                                            DPRB2350
      POL(I)=1.D0                                                       DPRB2360
      RR(I)=A                                                           DPRB2370
      RC(I)=0.D0                                                        DPRB2380
      IF(IER)39,39,41                                                   DPRB2390
   39 IF(A-EPS)41,41,40                                                 DPRB2400
C                                                                       DPRB2410
C        WARNING DUE TO POOR ACCURACY OF CALCULATED COEFFICIENT VECTOR  DPRB2420
   40 IER=-1                                                            DPRB2430
   41 RETURN                                                            DPRB2440
C                                                                       DPRB2450
C        ERROR EXIT DUE TO DEGENERATE POLYNOMIAL OR OVERFLOW IN         DPRB2460
C        NORMALIZATION                                                  DPRB2470
   42 IER=2                                                             DPRB2480
      IR=0                                                              DPRB2490
      RETURN                                                            DPRB2500
      END                                                               DPRB2510