home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / itrpapsm / dapch.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  4.7 KB  |  164 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DAPCH
  5. C
  6. C        PURPOSE
  7. C           SET UP NORMAL EQUATIONS OF LEAST SQUARES FIT IN TERMS OF
  8. C           CHEBYSHEV POLYNOMIALS FOR A GIVEN DISCRETE FUNCTION
  9. C
  10. C        USAGE
  11. C           CALL DAPCH(DATI,N,IP,XD,X0,WORK,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           DATI  - VECTOR OF DIMENSION 3*N (OR DIMENSION 2*N+1)
  15. C                   CONTAINING THE GIVEN ARGUMENTS, FOLLOWED BY THE
  16. C                   FUNCTION VALUES AND N (RESPECTIVELY 1) WEIGHT
  17. C                   VALUES. THE CONTENT OF VECTOR DATI REMAINS
  18. C                   UNCHANGED.
  19. C                   DATI MUST BE OF DOUBLE PRECISION
  20. C           N     - NUMBER OF GIVEN POINTS
  21. C           IP    - DIMENSION OF LEAST SQUARES FIT, I.E. NUMBER OF
  22. C                   CHEBYSHEV POLYNOMIALS USED AS FUNDAMENTAL FUNCTIONS
  23. C                   IP SHOULD NOT EXCEED N
  24. C           XD    - RESULTANT MULTIPLICATIVE CONSTANT FOR LINEAR
  25. C                   TRANSFORMATION OF ARGUMENT RANGE
  26. C                   XD MUST BE DOUBLE PRECISION
  27. C           X0    - RESULTANT ADDITIVE CONSTANT FOR LINEAR
  28. C                   TRANSFORMATION OF ARGUMENT RANGE
  29. C                   X0 MUST BE DOUBLE PRECISION
  30. C           WORK  - WORKING STORAGE OF DIMENSION (IP+1)*(IP+2)/2
  31. C                   ON RETURN WORK CONTAINS THE SYMMETRIC COEFFICIENT
  32. C                   MATRIX OF THE NORMAL EQUATIONS IN COMPRESSED FORM
  33. C                   FOLLOWED IMMEDIATELY BY RIGHT HAND SIDE
  34. C                   AND SQUARE SUM OF FUNCTION VALUES
  35. C                   WORK MUST BE OF DOUBLE PRECISION
  36. C           IER   - RESULTING ERROR PARAMETER
  37. C                   IER =-1 MEANS FORMAL ERRORS IN DIMENSION
  38. C                   IER = 0 MEANS NO ERRORS
  39. C                   IER = 1 MEANS COINCIDING ARGUMENTS
  40. C
  41. C        REMARKS
  42. C           NO WEIGHTS ARE USED IF THE VALUE OF DATI(2*N+1) IS
  43. C           NOT POSITIVE.
  44. C           EXECUTION OF SUBROUTINE DAPCH IS A PREPARATORY STEP FOR
  45. C           CALCULATION OF LEAST SQUARES FITS IN CHEBYSHEV POLYNOMIALS
  46. C           IT SHOULD BE FOLLOWED BY EXECUTION OF SUBROUTINE DAPFS
  47. C
  48. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  49. C           NONE
  50. C
  51. C        METHOD
  52. C           THE LEAST SQUARE FIT IS DETERMINED USING CHEBYSHEV
  53. C           POLYNOMIALS AS FUNDAMENTAL FUNCTION SYSTEM.
  54. C           THE METHOD IS DISCUSSED IN THE ARTICLE
  55. C           A.T.BERZTISS, LEAST SQUARES FITTING TO IRREGULARLY SPACED
  56. C           DATA, SIAM REVIEW, VOL.6, ISS.3, 1964, PP. 203-227.
  57. C
  58. C     ..................................................................
  59. C
  60.       SUBROUTINE DAPCH(DATI,N,IP,XD,X0,WORK,IER)
  61. C
  62. C
  63. C       DIMENSIONED DUMMY VARIABLES
  64.       DIMENSION DATI(1),WORK(1)
  65.       DOUBLE PRECISION DATI,WORK,XD,X0,XA,XE,XM,DF,T,SUM
  66. C
  67. C        CHECK FOR FORMAL ERRORS IN SPECIFIED DIMENSIONS
  68.       IF(N-1)19,20,1
  69.     1 IF(IP)19,19,2
  70. C
  71. C        SEARCH SMALLEST AND LARGEST ARGUMENT
  72.     2 IF(IP-N)3,3,19
  73.     3 XA=DATI(1)
  74.       X0=XA
  75.       XE=0.D0
  76.       DO 7 I=1,N
  77.       XM=DATI(I)
  78.       IF(XA-XM)5,5,4
  79.     4 XA=XM
  80.     5 IF(X0-XM)6,7,7
  81.     6 X0=XM
  82.     7 CONTINUE
  83. C
  84. C        INITIALIZE CALCULATION OF NORMAL EQUATIONS
  85.       XD=X0-XA
  86.       M=(IP*(IP+1))/2
  87.       IEND=M+IP+1
  88.       MT2=IP+IP
  89.       MT2M=MT2-1
  90. C
  91. C        SET WORKING STORAGE AND RIGHT HAND SIDE TO ZERO
  92.       DO 8 I=1,IP
  93.       J=MT2-I
  94.       WORK(J)=0.D0
  95.       WORK(I)=0.D0
  96.       K=M+I
  97.     8 WORK(K)=0.D0
  98. C
  99. C        CHECK FOR DEGENERATE ARGUMENT RANGE
  100.       IF(XD)20,20,9
  101. C
  102. C        CALCULATE CONSTANTS FOR REDUCTION OF ARGUMENTS
  103.     9 X0=-(X0+XA)/XD
  104.       XD=2.D0/XD
  105.       SUM=0.D0
  106. C
  107. C        START GREAT LOOP OVER ALL GIVEN POINTS
  108.       DO 15 I=1,N
  109.       T=DATI(I)*XD+X0
  110.       J=I+N
  111.       DF=DATI(J)
  112. C
  113. C        CALCULATE AND STORE VALUES OF CHEBYSHEV POLYNOMIALS
  114. C        FOR ARGUMENT T
  115.       XA=1.D0
  116.       XM=T
  117.       IF(DATI(2*N+1))11,11,10
  118.    10 J=J+N
  119.       XA=DATI(J)
  120.       XM=T*XA
  121.    11 T=T+T
  122.       SUM=SUM+DF*DF*XA
  123.       DF=DF+DF
  124.       J=1
  125.    12 K=M+J
  126.       WORK(K)=WORK(K)+DF*XA
  127.    13 WORK(J)=WORK(J)+XA
  128.       IF(J-MT2M)14,15,15
  129.    14 J=J+1
  130.       XE=T*XM-XA
  131.       XA=XM
  132.       XM=XE
  133.       IF(J-IP)12,12,13
  134.    15 CONTINUE
  135.       WORK(IEND)=SUM+SUM
  136. C
  137. C        CALCULATE MATRIX OF NORMAL EQUATIONS
  138.       LL=M
  139.       KK=MT2M
  140.       JJ=1
  141.       K=KK
  142.       DO 18 J=1,M
  143.       WORK(LL)=WORK(K)+WORK(JJ)
  144.       LL=LL-1
  145.       IF(K-JJ)16,16,17
  146.    16 KK=KK-2
  147.       K=KK
  148.       JJ=1
  149.       GOTO 18
  150.    17 JJ=JJ+1
  151.       K=K-1
  152.    18 CONTINUE
  153.       IER=0
  154.       RETURN
  155. C
  156. C        ERROR RETURN IN CASE OF FORMAL ERRORS
  157.    19 IER=-1
  158.       RETURN
  159. C
  160. C        ERROR RETURN IN CASE OF COINCIDING ARGUMENTS
  161.    20 IER=1
  162.       RETURN
  163.       END
  164.