!Homework 10
!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=1024, ntrial=10
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::A(n,n),rooty(1:n),y(1:n),x(1:n),u(1:n),t
real*8,allocatable,dimension(:,:)::temp
real*8,allocatable,dimension(:)::yb
call mpi_init(ierror)
call mpi_comm_size(mpi_comm_world,p,ierror)
call mpi_comm_rank(mpi_comm_world,rank,ierror)

m=n/p
allocate(temp(1:m,1:n+1))
allocate(yb(1:m))

if(rank.eq.0)then
   allocate(rootA(n,n))
endif

!Check to see if p divides n:
if(rank.eq.0) then
   print *,'   '
   print *,'This program assumes p divides n:'
   print *,'p=',p,' and n=',n,'m=n/p',m
endif
call random_number(A)
A=A+float(rank)
call random_number(x)
x=x+float(rank)

!initalize y on processor 0 and broadcast
if(rank.eq.0)then
   call random_number(y)
endif
call MPI_Bcast(y,n,MPI_REAL8,0,MPI_COMM_WORLD,ierror)
!store original number for checking
u=y
flush=0.0d0
do k=1,ntrial
   !preserve the inital value of y
   y=u
   flush=flush+0.1
   call MPI_Barrier(MPI_COMM_WORLD,ierror)
   t=MPI_Wtime()
   call DGEMV('n',n,n,1.0d0,rootA,n,x,1,1.0d0,y,1)
   
   array_time(ktrial)=MPI_Wtime()-t
   call MPI_Barrier(MPI_COMM_WORLD,ierror
   
   !prevent compiler from splitting flush
   flush(ktrial)=flush(ktrial)+y(ktrial)
enddo

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


!To check answers gather the A arrays to B, and gather the x arrays to v on root
call MPI_Gather(x,n,MPI_REAL8,v,n,MPI_REAL8,0,MPI_COMM_WORLD,ierror)


call MPI_Gather(A,n*n,MPI_REAL8,B,n*n,MPI_REAL8,0,MPI_COMM_WORLD,ierror)


if(rank==0)then
   print*,'error=',maxval(abs(y-matmul(B(,)))
   print*,'   '
   print*,'   '
   print*,'Time in seconds          MFLOPS'
   do ktrial=1,ntrial
      print*,max_time(ktrial),' ',array_perf(ktrial)
   enddo
   print*,'   '
   print*,'   '
endif

!To mess with the compiler optimizations
if(rank==0)then
   print*,'flush(2)+array_time(5)+y(7)=',flush(2)+array_time(5)+y(7)
endif

call MPI_Finalize(ierror)
end program