home *** CD-ROM | disk | FTP | other *** search
/ Frostbyte's 1980s DOS Shareware Collection / floppyshareware.zip / floppyshareware / DOOG / PCSSP1.ZIP / ITRPAPSM.ZIP / DARAT.FOR < prev    next >
Text File  |  1985-11-29  |  8KB  |  259 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DARAT
  5. C
  6. C        PURPOSE
  7. C           CALCULATE BEST RATIONAL APPROXIMATION OF A DISCRETE
  8. C           FUNCTION IN THE LEAST SQUARES SENSE
  9. C
  10. C        USAGE
  11. C           CALL DARAT(DATI,N,WORK,P,IP,IQ,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           DATI  - TWODIMENSIONAL ARRAY WITH 3 COLUMNS AND N ROWS
  15. C                   THE FIRST COLUMN MUST CONTAIN THE GIVEN ARGUMENTS,
  16. C                   THE SECOND COLUMN THE GIVEN FUNCTION VALUES AND
  17. C                   THE THIRD COLUMN THE GIVEN WEIGHTS IF ANY.
  18. C                   IF NO WEIGHTS ARE TO BE USED THEN THE THIRD
  19. C                   COLUMN MAY BE DROPPED , EXCEPT THE FIRST ELEMENT
  20. C                   WHICH MUST CONTAIN A NONPOSITIVE VALUE
  21. C                   DATI MUST BE OF DOUBLE PRECISION
  22. C           N     - NUMBER OF NODES OF THE GIVEN DISCRETE FUNCTION
  23. C           WORK  - WORKING STORAGE WHICH IS OF DIMENSION
  24. C                   (IP+IQ)*(IP+IQ+1)+4*N+1 AT LEAST.
  25. C                   ON RETURN THE VALUES OF THE NUMERATOR ARE CONTAINED
  26. C                   IN WORK(N+1) UP TO WORK(2*N), WHILE THE VALUES OF
  27. C                   THE DENOMINATOR ARE STORED IN WORK(2*N+1) UP TO
  28. C                   WORK(3*N)
  29. C                   WORK MUST BE OF DOUBLE PRECISION
  30. C           P     - RESULTANT COEFFICIENT VECTOR OF DENOMINATOR AND
  31. C                   NUMERATOR. THE DENOMINATOR IS STORED IN FIRST IQ
  32. C                   LOCATIONS, THE NUMERATOR IN THE FOLLOWING IP
  33. C                   LOCATIONS.
  34. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH.
  35. C                   P MUST BE OF DOUBLE PRECISION
  36. C           IP    - DIMENSION OF THE NUMERATOR   (INPUT VALUE)
  37. C           IQ    - DIMENSION OF THE DENOMINATOR (INPUT VALUE)
  38. C           IER   - RESULTANT ERROR PARAMETER
  39. C                   IER =-1 MEANS FORMAL ERRORS
  40. C                   IER = 0 MEANS NO ERRORS
  41. C                   IER = 1,2 MEANS POOR CONVERGENCE OF ITERATION
  42. C                   IER IS ALSO USED AS INPUT VALUE
  43. C                   A NONZERO INPUT VALUE INDICATES AVAILABILITY OF AN
  44. C                   INITIAL APPROXIMATION STORED IN P
  45. C
  46. C        REMARKS
  47. C           THE COEFFICIENT VECTORS OF THE DENOMINATOR AND NUMERATOR
  48. C           OF THE RATIONAL APPROXIMATION ARE BOTH STORED IN P
  49. C           STARTING WITH LOW POWERS (DENOMINATOR FIRST).
  50. C           IP+IQ MUST NOT EXCEED N, ALL THREE VALUES MUST BE POSITIVE.
  51. C           SINCE CHEBYSHEV POLYNOMIALS ARE USED AS FUNDAMENTAL
  52. C           FUNCTIONS, THE ARGUMENTS SHOULD BE REDUCED TO THE INTERVAL
  53. C           (-1,1). THIS CAN ALWAYS BE ACCOMPLISHED BY MEANS OF A LINEAR
  54. C           TRANSFORMATION OF THE ORIGINALLY GIVEN ARGUMENTS.
  55. C           IF A FIT IN OTHER FUNCTIONS IS REQUIRED, DCNP AND DCNPS MUST
  56. C           BE REPLACED BY SUBROUTINES WHICH ARE OF ANALOGOUS DESIGN.
  57. C
  58. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  59. C           DAPLL, DAPFS, DFRAT, DCNPS, DCNP
  60. C           DCNP IS REQUIRED WITHIN DFRAT
  61. C
  62. C        METHOD
  63. C           THE ITERATIVE SCHEME USED FOR CALCULATION OF THE
  64. C           APPROXIMATION IS REPEATED SOLUTION OF THE NORMAL EQUATIONS
  65. C           WHICH ARE OBTAINED BY LINEARIZATION.
  66. C           A REFINED TECHNIQUE OF THIS LINEAR LEAST SQUARES APPROACH
  67. C           IS USED WHICH GUARANTEES THAT THE DENOMINATOR IS FREE OF
  68. C           ZEROES WITHIN THE APPROXIMATION INTERVAL.
  69. C           FOR REFERENCE SEE
  70. C           D.BRAESS, UEBER DAEMPFUNG BEI MINIMALISIERUNGSVERFAHREN,
  71. C           COMPUTING(1966), VOL.1, ED.3, PP.264-272.
  72. C           D.W.MARQUARDT, AN ALGORITHM FOR LEAST-SQUARES ESTIMATION
  73. C           OF NONLINEAR PARAMETERS,
  74. C           JSIAM(1963), VOL.11, ED.2, PP.431-441.
  75. C
  76. C     ..................................................................
  77. C
  78.       SUBROUTINE DARAT(DATI,N,WORK,P,IP,IQ,IER)
  79. C
  80. C
  81.       EXTERNAL DFRAT
  82. C
  83. C        DIMENSIONED LOCAL VARIABLE
  84.       DIMENSION IERV(3)
  85. C
  86. C        DIMENSIONED DUMMY VARIABLES
  87.       DIMENSION DATI(1),WORK(1),P(1)
  88.       DOUBLE PRECISION DATI,WORK,P,T,OSUM,DIAG,RELAX,SUM,SSOE,SAVE
  89. C
  90. C        INITIALIZE TESTVALUES
  91.       LIMIT=20
  92.       ETA=1.E-29
  93.       EPS=1.E-14
  94. C
  95. C        CHECK FOR FORMAL ERRORS
  96.       IF(N)4,4,1
  97.     1 IF(IP)4,4,2
  98.     2 IF(IQ)4,4,3
  99.     3 IPQ=IP+IQ
  100.       IF(N-IPQ)4,5,5
  101. C
  102. C        ERROR RETURN IN CASE OF FORMAL ERRORS
  103.     4 IER=-1
  104.       RETURN
  105. C
  106. C        INITIALIZE ITERATION PROCESS
  107.     5 KOUNT=0
  108.       IERV(2)=IP
  109.       IERV(3)=IQ
  110.       NDP=N+N+1
  111.       NNE=NDP+NDP
  112.       IX=IPQ-1
  113.       IQP1=IQ+1
  114.       IRHS=NNE+IPQ*IX/2
  115.       IEND=IRHS+IX
  116. C
  117. C        TEST FOR AVAILABILITY OF AN INITIAL APPROXIMATION
  118.       IF(IER)8,6,8
  119. C
  120. C        INITIALIZE NUMERATOR AND DENOMINATOR
  121.     6 DO 7 I=2,IPQ
  122.     7 P(I)=0.D0
  123.       P(1)=1.D0
  124. C
  125. C        CALCULATE VALUES OF NUMERATOR AND DENOMINATOR FOR INITIAL
  126. C        APPROXIMATION
  127.     8 DO 9 J=1,N
  128.       T=DATI(J)
  129.       I=J+N
  130.       CALL DCNPS(WORK(I),T,P(IQP1),IP)
  131.       K=I+N
  132.     9 CALL DCNPS(WORK(K),T,P,IQ)
  133. C
  134. C        SET UP NORMAL EQUATIONS (MAIN LOOP OF ITERATION)
  135.    10 CALL DAPLL(DFRAT,N,IX,WORK,WORK(IEND+1),DATI,IERV)
  136. C
  137. C        CHECK FOR ZERO DENOMINATOR
  138.       IF(IERV(1))4,11,4
  139.    11 INCR=0
  140.       RELAX=2.D0
  141. C
  142. C        RESTORE MATRIX IN WORKING STORAGE
  143.    12 J=IEND
  144.       DO 13 I=NNE,IEND
  145.       J=J+1
  146.    13 WORK(I)=WORK(J)
  147.       IF(KOUNT)14,14,15
  148. C
  149. C        SAVE SQUARE SUM OF ERRORS
  150.    14 OSUM=WORK(IEND)
  151.       DIAG=OSUM*EPS
  152.       K=IQ
  153. C
  154. C        ADD CONSTANT TO DIAGONAL
  155.       IF(WORK(NNE))17,17,19
  156.    15 IF(INCR)19,19,16
  157.    16 K=IPQ
  158.    17 J=NNE-1
  159.       DO 18 I=1,K
  160.       WORK(J)=WORK(J)+DIAG
  161.    18 J=J+I
  162. C
  163. C        SOLVE NORMAL EQUATIONS
  164.    19 CALL DAPFS(WORK(NNE),IX,IRES,1,EPS,ETA,IER)
  165. C
  166. C        CHECK FOR FAILURE OF EQUATION SOLVER
  167.       IF(IRES)4,4,20
  168. C
  169. C        TEST FOR DEFECTIVE NORMALEQUATIONS
  170.    20 IF(IRES-IX)21,24,24
  171.    21 IF(INCR)22,22,23
  172.    22 DIAG=DIAG*0.125D0
  173.    23 DIAG=DIAG+DIAG
  174.       INCR=INCR+1
  175. C
  176. C        START WITH OVER RELAXATION
  177.       RELAX=8.D0
  178.       IF(INCR-LIMIT)12,45,45
  179. C
  180. C        CALCULATE VALUES OF CHANGE OF NUMERATOR AND DENOMINATOR
  181.    24 L=NDP
  182.       J=NNE+IRES*(IRES-1)/2-1
  183.       K=J+IQ
  184.       WORK(J)=0.D0
  185.       IRQ=IQ
  186.       IRP=IRES-IQ+1
  187.       IF(IRP)25,26,26
  188.    25 IRQ=IRES+1
  189.    26 DO 29 I=1,N
  190.       T=DATI(I)
  191.       WORK(I)=0.D0
  192.       CALL DCNPS(WORK(I),T,WORK(K),IRP)
  193.       M=L+N
  194.       CALL DCNPS(WORK(M),T,WORK(J),IRQ)
  195.       IF(WORK(M)*WORK(L))27,29,29
  196.    27 SUM=WORK(L)/WORK(M)
  197.       IF(RELAX+SUM)29,29,28
  198.    28 RELAX=-SUM
  199.    29 L=L+1
  200. C
  201. C        MODIFY RELAXATION FACTOR IF NECESSARY
  202.       SSOE=OSUM
  203.       ITER=LIMIT
  204.    30 SUM=0.D0
  205.       RELAX=RELAX*0.5D0
  206.       DO 32 I=1,N
  207.       M=I+N
  208.       K=M+N
  209.       L=K+N
  210.       SAVE=DATI(M)-(WORK(M)+RELAX*WORK(I))/(WORK(K)+RELAX*WORK(L))
  211.       SAVE=SAVE*SAVE
  212.       IF(DATI(NDP))32,32,31
  213.    31 SAVE=SAVE*DATI(K)
  214.    32 SUM=SUM+SAVE
  215.       IF(ITER)45,33,33
  216.    33 ITER=ITER-1
  217.       IF(SUM-OSUM)34,37,35
  218.    34 OSUM=SUM
  219.       GOTO 30
  220. C
  221. C        TEST FOR IMPROVEMENT
  222.    35 IF(OSUM-SSOE)36,30,30
  223.    36 RELAX=RELAX+RELAX
  224.    37 T=0.
  225.       SAVE=0.D0
  226.       K=IRES+1
  227.       DO 38 I=2,K
  228.       J=J+1
  229.       T=T+DABS(P(I))
  230.       P(I)=P(I)+RELAX*WORK(J)
  231.    38 SAVE=SAVE+DABS(P(I))
  232. C
  233. C        UPDATE CURRENT VALUES OF NUMERATOR AND DENOMINATOR
  234.       DO 39 I=1,N
  235.       J=I+N
  236.       K=J+N
  237.       L=K+N
  238.       WORK(J)=WORK(J)+RELAX*WORK(I)
  239.    39 WORK(K)=WORK(K)+RELAX*WORK(L)
  240. C
  241. C        TEST FOR CONVERGENCE
  242.       IF(INCR)40,40,42
  243.    40 IF(SSOE-OSUM-RELAX*OSUM*DBLE(EPS))46,46,41
  244.    41 IF(DABS(T-SAVE)-RELAX*SAVE*DBLE(EPS))46,46,42
  245.    42 IF(OSUM-SAVE*DBLE(ETA))46,46,43
  246.    43 KOUNT=KOUNT+1
  247.       IF(KOUNT-LIMIT)10,44,44
  248. C
  249. C        ERROR RETURN IN CASE OF POOR CONVERGENCE
  250.    44 IER=2
  251.       RETURN
  252.    45 IER=1
  253.       RETURN
  254. C
  255. C        NORMAL RETURN
  256.    46 IER=0
  257.       RETURN
  258.       END
  259.