home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / polyroot / dpqfb.for next >
Encoding:
Text File  |  1985-11-29  |  6.8 KB  |  241 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DPQFB
  5. C
  6. C        PURPOSE
  7. C           TO FIND AN APPROXIMATION Q(X)=Q1+Q2*X+X*X TO A QUADRATIC
  8. C           FACTOR OF A GIVEN POLYNOMIAL P(X) WITH REAL COEFFICIENTS.
  9. C
  10. C        USAGE
  11. C           CALL DPQFB(C,IC,Q,LIM,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           C   - DOUBLE PRECISION INPUT VECTOR CONTAINING THE
  15. C                 COEFFICIENTS OF P(X) - C(1) IS THE CONSTANT TERM
  16. C                 (DIMENSION IC)
  17. C           IC  - DIMENSION OF C
  18. C           Q   - DOUBLE PRECISION VECTOR OF DIMENSION 4 - ON INPUT Q(1)
  19. C                 AND Q(2) CONTAIN INITIAL GUESSES FOR Q1 AND Q2 - ON
  20. C                 RETURN Q(1) AND Q(2) CONTAIN THE REFINED COEFFICIENTS
  21. C                 Q1 AND Q2 OF Q(X), WHILE Q(3) AND Q(4) CONTAIN THE
  22. C                 COEFFICIENTS A AND B OF A+B*X, WHICH IS THE REMAINDER
  23. C                 OF THE QUOTIENT OF P(X) BY Q(X)
  24. C           LIM - INPUT VALUE SPECIFYING THE MAXIMUM NUMBER OF
  25. C                 ITERATIONS TO BE PERFORMED
  26. C           IER - RESULTING ERROR PARAMETER (SEE REMARKS)
  27. C                 IER= 0 - NO ERROR
  28. C                 IER= 1 - NO CONVERGENCE WITHIN LIM ITERATIONS
  29. C                 IER=-1 - THE POLYNOMIAL P(X) IS CONSTANT OR UNDEFINED
  30. C                          - OR OVERFLOW OCCURRED IN NORMALIZING P(X)
  31. C                 IER=-2 - THE POLYNOMIAL P(X) IS OF DEGREE 1
  32. C                 IER=-3 - NO FURTHER REFINEMENT OF THE APPROXIMATION TO
  33. C                          A QUADRATIC FACTOR IS FEASIBLE, DUE TO EITHER
  34. C                          DIVISION BY 0, OVERFLOW OR AN INITIAL GUESS
  35. C                          THAT IS NOT SUFFICIENTLY CLOSE TO A FACTOR OF
  36. C                          P(X)
  37. C
  38. C        REMARKS
  39. C           (1)  IF IER=-1 THERE IS NO COMPUTATION OTHER THAN THE
  40. C                POSSIBLE NORMALIZATION OF C.
  41. C           (2)  IF IER=-2 THERE IS NO COMPUTATION OTHER THAN THE
  42. C                NORMALIZATION OF C.
  43. C           (3)  IF IER =-3  IT IS SUGGESTED THAT A NEW INITIAL GUESS BE
  44. C                MADE FOR A QUADRATIC FACTOR.  Q, HOWEVER, WILL CONTAIN
  45. C                THE VALUES ASSOCIATED WITH THE ITERATION THAT YIELDED
  46. C                THE SMALLEST NORM OF THE MODIFIED LINEAR REMAINDER.
  47. C           (4)  IF IER=1, THEN, ALTHOUGH THE NUMBER OF ITERATIONS LIM
  48. C                WAS TOO SMALL TO INDICATE CONVERGENCE, NO OTHER PROB-
  49. C                LEMS HAVE BEEN DETECTED, AND Q WILL CONTAIN THE VALUES
  50. C                ASSOCIATED WITH THE ITERATION THAT YIELDED THE SMALLEST
  51. C                NORM OF THE MODIFIED LINEAR REMAINDER.
  52. C           (5)  FOR COMPLETE DETAIL SEE THE DOCUMENTATION FOR
  53. C                SUBROUTINES PQFB AND DPQFB.
  54. C
  55. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  56. C           NONE
  57. C
  58. C        METHOD
  59. C           COMPUTATION IS BASED ON BAIRSTOW'S ITERATIVE METHOD.  (SEE
  60. C           WILKINSON, J.H., THE EVALUATION OF THE ZEROS OF ILL-CON-
  61. C           DITIONED POLYNOMIALS (PART ONE AND TWO), NUMERISCHE MATHE-
  62. C           MATIK, VOL.1 (1959), PP. 150-180, OR HILDEBRAND, F.B.,
  63. C           INTRODUCTION TO NUMERICAL ANALYSIS, MC GRAW-HILL, NEW YORK/
  64. C           TORONTO/LONDON, 1956, PP. 472-476.)
  65. C
  66. C     ..................................................................
  67. C
  68.       SUBROUTINE DPQFB(C,IC,Q,LIM,IER)
  69. C
  70. C
  71.       DIMENSION C(1),Q(1)
  72.       DOUBLE PRECISION A,B,AA,BB,CA,CB,CC,CD,A1,B1,C1,H,HH,Q1,Q2,QQ1,
  73.      1                 QQ2,QQQ1,QQQ2,DQ1,DQ2,EPS,EPS1,C,Q
  74. C
  75. C        TEST ON LEADING ZERO COEFFICIENTS
  76.       IER=0
  77.       J=IC+1
  78.     1 J=J-1
  79.       IF(J-1)40,40,2
  80.     2 IF(C(J))3,1,3
  81. C
  82. C        NORMALIZATION OF REMAINING COEFFICIENTS
  83.     3 A=C(J)
  84.       IF(A-1.D0)4,6,4
  85.     4 DO 5 I=1,J
  86.       C(I)=C(I)/A
  87.       CALL OVERFL(N)
  88.       IF(N-2)40,5,5
  89.     5 CONTINUE
  90. C
  91. C        TEST ON NECESSITY OF BAIRSTOW ITERATION
  92.     6 IF(J-3)41,38,7
  93. C
  94. C        PREPARE BAIRSTOW ITERATION
  95.     7 EPS=1.D-14
  96.       EPS1=1.D-6
  97.       L=0
  98.       LL=0
  99.       Q1=Q(1)
  100.       Q2=Q(2)
  101.       QQ1=0.D0
  102.       QQ2=0.D0
  103.       AA=C(1)
  104.       BB=C(2)
  105.       CB=DABS(AA)
  106.       CA=DABS(BB)
  107.       IF(CB-CA)8,9,10
  108.     8 CC=CB+CB
  109.       CB=CB/CA
  110.       CA=1.D0
  111.       GO TO 11
  112.     9 CC=CA+CA
  113.       CA=1.D0
  114.       CB=1.D0
  115.       GO TO 11
  116.    10 CC=CA+CA
  117.       CA=CA/CB
  118.       CB=1.D0
  119.    11 CD=CC*.1D0
  120. C
  121. C        START BAIRSTOW ITERATION
  122. C        PREPARE NESTED MULTIPLICATION
  123.    12 A=0.D0
  124.       B=A
  125.       A1=A
  126.       B1=A
  127.       I=J
  128.       QQQ1=Q1
  129.       QQQ2=Q2
  130.       DQ1=HH
  131.       DQ2=H
  132. C
  133. C        START NESTED MULTIPLICATION
  134.    13 H=-Q1*B-Q2*A+C(I)
  135.       CALL OVERFL(N)
  136.       IF(N-2)42,14,14
  137.    14 B=A
  138.       A=H
  139.       I=I-1
  140.       IF(I-1)18,15,16
  141.    15 H=0.D0
  142.    16 H=-Q1*B1-Q2*A1+H
  143.       CALL OVERFL(N)
  144.       IF(N-2)42,17,17
  145.    17 C1=B1
  146.       B1=A1
  147.       A1=H
  148.       GO TO 13
  149. C        END OF NESTED MULTIPLICATION
  150. C
  151. C        TEST ON SATISFACTORY ACCURACY
  152.    18 H=CA*DABS(A)+CB*DABS(B)
  153.       IF(LL)19,19,39
  154.    19 L=L+1
  155.       IF(DABS(A)-EPS*DABS(C(1)))20,20,21
  156.    20 IF(DABS(B)-EPS*DABS(C(2)))39,39,21
  157. C
  158. C        TEST ON LINEAR REMAINDER OF MINIMUM NORM
  159.    21 IF(H-CC)22,22,23
  160.    22 AA=A
  161.       BB=B
  162.       CC=H
  163.       QQ1=Q1
  164.       QQ2=Q2
  165. C
  166. C        TEST ON LAST ITERATION STEP
  167.    23 IF(L-LIM)28,28,24
  168. C
  169. C        TEST ON RESTART OF BAIRSTOW ITERATION WITH ZERO INITIAL GUESS
  170.    24 IF(H-CD)43,43,25
  171.    25 IF(Q(1))27,26,27
  172.    26 IF(Q(2))27,42,27
  173.    27 Q(1)=0.D0
  174.       Q(2)=0.D0
  175.       GO TO 7
  176. C
  177. C        PERFORM ITERATION STEP
  178.    28 HH=DMAX1(DABS(A1),DABS(B1),DABS(C1))
  179.       IF(HH)42,42,29
  180.    29 A1=A1/HH
  181.       B1=B1/HH
  182.       C1=C1/HH
  183.       H=A1*C1-B1*B1
  184.       IF(H)30,42,30
  185.    30 A=A/HH
  186.       B=B/HH
  187.       HH=(B*A1-A*B1)/H
  188.       H=(A*C1-B*B1)/H
  189.       Q1=Q1+HH
  190.       Q2=Q2+H
  191. C        END OF ITERATION STEP
  192. C
  193. C        TEST ON SATISFACTORY RELATIVE ERROR OF ITERATED VALUES
  194.       IF(DABS(HH)-EPS*DABS(Q1))31,31,33
  195.    31 IF(DABS(H)-EPS*DABS(Q2))32,32,33
  196.    32 LL=1
  197.       GO TO 12
  198. C
  199. C        TEST ON DECREASING RELATIVE ERRORS
  200.    33 IF(L-1)12,12,34
  201.    34 IF(DABS(HH)-EPS1*DABS(Q1))35,35,12
  202.    35 IF(DABS(H)-EPS1*DABS(Q2))36,36,12
  203.    36 IF(DABS(QQQ1*HH)-DABS(Q1*DQ1))37,44,44
  204.    37 IF(DABS(QQQ2*H)-DABS(Q2*DQ2))12,44,44
  205. C        END OF BAIRSTOW ITERATION
  206. C
  207. C        EXIT IN CASE OF QUADRATIC POLYNOMIAL
  208.    38 Q(1)=C(1)
  209.       Q(2)=C(2)
  210.       Q(3)=0.D0
  211.       Q(4)=0.D0
  212.       RETURN
  213. C
  214. C        EXIT IN CASE OF SUFFICIENT ACCURACY
  215.    39 Q(1)=Q1
  216.       Q(2)=Q2
  217.       Q(3)=A
  218.       Q(4)=B
  219.       RETURN
  220. C
  221. C        ERROR EXIT IN CASE OF ZERO OR CONSTANT POLYNOMIAL
  222.    40 IER=-1
  223.       RETURN
  224. C
  225. C        ERROR EXIT IN CASE OF LINEAR POLYNOMIAL
  226.    41 IER=-2
  227.       RETURN
  228. C
  229. C        ERROR EXIT IN CASE OF NONREFINED QUADRATIC FACTOR
  230.    42 IER=-3
  231.       GO TO 44
  232. C
  233. C        ERROR EXIT IN CASE OF UNSATISFACTORY ACCURACY
  234.    43 IER=1
  235.    44 Q(1)=QQ1
  236.       Q(2)=QQ2
  237.       Q(3)=AA
  238.       Q(4)=BB
  239.       RETURN
  240.       END
  241.