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

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

      data am(1,1), am(2,3) / 2 * 1.718245836551854d0 /,
     &  am(1,2), am(2,2) / 2 * 1.d0 /, am(1,3), am(2,1)
     &  / 2 * -.2182458365518542d0 /

      h = ( b - a ) / dble ( float ( n ) )
      x(1) = a + .1127016653792583d0 * h
      x(2) = a + .8872983346207417d0 * h
      y(1) = a + h / 2.d0
      y(2) = a + h
      vint = 0.d0
      g(3) = gn ( a )
      do i = 1, n
        ag = fn ( y(1) )
        g(1) = g(3)
        do j = 1, 2
          f(j) = fn ( x(j) )
          g(j+1) = gn ( y(j) )
          x(j) = x(j) + h
          y(j) = y(j) + h
        end do
        vint = vint + ag * 4.d0 * g(2)
        do j = 1, 2
          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 / 9.d0

      return
      end