      MODULE m_qsf

!      interface qsf
!        module procedure qsfReal, qsfComplex
!      end interface

      CONTAINS
      
      SUBROUTINE qsf(h,y,z,ndim,isave)
c     ..................................................................
c     subroutine qsf
c     purpose
c        to integrate an equidistant table of function values.
c     description of parameters
c         h     - the increment of arguement values.
c         y     - the input vector of function values.  y must be
c                 dimensioned ndim.
c         z     - if isave = 0,z is the final value of the integration.
c                 if isave not 0,z is the resulting vector of integrals
c                 for each y value. in this case it must be dimensioned
c                 ndim.
c         ndim  - the number of points in the table.
c         isave - determines whether a single value or multiple results
c                 will be returned
c                  isave =  0  - only the final value of the integra-
c                                tion will be returned.
c                  isave ne 0  - the results of integration at each
c                                point in the table will be returned.
c     remarks
c        no action in case ndim less than 3.
c     subroutines and function subprograms required
c        none
c     method
c        integration is done by means of simpsons rule together with
c        newtons 3/8 rule or a combination of these two rules. trunca-
c        tion error is of order h**5 if ndim greater than 3. if ndim = 3
c        truncation error is of order h**4.
c     reference
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 mathematic fuer ingenieure und
c            physiker, springer, berlin/goettingen/heidelberg, 1963,
c            pp. 214-221.
c     ..................................................................
C     .. Scalar Arguments ..

      IMPLICIT NONE
      REAL h
      INTEGER isave,ndim
C     ..
C     .. Array Arguments ..
      REAL y(ndim),z(*)
C     ..
C     .. Local Scalars ..
      REAL aux,aux1,aux2,ht,sum1,sum2,val
      INTEGER i
      LOGICAL save
C     ..
      save = isave .NE. 0
      ht = h / 3.0
      IF ( ndim == 5 ) THEN
        GOTO 140
      ELSEIF ( ndim < 5 ) THEN
        GOTO 130
      ENDIF
c     ndim is greater than 5. preparations of integration loop
      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.875e0* (y(2)+y(5))+2.625e0* (y(3)+y(4))+y(6))
      IF (.NOT.save) GO TO 30
   20 z(1) = 0.e0
      sum2 = y(5) + y(5)
      sum2 = sum2 + sum2
      sum2 = aux2 - ht* (y(4)+sum2+y(6))
      aux = y(3) + y(3)
      aux = aux + aux
      z(2) = sum2 - ht* (y(2)+aux+y(4))
      z(3) = sum1
      z(4) = sum2
   30 CONTINUE
      IF ( ndim <= 6 ) THEN
         GOTO 70
      ENDIF
c     integration loop
   40 DO 60 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))
         IF (save) z(i-2) = sum1
         IF ( i >= ndim ) THEN
           GOTO 100
         ENDIF
         aux2 = y(i) + y(i)
         aux2 = aux2 + aux2
         aux2 = sum2 + ht* (y(i-1)+aux2+y(i+1))
         IF (save) z(i-1) = sum2
   60 CONTINUE
   70 IF (.NOT.save) GO TO 90
   80 z(ndim-1) = aux1
      z(ndim) = aux2
      RETURN
   90 z(1) = aux2
      RETURN
  100 IF (.NOT.save) GO TO 120
  110 z(ndim-1) = sum2
      z(ndim) = aux1
      RETURN
  120 z(1) = aux1
      RETURN
c     end of integration loop
  130 CONTINUE
      IF ( ndim == 3 ) THEN
        GOTO 210
      ELSEIF ( ndim < 3 ) THEN
        GOTO 230
      ENDIF
c     ndim is equal to 4 or 5
  140 sum2 = 1.125e0*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))
      IF (.NOT.save) GO TO 160
  150 z(1) = 0.e0
      aux1 = y(3) + y(3)
      aux1 = aux1 + aux1
      z(2) = sum2 - ht* (y(2)+aux1+y(4))
  160 CONTINUE
      IF ( ndim < 5 ) THEN
        GOTO 190
      ENDIF
  170 aux1 = y(4) + y(4)
      aux1 = aux1 + aux1
      val = sum1 + ht* (y(3)+aux1+y(5))
      IF (save) GO TO 180
      z(1) = val
      RETURN
  180 z(5) = val
      GO TO 200
  190 IF (save) GO TO 200
      z(1) = sum2
      RETURN
  200 z(3) = sum1
      z(4) = sum2
      RETURN
c     ndim is equal to 3
  210 sum2 = y(2) + y(2)
      sum2 = sum2 + sum2
      val = ht* (y(1)+sum2+y(3))
      IF (save) GO TO 220
      z(1) = val
      RETURN
  220 z(1) = 0.e0
      z(2) = ht* (1.25e0*y(1)+y(2)+y(2)-.25e0*y(3))
      z(3) = val
  230 RETURN
      END SUBROUTINE qsf

      subroutine qsfComplex( h, y, z, ndim, isave )

        implicit none

        REAL h
        INTEGER isave,ndim
        COMPLEX y(ndim),z(*)

        real zr(ndim), zi(ndim)

        call qsf( h,  real(y), zr, ndim, isave )
        call qsf( h, aimag(y), zi, ndim, isave )
        z(1:ndim) = cmplx( zr(1:ndim), zi(1:ndim) )

      end subroutine qsfComplex

      END MODULE m_qsf
