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