C
C     ..................................................................
C
C        SUBROUTINE QATR
C
C        PURPOSE
C           TO COMPUTE AN APPROXIMATION FOR INTEGRAL(FCT(X), SUMMED
C           OVER X FROM XL TO XU).
C
C        USAGE
C           CALL QATR (XL,XU,EPS,NDIM,FCT,Y,IER,AUX)
C           PARAMETER FCT REQUIRES AN EXTERNAL STATEMENT.
C
C        DESCRIPTION OF PARAMETERS
C           XL     - THE LOWER BOUND OF THE INTERVAL.
C           XU     - THE UPPER BOUND OF THE INTERVAL.
C           EPS    - THE UPPER BOUND OF THE ABSOLUTE ERROR.
C           NDIM   - THE DIMENSION OF THE AUXILIARY STORAGE ARRAY AUX.
C                    NDIM-1 IS THE MAXIMAL NUMBER OF BISECTIONS OF
C                    THE INTERVAL (XL,XU).
C           FCT    - THE NAME OF THE EXTERNAL FUNCTION SUBPROGRAM USED.
C           Y      - THE RESULTING APPROXIMATION FOR THE INTEGRAL VALUE.
C           IER    - A RESULTING ERROR PARAMETER.
C           AUX    - AN AUXILIARY STORAGE ARRAY WITH DIMENSION NDIM.
C
C        REMARKS
C           ERROR PARAMETER IER IS CODED IN THE FOLLOWING FORM
C           IER=0  - IT WAS POSSIBLE TO REACH THE REQUIRED ACCURACY.
C                    NO ERROR.
C           IER=1  - IT IS IMPOSSIBLE TO REACH THE REQUIRED ACCURACY
C                    BECAUSE OF ROUNDING ERRORS.
C           IER=2  - IT WAS IMPOSSIBLE TO CHECK ACCURACY BECAUSE NDIM
C                    IS LESS THAN 5, OR THE REQUIRED ACCURACY COULD NOT
C                    BE REACHED WITHIN NDIM-1 STEPS. NDIM SHOULD BE
C                    INCREASED.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           THE EXTERNAL FUNCTION SUBPROGRAM FCT(X) MUST BE CODED BY
C           THE USER. ITS ARGUMENT X SHOULD NOT BE DESTROYED.
C
C        METHOD
C           EVALUATION OF Y IS DONE BY MEANS OF TRAPEZOIDAL RULE IN
C           CONNECTION WITH ROMBERGS PRINCIPLE. ON RETURN Y CONTAINS
C           THE BEST POSSIBLE APPROXIMATION OF THE INTEGRAL VALUE AND
C           VECTOR AUX THE UPWARD DIAGONAL OF ROMBERG SCHEME.
C           COMPONENTS AUX(I) (I=1,2,...,IEND, WITH IEND LESS THAN OR
C           EQUAL TO NDIM) BECOME APPROXIMATIONS TO INTEGRAL VALUE WITH
C           DECREASING ACCURACY BY MULTIPLICATION WITH (XU-XL).
C           FOR REFERENCE, SEE
C           (1) FILIPPI, DAS VERFAHREN VON ROMBERG-STIEFEL-BAUER ALS
C               SPEZIALFALL DES ALLGEMEINEN PRINZIPS VON RICHARDSON,
C               MATHEMATIK-TECHNIK-WIRTSCHAFT, VOL.11, ISS.2 (1964),
C               PP.49-54.
C           (2) BAUER, ALGORITHM 60, CACM, VOL.4, ISS.6 (1961), PP.255.
C
C     ..................................................................
C
      SUBROUTINE QATR(XL,XU,EPS,NDIM,FCT,Y,IER,AUX)
C
C
      DIMENSION AUX(1)
C
C     PREPARATIONS OF ROMBERG-LOOP
      AUX(1)=.5*(FCT(XL)+FCT(XU))
      H=XU-XL
      IF(NDIM-1)8,8,1
    1 IF(H)2,10,2
C
C     NDIM IS GREATER THAN 1 AND H IS NOT EQUAL TO 0.
    2 HH=H
      E=EPS/ABS(H)
      DELT2=0.
      P=1.
      JJ=1
      DO 7 I=2,NDIM
      Y=AUX(1)
      DELT1=DELT2
      HD=HH
      HH=.5*HH
      P=.5*P
      X=XL+HH
      SM=0.
      DO 3 J=1,JJ
      SM=SM+FCT(X)
    3 X=X+HD
      AUX(I)=.5*AUX(I-1)+P*SM
C     A NEW APPROXIMATION OF INTEGRAL VALUE IS COMPUTED BY MEANS OF
C     TRAPEZOIDAL RULE.
C
C     START OF ROMBERGS EXTRAPOLATION METHOD.
      Q=1.
      JI=I-1
      DO 4 J=1,JI
      II=I-J
      Q=Q+Q
      Q=Q+Q
    4 AUX(II)=AUX(II+1)+(AUX(II+1)-AUX(II))/(Q-1.)
C     END OF ROMBERG-STEP
C
      DELT2=ABS(Y-AUX(1))
      IF(I-5)7,5,5
    5 IF(DELT2-E)10,10,6
    6 IF(DELT2-DELT1)7,11,11
    7 JJ=JJ+JJ
    8 IER=2
    9 Y=H*AUX(1)
      RETURN
   10 IER=0
      GO TO 9
   11 IER=1
      Y=H*Y
      RETURN
      END
