! workshare_ex.f90
include "mpiSim.f90"
program worksharing
include "useSim.f90"

!! calculate PI using collective communication with ..
!! .. more sophiticated work sharing

  implicit none

! include 'mpif.h'

  integer   ::         myrank, nprocs, n, i, islave, master=0, count
  integer   ::         chunk, rest, is, ie, ilen
  integer   ::         status(MPI_STATUS_SIZE)
  integer   ::         ierr, resultlen, tag
  character (LEN=8) :: hostname
  real*8    ::         t0, t1, t2
  integer, dimension(:), allocatable :: results
  real*8, parameter :: PI25DP = 3.141592653589793238462643d0
  ! integer chunk, ierr, ilen, rest, is, ie
  real*8 :: mypi
  double precision f, a, h, sum, x, pi

  f(a) = 4.d0 / (1.d0 + a*a) ! function to integrate

  ! Initialize MPI
  Call MPI_Init(Ierr)

  ! Get my rank id and number of MPI processes
  call MPI_Comm_Rank(MPI_COMM_WORLD, myrank, ierr)
  call MPI_Comm_Size(MPI_COMM_WORLD, nprocs, ierr)

  ! Get hostname
  call MPI_Get_processor_name(hostname, resultlen, ierr)
  write (*,*) "Rank ", myrank, " out of ", nprocs,  &
    " processes running on ", hostname(1:index(hostname,".")-1)


  ! Synchonization point (not really necessary here)
  t0 = MPI_Wtime()
  call MPI_Barrier (MPI_COMM_WORLD, ierr)
  t1 = MPI_Wtime()

  !-----------------------------------------------------------
  ! Do work here ...

  do

    ! The master reads an integer number form standard input and sends it to all slaves
    if ( myrank .eq. master ) then ! Master part
      write(*,*) 'Number of intervals ? (0 quits): '
      read (*,*) n
    end if
    call MPI_Bcast ( n, 1, MPI_INTEGER, master, MPI_COMM_WORLD, ierr )
    if ( n .le. 0 .or. n > 10000000 ) exit

     h = 1.0d0/dble(n)   ! even step size h as a function of partions
     sum  = 0.0d0        ! zeroise sum

     ! improved approach of work sharing
     rest = mod(n,nprocs)
     chunk = ( n - rest ) / nprocs
     if ( myrank .lt. rest ) then
       is = 1 + myrank * ( chunk + 1 )
       ie = is + chunk
     else
       is = 1 + myrank * chunk + rest
       ie = is + chunk - 1
     end if
     ilen = ie-is+1

     write(*,*) 'nprocs = ',nprocs,' myrank = ',myrank,' is = ',is,' ie = ',ie, &
          ' ilen = ',ilen

     do i = is, ie
        x = h * (dble(i) - 0.5d0)
        sum = sum + f(x)
     end do
     mypi = h * sum

     call MPI_Reduce(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0, &
          MPI_COMM_WORLD,ierr) ! collect all the sums

     if ( myrank .eq. 0 ) then
        write(*,'(a,F18.16,a,F18.16)') 'pi is approximately: ',pi, &
             ' Error is: ', abs(pi - PI25DP)
     end if

  end do

  !-----------------------------------------------------------

  ! Synchronization point (not really necessary here)
  call MPI_Barrier (MPI_COMM_WORLD, ierr)
  t2 = MPI_Wtime()
  write (*,*) "process ", myrank, ": ", t2-t1, " seconds"

  ! Shut Down MPI
  call MPI_Finalize(ierr )

end program worksharing