!Homework 9
!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=500, 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::rooty(1:n),y(1:n),x(1:n),z(1:n),t
real*8,allocatable,dimension(:,:)::temp,rootA,A
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))
allocate(A(m,n))

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
!**Begin master section for root processor
   call random_number(rootA);
   call random_number(x)
   call random_number(rooty)
   call DGEMV('n',n,n,1.0d0,rootA,n,x,1,1.0d0,z,1)
   z=z+rooty
 !print *,z
endif

call random_number(flush)

do ktrial=1,ntrial
   flush=flush+.0000123
   !start timer
   call MPI_Barrier(mpi_comm_world,ierror)
   t=mpi_wtime()
   !broadcast x
   call MPI_Bcast(x,n,MPI_REAL8,0,MPI_COMM_WORLD,ierror)
   
   !Send off all chunks of Ay
   if(rank.eq.0)then
      do i=2,p
         l=(i-1)*m
         !send off a block of m rows of Y
         do k=1,n
            do j=1,m
               temp(j,k)=rootA(j+l,k)
            enddo
         enddo
         do j=1,m
            temp(j,n+1)=rooty(j+l)
         enddo
            call MPI_Send(temp,m*(n+1),MPI_REAL8,i-1,1,MPI_COMM_WORLD,ierror)
      enddo
      !do same for root
      do k=1,n
         do j=1,m
            temp(j,k)=rootA(j,k)
         enddo
      enddo
      do j=1,m
         temp(j,n+1)=rooty(j)
      enddo
   else
      call MPI_Recv(temp,m*(n+1),MPI_REAL8,0, 1, MPI_COMM_WORLD,status,ierror)
   endif
   do i=1,m
      yb(i)=temp(i,n+1)
   enddo
   
   !calculate my portion of y=y+Ax
   
   !call DGEMV(trans,rows,cols,alpha,A,LDA,X,INCX,beta,Y,INCY)
   call DGEMV('n',m,n,1.0d0,temp,m,x,1,1.0d0,yb,1)

   !send result to root
   if(rank.ne.0)then
      call MPI_Send(yb,m,MPI_REAL8,0,1,MPI_COMM_WORLD,ierror)
   else
      do j=1,m
         y(j)=yb(j)
      enddo
      do i=2,p
         call MPI_Recv(yb,m,MPI_REAL8,MPI_ANY_SOURCE, 1, MPI_COMM_WORLD,status,ierror)
         do j=1,m
            y(j+(status(MPI_SOURCE)*m))=yb(j) 
         enddo
      enddo
   endif

   !call timer
   array_time(ktrial)=MPI_Wtime()-t
   call mpi_barrier(mpi_comm_world,ierror)
   flush(ktrial)=flush(ktrial)+y(ktrial)
   !Check for errors when k=1
   if(ktrial == 1)then
      if(rank==0)then
         print *,'Max absolute error =',maxval(abs(y-z))
      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

!!$    
!!$ This program assumes p divides n:
!!$ p=           4  and n=         500
!!$ Max absolute error =   5.68434188608080D-014
!!$   
!!$ Time in seconds 
!!$   2.23548412322998D-002
!!$   9.24777984619141D-003
!!$   9.32192802429199D-003
!!$   9.28401947021484D-003
!!$   9.24110412597656D-003
!!$   9.29188728332520D-003
!!$   9.25111770629883D-003
!!$   9.29212570190430D-003
!!$   9.26995277404785D-003
!!$   9.33289527893066D-003
!!$