c****************************************************************************
c   pi_hpf.hpf - compute pi by integrating f(x) = 4/(1 + x**2) from 0 to 1
c****************************************************************************
      program pi_hpf

c       This is the exact value of pi: 
        double precision  PI25DT

c       We declare all the variables here:
c       (their meaning is explained above)
        double precision     pi, h, gl_sum, x
        integer              n,numprocs,i

c       These real variables are used to store the elapsed time:
        real t(2),t_start,t_end

c       These are HPF directives that declare the processor arrays:
!hpf$   processors, dimension(number_of_processors()) :: all_procs
!hpf$   processors :: one_proc

c       Some of the variables are only needed on a single node
c       others are needed on all nodes (replicated):
!adp$   replicated onto all_procs ::  n
!adp$   single :: t_start,t_end,t,PI25DT,numprocs,pi

c       We must tell the compiler what kind of function ETIME is
c       in this case a serial fortran 77 routine:
        interface
          extrinsic(f77_serial) function etime(t)
            real :: etime,t(2)            
          end function etime
        end interface  

        interface
          extrinsic(f77_serial) subroutine printout()
          end subroutine printout
        end interface
        
        PI25DT = 3.141592653589793238462643D0

c       The number of processors is:        
        numprocs=number_of_processors()
        write(*,*) 'MPI initialized with ',numprocs,' processors' 
        write(*,*) 'on a processor grid of shape ',processors_shape()
                  
c       We let the user enter the number of intervals as many times
c       as he wants in an infinite (empty) DO loop: 
        do
       
           write(*,*) 'Enter the number of intervals (0 to quit) '
           call printout()
           read(*,*) n
             
c          We start timing the program on node 0 using ETIME
           t_start=etime(t)

c          Check for quit signal--the user entered 0 for n
           if(n.le.0) exit
           
c          The step size of the intervals is simply 1/n
           h = 1.0d0/n

c          Each of the nodes does its own sum of the rectangles
c          that belong to it (every numprocs rectangle): 
           gl_sum  = 0.0d0
c          We do this by using the independent
c          and reduction directive:           
!hpf$      independent, new(i,x),reduction(gl_sum)           
           do i = 1, n 
             x = h * (dble(i) - 0.5d0)
             gl_sum = gl_sum + 4.d0/(1.d0 + x*x)
           end do                       
           
           pi= h * gl_sum
           
           write(*,*) 'Approximation for pi is ', pi
           write(*,*) 'The error is thus', abs(pi - PI25DT)
c          We time the end of the run and print the elapsed time:
           t_end=ETIME(t)
           write(*,*) 'Elapsed time on node 0 is', t_end-t_start 
           call printout()         

c       We end the infinite DO loop:
        end do

      end program pi_hpf