C
C     ..................................................................
C
C        SUBROUTINE DAHI
C
C        PURPOSE
C           TO INTERPOLATE FUNCTION VALUE Y FOR A GIVEN ARGUMENT VALUE
C           X USING A GIVEN TABLE (ARG,VAL) OF ARGUMENT, FUNCTION, AND
C           DERIVATIVE VALUES.
C
C        USAGE
C           CALL DAHI (X,ARG,VAL,Y,NDIM,EPS,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X      - DOUBLE PRECISION ARGUMENT VALUE SPECIFIED BY INPUT.
C           ARG    - DOUBLE PRECISION INPUT VECTOR (DIMENSION NDIM) OF
C                    ARGUMENT VALUES OF THE TABLE (NOT DESTROYED).
C           VAL    - DOUBLE PRECISION INPUT VECTOR (DIMENSION 2*NDIM) OF
C                    FUNCTION AND DERIVATIVE VALUES OF THE TABLE (DES-
C                    TROYED). FUNCTION AND DERIVATIVE VALUES MUST BE
C                    STORED IN PAIRS, THAT MEANS BEGINNING WITH FUNCTION
C                    VALUE AT POINT ARG(1) EVERY FUNCTION VALUE MUST BE
C                    FOLLOWED BY THE DERIVATIVE VALUE AT THE SAME POINT.
C           Y      - RESULTING INTERPOLATED DOUBLE PRECISION FUNCTION
C                    VALUE.
C           NDIM   - AN INPUT VALUE WHICH SPECIFIES THE NUMBER OF
C                    POINTS IN TABLE (ARG,VAL).
C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED AS
C                    UPPER BOUND FOR THE ABSOLUTE ERROR.
C           IER    - A RESULTING ERROR PARAMETER.
C
C        REMARKS
C           (1) TABLE (ARG,VAL) SHOULD REPRESENT A SINGLE-VALUED
C               FUNCTION AND SHOULD BE STORED IN SUCH A WAY, THAT THE
C               DISTANCES ABS(ARG(I)-X) INCREASE WITH INCREASING
C               SUBSCRIPT I. TO GENERATE THIS ORDER IN TABLE (ARG,VAL),
C               SUBROUTINES DATSG, DATSM OR DATSE COULD BE USED IN A
C               PREVIOUS STAGE.
C           (2) NO ACTION BESIDES ERROR MESSAGE IN CASE NDIM LESS
C               THAN 1.
C           (3) INTERPOLATION IS TERMINATED EITHER IF THE DIFFERENCE
C               BETWEEN TWO SUCCESSIVE INTERPOLATED VALUES IS
C               ABSOLUTELY LESS THAN TOLERANCE EPS, OR IF THE ABSOLUTE
C               VALUE OF THIS DIFFERENCE STOPS DIMINISHING, OR AFTER
C               (2*NDIM-2) STEPS. FURTHER IT IS TERMINATED IF THE
C               PROCEDURE DISCOVERS TWO ARGUMENT VALUES IN VECTOR ARG
C               WHICH ARE IDENTICAL. DEPENDENT ON THESE FOUR CASES,
C               ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
C                IER=0 - IT WAS POSSIBLE TO REACH THE REQUIRED
C                        ACCURACY (NO ERROR).
C                IER=1 - IT WAS IMPOSSIBLE TO REACH THE REQUIRED
C                        ACCURACY BECAUSE OF ROUNDING ERRORS.
C                IER=2 - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE
C                        NDIM IS LESS THAN 2, OR THE REQUIRED ACCURACY
C                        COULD NOT BE REACHED BY MEANS OF THE GIVEN
C                        TABLE. NDIM SHOULD BE INCREASED.
C                IER=3 - THE PROCEDURE DISCOVERED TWO ARGUMENT VALUES
C                        IN VECTOR ARG WHICH ARE IDENTICAL.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           INTERPOLATION IS DONE BY MEANS OF AITKENS SCHEME OF
C           HERMITE INTERPOLATION. ON RETURN Y CONTAINS AN INTERPOLATED
C           FUNCTION VALUE AT POINT X, WHICH IS IN THE SENSE OF REMARK
C           (3) OPTIMAL WITH RESPECT TO GIVEN TABLE. FOR REFERENCE, SEE
C           F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
C           MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.314-317, AND
C           GERSHINSKY/LEVINE, AITKEN-HERMITE INTERPOLATION,
C           JACM, VOL.11, ISS.3 (1964), PP.352-356.
C
C     ..................................................................
C
      SUBROUTINE DAHI(X,ARG,VAL,Y,NDIM,EPS,IER)
C
C
      DIMENSION ARG(1),VAL(1)
      DOUBLE PRECISION ARG,VAL,X,Y,H,H1,H2
      IER=2
      H2=X-ARG(1)
      IF(NDIM-1)2,1,3
    1 Y=VAL(1)+VAL(2)*H2
    2 RETURN
C
C     VECTOR ARG HAS MORE THAN 1 ELEMENT.
C     THE FIRST STEP PREPARES VECTOR VAL SUCH THAT AITKEN SCHEME CAN BE
C     USED.
    3 I=1
      DO 5 J=2,NDIM
      H1=H2
      H2=X-ARG(J)
      Y=VAL(I)
      VAL(I)=Y+VAL(I+1)*H1
      H=H1-H2
      IF(H)4,13,4
    4 VAL(I+1)=Y+(VAL(I+2)-Y)*H1/H
    5 I=I+2
      VAL(I)=VAL(I)+VAL(I+1)*H2
C     END OF FIRST STEP
C
C     PREPARE AITKEN SCHEME
      DELT2=0.
      IEND=I-1
C
C     START AITKEN-LOOP
      DO 9 I=1,IEND
      DELT1=DELT2
      Y=VAL(1)
      M=(I+3)/2
      H1=ARG(M)
      DO 6 J=1,I
      K=I+1-J
      L=(K+1)/2
      H=ARG(L)-H1
      IF(H)6,14,6
    6 VAL(K)=(VAL(K)*(X-H1)-VAL(K+1)*(X-ARG(L)))/H
      DELT2=DABS(Y-VAL(1))
      IF(DELT2-EPS)11,11,7
    7 IF(I-8)9,8,8
    8 IF(DELT2-DELT1)9,12,12
    9 CONTINUE
C     END OF AITKEN-LOOP
C
   10 Y=VAL(1)
      RETURN
C
C     THERE IS SUFFICIENT ACCURACY WITHIN 2*NDIM-2 ITERATION STEPS
   11 IER=0
      GOTO 10
C
C     TEST VALUE DELT2 STARTS OSCILLATING
   12 IER=1
      RETURN
C
C     THERE ARE TWO IDENTICAL ARGUMENT VALUES IN VECTOR ARG
   13 Y=VAL(1)
   14 IER=3
      RETURN
      END