program cell_mpi

  use mpi

  implicit none

  integer :: Ndim,Niter
  integer :: nlocal,offset
  integer :: x,n,maxval=0,local_maxval=0
  real*8 :: sum=0,local_sum=0
  integer :: numprocs,myrank,ierr
  integer :: rank_left,rank_right
  integer :: tag1=1,tag2=2
  character(len=80) :: paramin
  integer,pointer :: tmp(:),buffer(:),nextbuffer(:)

  integer :: status(MPI_STATUS_SIZE)
  integer :: request1,request2

  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,numprocs,ierr)
  call mpi_comm_rank(mpi_comm_world,myrank,ierr)

  ! Exit if # of arguments < 2
  if (command_argument_count().lt.2) then
     if (myrank.eq.0) then
        write(*,*) "Usage: cell_mpi <Ndim> <Niter>"
     endif
     call mpi_finalize(ierr)
     stop
  endif

  call getarg(1,paramin)
  read(paramin,*) Ndim
  call getarg(2,paramin)
  read(paramin,*) Niter

  if (myrank.eq.0) then
     write(*,*)
     write(*,'(A,T24,I8)') "Size of the array: ",Ndim
     write(*,'(A,T24,I8)') "Number of iterations: ",Niter
  endif

  ! Calculate the size of the local array.
  nlocal=ndim/numprocs
  offset=nlocal*myrank

  allocate(buffer(nlocal+2))
  allocate(nextbuffer(nlocal+2))
  allocate(tmp(nlocal+2))

  ! Initialize the array.
  do x=1,nlocal+2
     buffer(x)=mod(x+offset-2,10)
  enddo
  if (myrank.eq.0) then
     buffer(1)=mod(ndim-1,10)
  elseif (myrank.eq.numprocs-1) then
     buffer(nlocal+2)=0
  endif

  rank_left=myrank-1
  rank_right=myrank+1

  ! Main loop
  do n=1,Niter
     do x=2,nlocal+1
        nextbuffer(x)=mod(buffer(x-1)+buffer(x+1),10)
     enddo
     ! Swap the pointers.
     tmp => buffer
     buffer => nextbuffer
     nextbuffer => tmp
     ! Exchange boundary elements with neighbors.
     call mpi_isend(buffer(2),1,MPI_INTEGER,rank_left,tag1,mpi_comm_world,request1,ierr)
     call mpi_isend(buffer(nlocal+1),1,MPI_INTEGER,rank_right,tag2,mpi_comm_world,request2,ierr)
     call mpi_recv(buffer(1),1,MPI_INTEGER,rank_left,tag1,mpi_comm_world,status,ierr)
     call mpi_recv(buffer(nlocal+2),1,MPI_INTEGER,rank_right,tag2,mpi_comm_world,status,ierr)
     call mpi_wait(request1,status,ierr)
     call mpi_wait(request2,status,ierr)
  enddo

  ! Find the local maximum and sum.
  do x=2,nlocal+1
     local_sum=local_sum+buffer(x)
     if (buffer(x).gt.local_maxval) local_maxval=buffer(x)
  enddo

  ! Reduce the result to the root process
  call mpi_reduce(local_sum,sum,1,MPI_INTEGER,MPI_SUM,0,mpi_comm_world,ierr)
  call mpi_reduce(local_maxval,maxval,1,MPI_INTEGER,MPI_MAX,0,mpi_comm_world,ierr)
  
  if (myrank.eq.0) then
     write(*,*)
     write(*,'(A,T24,I4)') "Maximum: ",maxval
     write(*,'(A,T24,F10.6)') "Average: ",sum/float(Ndim)
     write(*,*)
  endif

  tmp => null()
  buffer => null()
  nextbuffer => null()

  call mpi_finalize(ierr)

  stop
end program cell_mpi