c-----This program calculates the value of pi, using numerical integration
c-----and clocks the solution time.
c-----The user selects the number of points of integration and the limits a & b
c-----PI program evaluates the integral of 4/(1+x*x) from a to b
c-----by dividing the area under the curve [(b - a) / #points] slices.

c-----PI program loops over slices in interval, summing
c-----   area of each slice

c   Variables:
c
c       points          Number of points used in the integration. Input by the user
c       a               Lower limit of integration
c       b               Upper limit of integration
c       x               Independent variable of integration
c       f               Function to integrate (external)
c       slicesize       Size of integration slice
c       partialint      Partial integral calculated on node
c       integral        Calculated integral

c
c   Calls:
c
c       userinpt        to get number of points of integration from user
c       pi              calculate the value of pi


	PROGRAM PI

        external  userinpt
        integer*4 userinpt

        external  f
        double precision f

        integer*4 i,j, points
        double precision a, b
	double precision slicesize, x, integral
        double precision mflops

        real tt1(2),tt2(2)

c
c       Receive points of integration from the user
c
 120    j = userinpt(a, b, points)
        if (j.le.0) goto 200

        toh1 = etime(tt1)
        toh2 = etime(tt2)
        toh  = toh2 - toh1
c
c------Compute width of cut
c
        integral = 0.d0

c       calculate slice size
        slicesize = (b - a) / points

c
c------Loop over interval, suming areas
c
        t1 = etime(tt1)
        do 100 i=1,points
         x = a + (i - 0.5)*slicesize
         integral = integral+ slicesize*f(x)
 100    continue
c
c------Finish overall timing and write results
c
        t2 = etime(tt2)
c------A multiply and/or an add counts as one FLOP. Divide counts as 3 FLOPs.
        mflops = points*7/((t2-t1-toh)*10**6)
        write(6,*) 'Time in main = ',t2-t1-toh,' MFLOPS = ',mflops
        write(6,*) 'Sum = ',integral,
     &   ', Error = ',integral-3.14159265358979323846d0
        
       goto 120

 200   stop
       end

