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