home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / fortran / library / ssp / itrpapsm / dahi.for < prev    next >
Encoding:
Text File  |  1985-11-29  |  5.1 KB  |  141 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DAHI
  5. C
  6. C        PURPOSE
  7. C           TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE
  8. C           X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT, FUNCTION, AND
  9. C           DERIVATIVE VALUES.
  10. C
  11. C        USAGE
  12. C           CALL DAHI (X,ARG,VAL,Y,NDIM,EPS,IER)
  13. C
  14. C        DESCRIPTION OF PARAMETERS
  15. C           X      - DOUBLE PRECISION ARGUMENT VALUE SPECIFIED BY INPUT.
  16. C           ARG    - DOUBLE PRECISION INPUT VECTOR (DIMENSION NDIM) OF
  17. C                    ARGUMENT VALUES OF THE TABLE (NOT DESTROYED).
  18. C           VAL    - DOUBLE PRECISION INPUT VECTOR (DIMENSION 2*NDIM) OF
  19. C                    FUNCTION AND DERIVATIVE VALUES OF THE TABLE (DES-
  20. C                    TROYED). FUNCTION AND DERIVATIVE VALUES MUST BE
  21. C                    STORED IN PAIRS, THAT MEANS BEGINNING WITH FUNCTION
  22. C                    VALUE AT POINT ARG(1) EVERY FUNCTION VALUE MUST BE
  23. C                    FOLLOWED BY THE DERIVATIVE VALUE AT THE SAME POINT.
  24. C           Y      - RESULTING INTERPOLATED DOUBLE PRECISION FUNCTION
  25. C                    VALUE.
  26. C           NDIM   - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
  27. C                    POINTS IN TABLE (ARG,VAL).
  28. C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS
  29. C                    UPPER BOUND FOR THE ABSOLUTE ERROR.
  30. C           IER    - A RESULTING ERROR PARAMETER.
  31. C
  32. C        REMARKS
  33. C           (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
  34. C               FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
  35. C               DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
  36. C               SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
  37. C               SUBROUTINES DATSG, DATSM OR DATSE COULD BE USED IN A
  38. C               PREVIOUS STAGE.
  39. C           (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
  40. C               THAN 1.
  41. C           (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
  42. C               BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
  43. C               ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
  44. C               VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
  45. C               (2*NDIM-2) STEPS. FURTHER IT IS TERMINATED IF THE
  46. C               PROCEDURE DISCOVERS TWO ARGUMENT VALUES IN VECTOR ARG
  47. C               WHICH ARE IDENTICAL. DEPENDENT ON THESE FOUR CASES,
  48. C               ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
  49. C                IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
  50. C                        ACCURACY (NO ERROR).
  51. C                IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
  52. C                        ACCURACY BECAUSE OF ROUNDING ERRORS.
  53. C                IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
  54. C                        NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
  55. C                        COULD NOT BE REACHED BY MEANS OF THE GIVEN
  56. C                        TABLE. NDIM SHOULD BE INCREASED.
  57. C                IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
  58. C                        IN VECTOR ARG WHICH ARE IDENTICAL.
  59. C
  60. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  61. C           NONE
  62. C
  63. C        METHOD
  64. C           INTERPOLATION IS DONE BY MEANS OF AITKENS SCHEME OF
  65. C           HERMITE INTERPOLATION. ON RETURN Y CONTAINS AN INTERPOLATED
  66. C           FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
  67. C           (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
  68. C           F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
  69. C           MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.314-317, AND
  70. C           GERSHINSKY/LEVINE, AITKEN-HERMITE INTERPOLATION,
  71. C           JACM, VOL.11, ISS.3 (1964), PP.352-356.
  72. C
  73. C     ..................................................................
  74. C
  75.       SUBROUTINE DAHI(X,ARG,VAL,Y,NDIM,EPS,IER)
  76. C
  77. C
  78.       DIMENSION ARG(1),VAL(1)
  79.       DOUBLE PRECISION ARG,VAL,X,Y,H,H1,H2
  80.       IER=2
  81.       H2=X-ARG(1)
  82.       IF(NDIM-1)2,1,3
  83.     1 Y=VAL(1)+VAL(2)*H2
  84.     2 RETURN
  85. C
  86. C     VECTOR ARG HAS MORE THAN 1 ELEMENT.
  87. C     THE FIRST STEP PREPARES VECTOR VAL SUCH THAT AITKEN SCHEME CAN BE
  88. C     USED.
  89.     3 I=1
  90.       DO 5 J=2,NDIM
  91.       H1=H2
  92.       H2=X-ARG(J)
  93.       Y=VAL(I)
  94.       VAL(I)=Y+VAL(I+1)*H1
  95.       H=H1-H2
  96.       IF(H)4,13,4
  97.     4 VAL(I+1)=Y+(VAL(I+2)-Y)*H1/H
  98.     5 I=I+2
  99.       VAL(I)=VAL(I)+VAL(I+1)*H2
  100. C     END OF FIRST STEP
  101. C
  102. C     PREPARE AITKEN SCHEME
  103.       DELT2=0.
  104.       IEND=I-1
  105. C
  106. C     START AITKEN-LOOP
  107.       DO 9 I=1,IEND
  108.       DELT1=DELT2
  109.       Y=VAL(1)
  110.       M=(I+3)/2
  111.       H1=ARG(M)
  112.       DO 6 J=1,I
  113.       K=I+1-J
  114.       L=(K+1)/2
  115.       H=ARG(L)-H1
  116.       IF(H)6,14,6
  117.     6 VAL(K)=(VAL(K)*(X-H1)-VAL(K+1)*(X-ARG(L)))/H
  118.       DELT2=DABS(Y-VAL(1))
  119.       IF(DELT2-EPS)11,11,7
  120.     7 IF(I-8)9,8,8
  121.     8 IF(DELT2-DELT1)9,12,12
  122.     9 CONTINUE
  123. C     END OF AITKEN-LOOP
  124. C
  125.    10 Y=VAL(1)
  126.       RETURN
  127. C
  128. C     THERE IS SUFFICIENT ACCURACY WITHIN 2*NDIM-2 ITERATION STEPS
  129.    11 IER=0
  130.       GOTO 10
  131. C
  132. C     TEST VALUE DELT2 STARTS OSCILLATING
  133.    12 IER=1
  134.       RETURN
  135. C
  136. C     THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
  137.    13 Y=VAL(1)
  138.    14 IER=3
  139.       RETURN
  140.       END
  141.