subroutine psimp ( a, b, n, fn, gn, vint )

c*********************************************************************72
c
c  this subroutine uses the product type simpson rule
c  compounded n times to approximate the integral from a to b
c  of the function fn(x) * gn(x).  fn and gn are function
c  subprograms which must be supplied by the user.  the
c  result is stored in vint.
c
      double precision a, ag, am(3,3), b, f(3), fn, g(3),
     &  gn, h, vint, x(2), dble

      data am(1,1), am(3,3) /2 * 4.d0 /, am(1,2), am(2,1),
     &  am(2,3), am(3,2) / 4 * 2.d0/, am(1,3), am(3,1)
     &  /2 * -1.d0/, am(2,2) / 16.0d0 /

      h = ( b - a ) / dble ( float ( n ) )
      x(1) = a + h / 2.d0
      x(2) = a + h
      vint = 0.d0
      f(3) = fn ( a )
      g(3) = gn ( a )

      do i = 1, n

        f(1) = f(3)
        g(1) = g(3)
        do j = 1, 2
          f(j+1) = fn ( x(j) )
          g(j+1) = gn ( x(j) )
          x(j) = x(j) + h
        end do

        do j = 1, 3
          ag = 0.d0
          do k = 1, 3
            ag = ag + am(j,k) * g(k)
          end do
          vint = vint + f(j) * ag
        end do

      end do

      vint = h * vint / 30.d0

      return
      end