Power-Programmierung/CD1.mdf/fortran/library/ssp/polyroot/dprbm.for

Text File | 1985-11-29 | 8.1 KB | 252 lines

C
C     ..................................................................
C
C        SUBROUTINE DPRBM
C
C        PURPOSE
C           TO CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN
C           POLYNOMIAL WITH REAL COEFFICIENTS.
C
C        USAGE
C           CALL DPRBM (C,IC,RR,RC,POL,IR,IER)
C
C        DESCRIPTION OF PARAMETERS
C           C      - DOUBLE PRECISION INPUT VECTOR CONTAINING THE
C                    COEFFICIENTS OF THE GIVEN POLYNOMIAL. COEFFICIENTS
C                    ARE ORDERED FROM LOW TO HIGH. ON RETURN COEFFI-
C                    CIENTS ARE DIVIDED BY THE LAST NONZERO TERM.
C           IC     - DIMENSION OF VECTORS C, RR, RC, AND POL.
C           RR     - RESULTANT DOUBLE PRECISION VECTOR OF REAL PARTS
C                    OF THE ROOTS.
C           RC     - RESULTANT DOUBLE PRECISION VECTOR OF COMPLEX PARTS
C                    OF THE ROOTS.
C           POL    - RESULTANT DOUBLE PRECISION VECTOR OF COEFFICIENTS
C                    OF THE POLYNOMIAL WITH CALCULATED ROOTS.
C                    COEFFICIENTS ARE ORDERED FROM LOW TO HIGH (SEE
C                    REMARK 4).
C           IR     - OUTPUT VALUE SPECIFYING THE NUMBER OF CALCULATED
C                    ROOTS. NORMALLY IR IS EQUAL TO IC-1.
C           IER    - RESULTANT ERROR PARAMETER CODED AS FOLLOWS
C                     IER=0  - NO ERROR,
C                     IER=1  - SUBROUTINE DPQFB RECORDS POOR CONVERGENCE
C                              AT SOME QUADRATIC FACTORIZATION WITHIN
C                              100 ITERATION STEPS,
C                     IER=2  - POLYNOMIAL IS DEGENERATE, I.E. ZERO OR
C                              CONSTANT,
C                              OR OVERFLOW IN NORMALIZATION OF GIVEN
C                              POLYNOMIAL,
C                     IER=3  - THE SUBROUTINE IS BYPASSED DUE TO
C                              SUCCESSIVE ZERO DIVISORS OR OVERFLOWS
C                              IN QUADRATIC FACTORIZATION OR DUE TO
C                              COMPLETELY UNSATISFACTORY ACCURACY,
C                     IER=-1 - CALCULATED COEFFICIENT VECTOR HAS LESS
C                              THAN SIX CORRECT SIGNIFICANT DIGITS.
C                              THIS REVEALS POOR ACCURACY OF CALCULATED
C                              ROOTS.
C
C        REMARKS
C           (1) REAL PARTS OF THE ROOTS ARE STORED IN RR(1) UP TO RR(IR)
C               AND CORRESPONDING COMPLEX PARTS IN RC(1) UP TO RC(IR).
C           (2) ERROR MESSAGE IER=1 INDICATES POOR CONVERGENCE WITHIN
C               100 ITERATION STEPS AT SOME QUADRATIC FACTORIZATION
C               PERFORMED BY SUBROUTINE DPQFB.
C           (3) NO ACTION BESIDES ERROR MESSAGE IER=2 IN CASE OF A ZERO
C               OR CONSTANT POLYNOMIAL. THE SAME ERROR MESSAGE IS GIVEN
C               IN CASE OF AN OVERFLOW IN NORMALIZATION OF GIVEN
C               POLYNOMIAL.
C           (4) ERROR MESSAGE IER=3 INDICATES SUCCESSIVE ZERO DIVISORS
C               OR OVERFLOWS OR COMPLETELY UNSATISFACTORY ACCURACY AT
C               ANY QUADRATIC FACTORIZATION PERFORMED BY
C               SUBROUTINE DPQFB. IN THIS CASE CALCULATION IS BYPASSED.
C               IR RECORDS THE NUMBER OF CALCULATED ROOTS.
C               POL(1),...,POL(J-IR) ARE THE COEFFICIENTS OF THE
C               REMAINING POLYNOMIAL, WHERE J IS THE ACTUAL NUMBER OF
C               COEFFICIENTS IN VECTOR C (NORMALLY J=IC).
C           (5) IF CALCULATED COEFFICIENT VECTOR HAS LESS THAN SIX
C               CORRECT SIGNIFICANT DIGITS THOUGH ALL QUADRATIC
C               FACTORIZATIONS SHOWED SATISFACTORY ACCURACY, THE ERROR
C               MESSAGE IER=-1 IS GIVEN.
C           (6) THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED
C               COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE
C               BEEN CALCULATED. IN THIS CASE THE NUMBER OF ROOTS IR IS
C               EQUAL TO THE ACTUAL DEGREE OF THE POLYNOMIAL (NORMALLY
C               IR=IC-1). THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT
C               VECTOR IS RECORDED IN RR(IR+1).
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           SUBROUTINE DPQFB    QUADRATIC FACTORIZATION OF A POLYNOMIAL
C                               BY BAIRSTOW ITERATION.
C
C        METHOD
C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF
C           SUCCESSIVE QUADRATIC FACTORIZATION PERFORMED BY BAIRSTOW
C           ITERATION. X**2 IS USED AS INITIAL GUESS FOR THE FIRST
C           QUADRATIC FACTOR, AND FURTHER EACH CALCULATED QUADRATIC
C           FACTOR IS USED AS INITIAL GUESS FOR THE NEXT ONE. AFTER
C           COMPUTATION OF ALL ROOTS THE COEFFICIENT VECTOR IS
C           CALCULATED AND COMPARED WITH THE GIVEN ONE.
C           FOR REFERENCE, SEE J. H. WILKINSON, THE EVALUATION OF THE
C           ZEROS OF ILL-CONDITIONED POLYNOMIALS (PART ONE AND TWO),
C           NUMERISCHE MATHEMATIK, VOL.1 (1959), PP.150-180.
C
C     ..................................................................
C
      SUBROUTINE DPRBM(C,IC,RR,RC,POL,IR,IER)
