!Homework 11
!Programmer: Chad Brewbaker
!Email address: crb002@iastate.edu
!Date: March 2, 2004
!Machine used: hpc-class.iastate.edu (Intel Xenon/Myrinet cluster)
!Compiler options used:

program
implicit none
integer, parameter :: cache=96,ncache=cache*1024/8,nflush=ncache*4
real*8::flush(nflush)
integer, parameter :: n=50, ntrial=10
real*8,parameter::a=0.0,b=1.0
include "mpif.h"
real*8,dimension(ntrial)::max_time,array_time,array_perf
integer::irank,i,j,k,l,ktrial,m,p,rank,status(mpi_status_size)
integer::ier,ierror
real*8::t,xl,xr,time,sum,psum,h,x,temp,delta,exactval
real*8,allocatable,dimension(:)::procsums
call mpi_init(ierror)
call mpi_comm_size(mpi_comm_world,p,ierror)
call mpi_comm_rank(mpi_comm_world,rank,ierror)

if(rank.eq.0)then
   allocate(procsums(p))
   print *,'   '
   print *,'a=',a,' b=',b
   print *,'p=',p,' and n=',n
   psum =0.0d0
   do i=1,10000
      delta =abs(a-b)/dble(10000)
      x= a+(delta*dble(i))
      psum=psum+(1.0d0/(1.0d0+(x*x)))
   enddo
   exactval=psum/dble(10000)
   print *,'The exact value of the integral=',exactval
endif

do ktrial=1,ntrial

   flush=flush+.0000123
   !start timer
   call MPI_Barrier(mpi_comm_world,ierror)
   t=mpi_wtime()


   !temp is where I start my part of the integral
   temp=((abs(a-b)/dble(p))*dble(rank))+a
   delta=(abs(a-b)/dble(p))/dble(n)
   psum=0.0d0
   do i=0,n-1
      x=(dble(i)*delta)+temp
      sum=1.0d0/(1.0d0+(x*x))
      x=(dble(i+1)*delta)+temp
      sum=sum+(1.0d0/(1.0d0+(x*x)))
      sum=sum/2
      psum=psum+sum
   enddo
   psum=psum/dble(n)
   
   !send my result to root
   call MPI_Gather(psum,1,MPI_REAL8,procsums,1,MPI_REAL8,0,MPI_COMM_WORLD,ierror)
   
   
   if(rank.eq.0)then
      psum=0.0d0
      do i=1,p
         psum=psum+procsums(i)
      enddo
      psum=psum/dble(p)
   endif
   
   array_time(ktrial)=MPI_Wtime()-t
   call mpi_barrier(mpi_comm_world,ierror)
   flush(ktrial)=flush(ktrial)+56
   !Check for errors when k=1
   if(ktrial == 1)then
      if(rank==0)then
         print *,'Computed value=',psum
         print *,'Error (result-exact)=', psum-exactval
      endif
   endif
enddo




call mpi_reduce(array_time,max_time,ntrial,mpi_real8,mpi_max,0,mpi_comm_world,ierror)

if(rank==0)then
   print *,'  '
   print *,'Time in seconds '
   do ktrial=1,ntrial
      print*,max_time(ktrial)
   enddo
   print *,'   '
   print *,'   '
endif

!!Quit the program
call MPI_Finalize(ierror)
end program

!!$vincent% mpirun -np 8 hw11.exe
!!$    
!!$ a=   0.00000000000000D+000  b=   1.00000000000000     
!!$ p=           8  and n=          50
!!$ The exact value of the integral=  0.785373162980782     
!!$ Computed value=  0.785397902980782     
!!$ Error (result-exact)=   2.47399999999676D-005
!!$   
!!$ Time in seconds 
!!$   9.98973846435547D-005
!!$   4.79221343994141D-005
!!$   4.69684600830078D-005
!!$   4.50611114501953D-005
!!$   4.50611114501953D-005
!!$   4.31537628173828D-005
!!$   4.79221343994141D-005
!!$   4.79221343994141D-005
!!$   4.91142272949219D-005
!!$   4.48226928710938D-005
!!$