C
C     ..................................................................
C
C        SUBROUTINE DQSF
C
C        PURPOSE
C           TO COMPUTE THE VECTOR OF INTEGRAL VALUES FOR A GIVEN
C           EQUIDISTANT TABLE OF FUNCTION VALUES.
C
C        USAGE
C           CALL DQSF (H,Y,Z,NDIM)
C
C        DESCRIPTION OF PARAMETERS
C           H      - DOUBLE PRECISION INCREMENT OF ARGUMENT VALUES.
C           Y      - DOUBLE PRECISION INPUT VECTOR OF FUNCTION VALUES.
C           Z      - RESULTING DOUBLE PRECISION VECTOR OF INTEGRAL
C                    VALUES. Z MAY BE IDENTICAL WITH Y.
C           NDIM   - THE DIMENSION OF VECTORS Y AND Z.
C
C        REMARKS
C           NO ACTION IN CASE NDIM LESS THAN 3.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           BEGINNING WITH Z(1)=0, EVALUATION OF VECTOR Z IS DONE BY
C           MEANS OF SIMPSONS RULE TOGETHER WITH NEWTONS 3/8 RULE OR A
C           COMBINATION OF THESE TWO RULES. TRUNCATION ERROR IS OF
C           ORDER H**5 (I.E. FOURTH ORDER METHOD). ONLY IN CASE NDIM=3
C           TRUNCATION ERROR OF Z(2) IS OF ORDER H**4.
C           FOR REFERENCE, SEE
C           (1) F.B.HILDEBRAND, INTRODUCTION TO NUMERICAL ANALYSIS,
C               MCGRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP.71-76.
C           (2) R.ZURMUEHL, PRAKTISCHE MATHEMATIK FUER INGENIEURE UND
C               PHYSIKER, SPRINGER, BERLIN/GOETTINGEN/HEIDELBERG, 1963,
C               PP.214-221.
C
C     ..................................................................
C
      SUBROUTINE DQSF(H,Y,Z,NDIM)
C
C
      DIMENSION Y(1),Z(1)
      DOUBLE PRECISION Y,Z,H,HT,SUM1,SUM2,AUX,AUX1,AUX2
C
      HT=.33333333333333333D0*H
      IF(NDIM-5)7,8,1
C
C     NDIM IS GREATER THAN 5. PREPARATIONS OF INTEGRATION LOOP
    1 SUM1=Y(2)+Y(2)
      SUM1=SUM1+SUM1
      SUM1=HT*(Y(1)+SUM1+Y(3))
      AUX1=Y(4)+Y(4)
      AUX1=AUX1+AUX1
      AUX1=SUM1+HT*(Y(3)+AUX1+Y(5))
      AUX2=HT*(Y(1)+3.875D0*(Y(2)+Y(5))+2.625D0*(Y(3)+Y(4))+Y(6))
      SUM2=Y(5)+Y(5)
      SUM2=SUM2+SUM2
      SUM2=AUX2-HT*(Y(4)+SUM2+Y(6))
      Z(1)=0.D0
      AUX=Y(3)+Y(3)
      AUX=AUX+AUX
      Z(2)=SUM2-HT*(Y(2)+AUX+Y(4))
      Z(3)=SUM1
      Z(4)=SUM2
      IF(NDIM-6)5,5,2
C
C     INTEGRATION LOOP
    2 DO 4 I=7,NDIM,2
      SUM1=AUX1
      SUM2=AUX2
      AUX1=Y(I-1)+Y(I-1)
      AUX1=AUX1+AUX1
      AUX1=SUM1+HT*(Y(I-2)+AUX1+Y(I))
      Z(I-2)=SUM1
      IF(I-NDIM)3,6,6
    3 AUX2=Y(I)+Y(I)
      AUX2=AUX2+AUX2
      AUX2=SUM2+HT*(Y(I-1)+AUX2+Y(I+1))
    4 Z(I-1)=SUM2
    5 Z(NDIM-1)=AUX1
      Z(NDIM)=AUX2
      RETURN
    6 Z(NDIM-1)=SUM2
      Z(NDIM)=AUX1
      RETURN
C     END OF INTEGRATION LOOP
C
    7 IF(NDIM-3)12,11,8
C
C     NDIM IS EQUAL TO 4 OR 5
    8 SUM2=1.125D0*HT*(Y(1)+Y(2)+Y(2)+Y(2)+Y(3)+Y(3)+Y(3)+Y(4))
      SUM1=Y(2)+Y(2)
      SUM1=SUM1+SUM1
      SUM1=HT*(Y(1)+SUM1+Y(3))
      Z(1)=0.D0
      AUX1=Y(3)+Y(3)
      AUX1=AUX1+AUX1
      Z(2)=SUM2-HT*(Y(2)+AUX1+Y(4))
      IF(NDIM-5)10,9,9
    9 AUX1=Y(4)+Y(4)
      AUX1=AUX1+AUX1
      Z(5)=SUM1+HT*(Y(3)+AUX1+Y(5))
   10 Z(3)=SUM1
      Z(4)=SUM2
      RETURN
C
C     NDIM IS EQUAL TO 3
   11 SUM1=HT*(1.25D0*Y(1)+Y(2)+Y(2)-.25D0*Y(3))
      SUM2=Y(2)+Y(2)
      SUM2=SUM2+SUM2
      Z(3)=HT*(Y(1)+SUM2+Y(3))
      Z(1)=0.D0
      Z(2)=SUM1
   12 RETURN
      END