home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / polyroot / dprqd.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  10.1 KB  |  393 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DPRQD
  5. C
  6. C        PURPOSE
  7. C           CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN POLYNOMIAL
  8. C           WITH REAL COEFFICIENTS.
  9. C
  10. C        USAGE
  11. C           CALL DPRQD(C,IC,Q,E,POL,IR,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           C     - COEFFICIENT VECTOR OF GIVEN POLYNOMIAL
  15. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  16. C                   THE GIVEN COEFFICIENT VECTOR GETS DIVIDED BY THE
  17. C                   LAST NONZERO TERM
  18. C                   DOUBLE PRECISION ARRAY
  19. C           IC    - DIMENSION OF VECTOR C
  20. C           Q     - WORKING STORAGE OF DIMENSION IC
  21. C                   ON RETURN Q CONTAINS REAL PARTS OF ROOTS
  22. C                   DOUBLE PRECISION ARRAY
  23. C           E     - WORKING STORAGE OF DIMENSION IC
  24. C                   ON RETURN E CONTAINS COMPLEX PARTS OF ROOTS
  25. C                   DOUBLE PRECISION ARRAY
  26. C           POL   - WORKING STORAGE OF DIMENSION IC
  27. C                   ON RETURN POL CONTAINS THE COEFFICIENTS OF THE
  28. C                   POLYNOMIAL WITH CALCULATED ROOTS
  29. C                   THIS RESULTING COEFFICIENT VECTOR HAS DIMENSION IR+1
  30. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  31. C                   DOUBLE PRECISION ARRAY
  32. C           IR    - NUMBER OF CALCULATED ROOTS
  33. C                   NORMALLY IR IS EQUAL TO DIMENSION IC MINUS ONE
  34. C           IER   - RESULTING ERROR PARAMETER. SEE REMARKS
  35. C
  36. C        REMARKS
  37. C           THE REAL PART OF THE ROOTS IS STORED IN Q(1) UP TO Q(IR)
  38. C           CORRESPONDING COMPLEX PARTS ARE STORED IN E(1) UP TO E(IR).
  39. C           IER = 0 MEANS NO ERRORS
  40. C           IER = 1 MEANS NO CONVERGENCE WITH FEASIBLE TOLERANCE
  41. C           IER = 2 MEANS POLYNOMIAL IS DEGENERATE (CONSTANT OR ZERO)
  42. C           IER = 3 MEANS SUBROUTINE WAS ABANDONED DUE TO ZERO DIVISOR
  43. C           IER = 4 MEANS THERE EXISTS NO S-FRACTION
  44. C           IER =-1 MEANS CALCULATED COEFFICIENT VECTOR REVEALS POOR
  45. C                   ACCURACY OF THE CALCULATED ROOTS.
  46. C                   THE CALCULATED COEFFICIENT VECTOR HAS LESS THAN
  47. C                   6 CORRECT DIGITS.
  48. C           THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED
  49. C           COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE BEEN
  50. C           CALCULATED.
  51. C           THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT VECTOR IS
  52. C           RECORDED IN Q(IR+1).
  53. C
  54. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  55. C           NONE
  56. C
  57. C        METHOD
  58. C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF
  59. C           THE QUOTIENT-DIFFERENCE ALGORITHM WITH DISPLACEMENT.
  60. C           REFERENCE
  61. C           H.RUTISHAUSER, DER QUOTIENTEN-DIFFERENZEN-ALGORITHMUS,
  62. C           BIRKHAEUSER, BASEL/STUTTGART, 1957.
  63. C
  64. C     ..................................................................
  65. C
  66.       SUBROUTINE DPRQD(C,IC,Q,E,POL,IR,IER)
  67. C
  68. C      DIMENSIONED DUMMY VARIABLES
  69.       DIMENSION E(1),Q(1),C(1),POL(1)
  70.       DOUBLE PRECISION Q,E,O,P,T,EXPT,ESAV,U,V,W,C,POL,EPS
  71. C
  72. C        NORMALIZATION OF GIVEN POLYNOMIAL
  73. C           TEST OF DIMENSION
  74. C        IR CONTAINS INDEX OF HIGHEST COEFFICIENT
  75.       IR=IC
  76.       IER=0
  77.       EPS=1.D-16
  78.       TOL=1.E-6
  79.       LIMIT=10*IC
  80.       KOUNT=0
  81.     1 IF(IR-1)79,79,2
  82. C
  83. C        DROP TRAILING ZERO COEFFICIENTS
  84.     2 IF(C(IR))4,3,4
  85.     3 IR=IR-1
  86.       GOTO 1
  87. C
  88. C           REARRANGEMENT OF GIVEN POLYNOMIAL
  89. C        EXTRACTION OF ZERO ROOTS
  90.     4 O=1.0D0/C(IR)
  91.       IEND=IR-1
  92.       ISTA=1
  93.       NSAV=IR+1
  94.       JBEG=1
  95. C
  96. C        Q(J)=1.
  97. C        Q(J+I)=C(IR-I)/C(IR)
  98. C        Q(IR)=C(J)/C(IR)
  99. C        WHERE J IS THE INDEX OF THE LOWEST NONZERO COEFFICIENT
  100.       DO 9 I=1,IR
  101.       J=NSAV-I
  102.       IF(C(I))7,5,7
  103.     5 GOTO(6,8),JBEG
  104.     6 NSAV=NSAV+1
  105.       Q(ISTA)=0.D0
  106.       E(ISTA)=0.D0
  107.       ISTA=ISTA+1
  108.       GOTO 9
  109.     7 JBEG=2
  110.     8 Q(J)=C(I)*O
  111.       C(I)=Q(J)
  112.     9 CONTINUE
  113. C
  114. C           INITIALIZATION
  115.       ESAV=0.D0
  116.       Q(ISTA)=0.D0
  117.    10 NSAV=IR
  118. C
  119. C        COMPUTATION OF DERIVATIVE
  120.       EXPT=IR-ISTA
  121.       E(ISTA)=EXPT
  122.       DO 11 I=ISTA,IEND
  123.       EXPT=EXPT-1.0D0
  124.       POL(I+1)=EPS*DABS(Q(I+1))+EPS
  125.    11 E(I+1)=Q(I+1)*EXPT
  126. C
  127. C        TEST OF REMAINING DIMENSION
  128.       IF(ISTA-IEND)12,20,60
  129.    12 JEND=IEND-1
  130. C
  131. C        COMPUTATION OF S-FRACTION
  132.       DO 19 I=ISTA,JEND
  133.       IF(I-ISTA)13,16,13
  134.    13 IF(DABS(E(I))-POL(I+1))14,14,16
  135. C
  136. C        THE GIVEN POLYNOMIAL HAS MULTIPLE ROOTS, THE COEFFICIENTS OF
  137. C        THE COMMON FACTOR ARE STORED FROM Q(NSAV) UP TO Q(IR)
  138.    14 NSAV=I
  139.       DO 15 K=I,JEND
  140.       IF(DABS(E(K))-POL(K+1))15,15,80
  141.    15 CONTINUE
  142.       GOTO 21
  143. C
  144. C           EUCLIDEAN ALGORITHM
  145.    16 DO 19 K=I,IEND
  146.       E(K+1)=E(K+1)/E(I)
  147.       Q(K+1)=E(K+1)-Q(K+1)
  148.       IF(K-I)18,17,18
  149. C
  150. C        TEST FOR SMALL DIVISOR
  151.    17 IF(DABS(Q(I+1))-POL(I+1))80,80,19
  152.    18 Q(K+1)=Q(K+1)/Q(I+1)
  153.       POL(K+1)=POL(K+1)/DABS(Q(I+1))
  154.       E(K)=Q(K+1)-E(K)
  155.    19 CONTINUE
  156.    20 Q(IR)=-Q(IR)
  157. C
  158. C           THE DISPLACEMENT EXPT IS SET TO 0 AUTOMATICALLY.
  159. C           E(ISTA)=0.,Q(ISTA+1),...,E(NSAV-1),Q(NSAV),E(NSAV)=0.,
  160. C           FORM A DIAGONAL OF THE QD-ARRAY.
  161. C        INITIALIZATION OF BOUNDARY VALUES
  162.    21 E(ISTA)=0.D0
  163.       NRAN=NSAV-1
  164.    22 E(NRAN+1)=0.D0
  165. C
  166. C           TEST FOR LINEAR OR CONSTANT FACTOR
  167. C        NRAN-ISTA IS DEGREE-1
  168.       IF(NRAN-ISTA)24,23,31
  169. C
  170. C        LINEAR FACTOR
  171.    23 Q(ISTA+1)=Q(ISTA+1)+EXPT
  172.       E(ISTA+1)=0.D0
  173. C
  174. C        TEST FOR UNFACTORED COMMON DIVISOR
  175.    24 E(ISTA)=ESAV
  176.       IF(IR-NSAV)60,60,25
  177. C
  178. C        INITIALIZE QD-ALGORITHM FOR COMMON DIVISOR
  179.    25 ISTA=NSAV
  180.       ESAV=E(ISTA)
  181.       GOTO 10
  182. C
  183. C        COMPUTATION OF ROOT PAIR
  184.    26 P=P+EXPT
  185. C
  186. C        TEST FOR REALITY
  187.       IF(O)27,28,28
  188. C
  189. C        COMPLEX ROOT PAIR
  190.    27 Q(NRAN)=P
  191.       Q(NRAN+1)=P
  192.       E(NRAN)=T
  193.       E(NRAN+1)=-T
  194.       GOTO 29
  195. C
  196. C        REAL ROOT PAIR
  197.    28 Q(NRAN)=P-T
  198.       Q(NRAN+1)=P+T
  199.       E(NRAN)=0.D0
  200. C
  201. C           REDUCTION OF DEGREE BY 2 (DEFLATION)
  202.    29 NRAN=NRAN-2
  203.       GOTO 22
  204. C
  205. C        COMPUTATION OF REAL ROOT
  206.    30 Q(NRAN+1)=EXPT+P
  207. C
  208. C           REDUCTION OF DEGREE BY 1 (DEFLATION)
  209.       NRAN=NRAN-1
  210.       GOTO 22
  211. C
  212. C        START QD-ITERATION
  213.    31 JBEG=ISTA+1
  214.       JEND=NRAN-1
  215.       TEPS=EPS
  216.       TDELT=1.E-2
  217.    32 KOUNT=KOUNT+1
  218.       P=Q(NRAN+1)
  219.       R=ABS(SNGL(E(NRAN)))
  220. C
  221. C           TEST FOR CONVERGENCE
  222.       IF(R-TEPS)30,30,33
  223.    33 S=ABS(SNGL(E(JEND)))
  224. C
  225. C        IS THERE A REAL ROOT NEXT
  226.       IF(S-R)38,38,34
  227. C
  228. C        IS DISPLACEMENT SMALL ENOUGH
  229.    34 IF(R-TDELT)36,35,35
  230.    35 P=0.D0
  231.    36 O=P
  232.       DO 37 J=JBEG,NRAN
  233.       Q(J)=Q(J)+E(J)-E(J-1)-O
  234. C
  235. C           TEST FOR SMALL DIVISOR
  236.       IF(DABS(Q(J))-POL(J))81,81,37
  237.    37 E(J)=Q(J+1)*E(J)/Q(J)
  238.       Q(NRAN+1)=-E(NRAN)+Q(NRAN+1)-O
  239.       GOTO 54
  240. C
  241. C        CALCULATE DISPLACEMENT FOR DOUBLE ROOTS
  242. C           QUADRATIC EQUATION FOR DOUBLE ROOTS
  243. C           X**2-(Q(NRAN)+Q(NRAN+1)+E(NRAN))*X+Q(NRAN)*Q(NRAN+1)=0
  244.    38 P=0.5D0*(Q(NRAN)+E(NRAN)+Q(NRAN+1))
  245.       O=P*P-Q(NRAN)*Q(NRAN+1)
  246.       T=DSQRT(DABS(O))
  247. C
  248. C        TEST FOR CONVERGENCE
  249.       IF(S-TEPS)26,26,39
  250. C
  251. C        ARE THERE COMPLEX ROOTS
  252.    39 IF(O)43,40,40
  253.    40 IF(P)42,41,41
  254.    41 T=-T
  255.    42 P=P+T
  256.       R=S
  257.       GOTO 34
  258. C
  259. C        MODIFICATION FOR COMPLEX ROOTS
  260. C        IS DISPLACEMENT SMALL ENOUGH
  261.    43 IF(S-TDELT)44,35,35
  262. C
  263. C           INITIALIZATION
  264.    44 O=Q(JBEG)+E(JBEG)-P
  265. C
  266. C           TEST FOR SMALL DIVISOR
  267.       IF(DABS(O)-POL(JBEG))81,81,45
  268.    45 T=(T/O)**2
  269.       U=E(JBEG)*Q(JBEG+1)/(O*(1.0D0+T))
  270.       V=O+U
  271. C
  272. C        THREEFOLD LOOP FOR COMPLEX DISPLACEMENT
  273.       KOUNT=KOUNT+2
  274.       DO 53 J=JBEG,NRAN
  275.       O=Q(J+1)+E(J+1)-U-P
  276. C
  277. C           TEST FOR SMALL DIVISOR
  278.       IF(DABS(V)-POL(J))46,46,49
  279.    46 IF(J-NRAN)81,47,81
  280.    47 EXPT=EXPT+P
  281.       IF(ABS(SNGL(E(JEND)))-TOL)48,48,81
  282.    48 P=0.5D0*(V+O-E(JEND))
  283.       O=P*P-(V-U)*(O-U*T-O*W*(1.D0+T)/Q(JEND))
  284.       T=DSQRT(DABS(O))
  285.       GOTO 26
  286. C
  287. C           TEST FOR SMALL DIVISOR
  288.    49 IF(DABS(O)-POL(J+1))46,46,50
  289.    50 W=U*O/V
  290.       T=T*(V/O)**2
  291.       Q(J)=V+W-E(J-1)
  292.       U=0.D0
  293.       IF(J-NRAN)51,52,52
  294.    51 U=Q(J+2)*E(J+1)/(O*(1.D0+T))
  295.    52 V=O+U-W
  296. C
  297. C           TEST FOR SMALL DIVISOR
  298.       IF(DABS(Q(J))-POL(J))81,81,53
  299.    53 E(J)=W*V*(1.0D0+T)/Q(J)
  300.       Q(NRAN+1)=V-E(NRAN)
  301.    54 EXPT=EXPT+P
  302.       TEPS=TEPS*1.1
  303.       TDELT=TDELT*1.1
  304.       IF(KOUNT-LIMIT)32,55,55
  305. C
  306. C        NO CONVERGENCE WITH FEASIBLE TOLERANCE
  307. C           ERROR RETURN IN CASE OF UNSATISFACTORY CONVERGENCE
  308.    55 IER=1
  309. C
  310. C        REARRANGE CALCULATED ROOTS
  311.    56 IEND=NSAV-NRAN-1
  312.       E(ISTA)=ESAV
  313.       IF(IEND)59,59,57
  314.    57 DO 58 I=1,IEND
  315.       J=ISTA+I
  316.       K=NRAN+1+I
  317.       E(J)=E(K)
  318.    58 Q(J)=Q(K)
  319.    59 IR=ISTA+IEND
  320. C
  321. C        NORMAL RETURN
  322.    60 IR=IR-1
  323.       IF(IR)78,78,61
  324. C
  325. C        REARRANGE CALCULATED ROOTS
  326.    61 DO 62 I=1,IR
  327.       Q(I)=Q(I+1)
  328.    62 E(I)=E(I+1)
  329. C
  330. C        CALCULATE COEFFICIENT VECTOR FROM ROOTS
  331.       POL(IR+1)=1.D0
  332.       IEND=IR-1
  333.       JBEG=1
  334.       DO 69 J=1,IR
  335.       ISTA=IR+1-J
  336.       O=0.D0
  337.       P=Q(ISTA)
  338.       T=E(ISTA)
  339.       IF(T)65,63,65
  340. C
  341. C        MULTIPLY WITH LINEAR FACTOR
  342.    63 DO 64 I=ISTA,IR
  343.       POL(I)=O-P*POL(I+1)
  344.    64 O=POL(I+1)
  345.       GOTO 69
  346.    65 GOTO(66,67),JBEG
  347.    66 JBEG=2
  348.       POL(ISTA)=0.D0
  349.       GOTO 69
  350. C
  351. C        MULTIPLY WITH QUADRATIC FACTOR
  352.    67 JBEG=1
  353.       U=P*P+T*T
  354.       P=P+P
  355.       DO 68 I=ISTA,IEND
  356.       POL(I)=O-P*POL(I+1)+U*POL(I+2)
  357.    68 O=POL(I+1)
  358.       POL(IR)=O-P
  359.    69 CONTINUE
  360.       IF(IER)78,70,78
  361. C
  362. C        COMPARISON OF COEFFICIENT VECTORS, IE. TEST OF ACCURACY
  363.    70 P=0.D0
  364.       DO 75 I=1,IR
  365.       IF(C(I))72,71,72
  366.    71 O=DABS(POL(I))
  367.       GOTO 73
  368.    72 O=DABS((POL(I)-C(I))/C(I))
  369.    73 IF(P-O)74,75,75
  370.    74 P=O
  371.    75 CONTINUE
  372.       IF(SNGL(P)-TOL)77,76,76
  373.    76 IER=-1
  374.    77 Q(IR+1)=P
  375.       E(IR+1)=0.D0
  376.    78 RETURN
  377. C
  378. C        ERROR RETURNS
  379. C           ERROR RETURN FOR POLYNOMIALS OF DEGREE LESS THAN 1
  380.    79 IER=2
  381.       IR=0
  382.       RETURN
  383. C
  384. C           ERROR RETURN IF THERE EXISTS NO S-FRACTION
  385.    80 IER=4
  386.       IR=ISTA
  387.       GOTO 60
  388. C
  389. C           ERROR RETURN IN CASE OF INSTABLE QD-ALGORITHM
  390.    81 IER=3
  391.       GOTO 56
  392.       END
  393.