C
C
      DIMENSION C(1),RR(1),RC(1),POL(1),Q(4)
      DOUBLE PRECISION C,RR,RC,POL,Q,EPS,A,B,H,Q1,Q2
C
C        TEST ON LEADING ZERO COEFFICIENTS
      EPS=1.D-6
      LIM=100
      IR=IC+1
    1 IR=IR-1
      IF(IR-1)42,42,2
    2 IF(C(IR))3,1,3
C
C        WORK UP ZERO ROOTS AND NORMALIZE REMAINING POLYNOMIAL
    3 IER=0
      J=IR
      L=0
      A=C(IR)
      DO 8 I=1,IR
      IF(L)4,4,7
    4 IF(C(I))6,5,6
    5 RR(I)=0.D0
      RC(I)=0.D0
      POL(J)=0.D0
      J=J-1
      GO TO 8
    6 L=1
      IST=I
      J=0
    7 J=J+1
      C(I)=C(I)/A
      POL(J)=C(I)
      CALL OVERFL(N)
      IF(N-2)42,8,8
    8 CONTINUE
C
C        START BAIRSTOW ITERATION
      Q1=0.D0
      Q2=0.D0
    9 IF(J-2)33,10,14
C
C        DEGREE OF RESTPOLYNOMIAL IS EQUAL TO ONE
   10 A=POL(1)
      RR(IST)=-A
      RC(IST)=0.D0
      IR=IR-1
      Q2=0.D0
      IF(IR-1)13,13,11
   11 DO 12 I=2,IR
      Q1=Q2
      Q2=POL(I+1)
   12 POL(I)=A*Q2+Q1
   13 POL(IR+1)=A+Q2
      GO TO 34
C        THIS IS BRANCH TO COMPARISON OF COEFFICIENT VECTORS C AND POL
C
C        DEGREE OF RESTPOLYNOMIAL IS GREATER THAN ONE
   14 DO 22 L=1,10
      N=1
   15 Q(1)=Q1
      Q(2)=Q2
      CALL DPQFB(POL,J,Q,LIM,I)
      IF(I)16,24,23
   16 IF(Q1)18,17,18
   17 IF(Q2)18,21,18
   18 GO TO (19,20,19,21),N
   19 Q1=-Q1
      N=N+1
      GO TO 15
   20 Q2=-Q2
      N=N+1
      GO TO 15
   21 Q1=1.D0+Q1
   22 Q2=1.D0-Q2
C
C        ERROR EXIT DUE TO UNSATISFACTORY RESULTS OF FACTORIZATION
      IER=3
      IR=IR-J
      RETURN
C
C        WORK UP RESULTS OF QUADRATIC FACTORIZATION
   23 IER=1
   24 Q1=Q(1)
      Q2=Q(2)
C
C        PERFORM DIVISION OF FACTORIZED POLYNOMIAL BY QUADRATIC FACTOR
      B=0.D0
      A=0.D0
      I=J
   25 H=-Q1*B-Q2*A+POL(I)
      POL(I)=B
      B=A
      A=H
      I=I-1
      IF(I-2)26,26,25
   26 POL(2)=B
      POL(1)=A
C
C        MULTIPLY POLYNOMIAL WITH CALCULATED ROOTS BY QUADRATIC FACTOR
      L=IR-1
      IF(J-L)27,27,29
   27 DO 28 I=J,L
   28 POL(I-1)=POL(I-1)+POL(I)*Q2+POL(I+1)*Q1
   29 POL(L)=POL(L)+POL(L+1)*Q2+Q1
      POL(IR)=POL(IR)+Q2
C
C        CALCULATE ROOT-PAIR FROM QUADRATIC FACTOR X*X+Q2*X+Q1
      H=-.5D0*Q2
      A=H*H-Q1
      B=DSQRT(DABS(A))
      IF(A)30,30,31
   30 RR(IST)=H
      RC(IST)=B
      IST=IST+1
      RR(IST)=H
      RC(IST)=-B
      GO TO 32
   31 B=H+DSIGN(B,H)
      RR(IST)=Q1/B
      RC(IST)=0.D0
      IST=IST+1
      RR(IST)=B
      RC(IST)=0.D0
   32 IST=IST+1
      J=J-2
      GO TO 9
C
C        SHIFT BACK ELEMENTS OF POL BY 1 AND COMPARE VECTORS POL AND C
   33 IR=IR-1
   34 A=0.D0
      DO 38 I=1,IR
      Q1=C(I)
      Q2=POL(I+1)
      POL(I)=Q2
      IF(Q1)35,36,35
   35 Q2=(Q1-Q2)/Q1
   36 Q2=DABS(Q2)
      IF(Q2-A)38,38,37
   37 A=Q2
   38 CONTINUE
      I=IR+1
      POL(I)=1.D0
      RR(I)=A
      RC(I)=0.D0
      IF(IER)39,39,41
   39 IF(A-EPS)41,41,40
C
C        WARNING DUE TO POOR ACCURACY OF CALCULATED COEFFICIENT VECTOR
   40 IER=-1
   41 RETURN
C
C        ERROR EXIT DUE TO DEGENERATE POLYNOMIAL OR OVERFLOW IN
C        NORMALIZATION
   42 IER=2
      IR=0
      RETURN
      END