C    -----------------------------------------------------------------------
C    This program implements the concurrent wave equation described 
C    in Chapter 5 of Fox et al., 1988, Solving Problems on Concurrent
C    Processors, vol 1.  
C
C    A vibrating string is decomposed into points.  Each processor is 
C    responsible for updating the amplitude of a number of points over
C    time.
C	
C    At each iteration, each processor exchanges boundary points with
C    nearest neighbors.  This version uses low level sends and receives
C    to exchange boundary points.
C
C  AUTHOR: Roslyn Leibensperger (C program for MPL)
C  REVISED: 06/07/93 R. Leibensperger     Ported to Fortran
C  CONVERTED to MPI: 11/12/94 by Xianneng Shen
C  ------------------------------------------------------------------------ 
C  ------------------------------------------------------------------------ 
C  Explanation of constants and variables used in common blocks and
C  include files
C    MASTER           = task ID of master
C    INTSIZE          = size of integer in bytes
C    REAL8SIZE        = size of real*8 in bytes
C    E_OUT1, E_OUT2   = message types
C    taskid           = task ID
C    nproc            = number of tasks
C    tpoints          = total points along wave
C    nsteps           = number of time steps
C    npoints          = number of points handled by this task
C    first            = index of first point handled by this task
C    values(0:1001)   = values at time t
C    oldval(0:1001)   = values at time (t-dt)
C    newval(0:1001)   = values at time (t+dt)
C  ------------------------------------------------------------------------ 

      program wave
      implicit none
      include 'mpif.h'
      include 'parameters.h'
      integer ierr
      integer taskid, nproc
      common/config/taskid, nproc

      integer left, right, nbuf(4)

C     Learn number of tasks and taskid
      call mpi_init(ierr)
      call mpi_comm_rank(MPI_COMM_WORLD, taskid, ierr)
      call mpi_comm_size(MPI_COMM_WORLD, nproc, ierr)
      write (*,5) taskid
 5    format (I5, ': Wave Program running')

C     Learn value of allgrp and dontcare

C     Determine left and right neighbors
      if (taskid .eq. nproc-1) then
         right = 0
      else 
         right = taskid + 1
      end if

      if (taskid .eq. 0) then
         left = nproc - 1
      else 
         left = taskid - 1
      end if

C     Get program parameters and initialize wave values
      if (taskid .eq. MASTER) then
         call init_master
      else
         call init_workers
      end if
      call init_line

C     Update values along the wave for nstep time steps
      call update(left, right)

C     Master collects results from workers and prints
      if (taskid .eq. MASTER) then
         call out_master
      else 
         call out_workers
      end if
      call mpi_finalize(ierr)
      end

C  ------------------------------------------------------------------------
C     Master obtains input values from user
C  ------------------------------------------------------------------------ 

      subroutine init_master

      implicit none
      include 'parameters.h'
      include 'mpif.h'
      integer ierr
      integer taskid, nproc
      common/config/taskid, nproc

      integer tpoints, nsteps
      common/inputs/tpoints, nsteps

      integer MAXPOINTS, MAXSTEPS
      parameter (MAXPOINTS = 1000)
      parameter (MAXSTEPS = 10000)
      integer buffer(2), msglen

      tpoints = 0
      nsteps = 0
 
      do while ((tpoints .lt. nproc) .or. (tpoints .gt. MAXPOINTS))
         write (*,*)'Enter number of points along vibrating string'
         read (*,*) tpoints
         if ((tpoints .lt. nproc) .or. (tpoints .gt. MAXPOINTS))
     &     write (*,*) 'enter value between ',nproc,' and ',MAXPOINTS
      end do

      do while ((nsteps .lt. 1) .or. (nsteps .gt. MAXSTEPS))
         write (*,*) 'Enter number of time steps'
         read (*,*) nsteps
         if ((nsteps .lt. 1) .or. (nsteps .gt. MAXSTEPS))
     &        write (*,*) 'enter value between 1 and ', MAXSTEPS
      end do

      write (*,10) taskid, tpoints, nsteps
 10   format(I5, ': points = ', I5, ' steps = ', I5)

C     Broadcast total points, time steps
      buffer(1) = tpoints
      buffer(2) = nsteps
      call mpi_bcast(buffer, 2, MPI_INTEGER, 0, MPI_COMM_WORLD,ierr)
      end

C  -------------------------------------------------------------------------
C     Workers receive input values from master
C  -------------------------------------------------------------------------
 
      subroutine init_workers

      implicit none
      include 'parameters.h'
      include 'mpif.h'
      integer ierr
      integer taskid, nproc
      common/config/taskid, nproc

      integer tpoints, nsteps
      common/inputs/tpoints, nsteps

      integer buffer(2), msglen

C     Receive time advance parameter, total points, time steps
      call mpi_bcast(buffer, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
      tpoints = buffer(1)
      nsteps =  buffer(2)

      end

C  ------------------------------------------------------------------------
C     Initialize points on line
C  -----------------------------------------------------------------------

      subroutine init_line

      implicit none

      integer taskid, nproc
      common/config/taskid, nproc

      integer tpoints, nsteps
      common/inputs/tpoints, nsteps

      integer npoints, first
      common/decomp/npoints, first

      real*8 values(0:1001), oldval(0:1001), newval(0:1001)
      common/data/values, oldval, newval

      real*8 PI
      parameter (PI = 3.14159265)
      integer nmin, nleft, npts, i, j, k
      real*8 x, fac

C     Calculate initial values based on sine curve
      nmin = tpoints/nproc
      nleft = mod(tpoints, nproc)
      fac = 2.0 * PI

      k = 0
      do i = 0, nproc-1
         if (i .lt. nleft) then
            npts = nmin + 1
         else
            npts = nmin
         endif
         if (taskid .eq. i) then
            first = k + 1
            npoints = npts
	    write (*,15) taskid, first, npts
 15         format (I5, ': first = ', I5, ' npoints = ', I5)
            do j = 1, npts
               x = float(k)/float(tpoints - 1)
               values(j) = sin (fac * x)
               k = k + 1
            end do
         else 
            k = k + npts
         end if
      end do
      do i = 1, npoints
         oldval(i) = values(i)
      end do
      end

C  -------------------------------------------------------------------------
C     Calculate new values using wave equation
C  -------------------------------------------------------------------------
 
      subroutine do_math(i)
      implicit none
      integer i

      integer tpoints, nsteps
      common/inputs/tpoints, nsteps

      real*8 values(0:1001), oldval(0:1001), newval(0:1001)
      common/data/values, oldval, newval

      real*8 dtime, c, dx, tau, sqtau

      dtime = 0.3
      c = 1.0
      dx = 1.0
      tau = (c * dtime / dx)
      sqtau = tau * tau
      newval(i) = (2.0 * values(i)) - oldval(i) 
     &     + (sqtau * (values(i-1) - (2.0 * values(i)) + values(i+1)))
      end

C  -------------------------------------------------------------------------
C     Update all values along line a specified number of times 
C  -------------------------------------------------------------------------
 
      subroutine update(left, right)
      implicit none
      include 'mpif.h'
      integer left, right
      integer ierr, status(MPI_STATUS_SIZE), request
      include 'parameters.h'

      integer npoints, first
      common/decomp/npoints, first

      integer tpoints, nsteps
      common/inputs/tpoints, nsteps

      real*8 values(0:1001), oldval(0:1001), newval(0:1001)
      common/data/values, oldval, newval

      integer E_RtoL, E_LtoR
      parameter (E_RtoL = 10)
      parameter (E_LtoR = 20)
      integer i, j, id_rtol, id_ltor, nbytes, msglen

C     Update values for each point along string
      do i = 1, nsteps

C        Exchange data with "left-hand" neighbor
         if (first .ne. 1) then
            call mpi_isend(values(1), 1, MPI_DOUBLE_PRECISION, left,
     .       E_RtoL, MPI_COMM_WORLD, request, ierr)
            call mpi_wait(request, status, ierr)
            call mpi_recv(values(0), 1, MPI_DOUBLE_PRECISION, left,
     .       E_LtoR, MPI_COMM_WORLD, status, ierr)
         end if
C        Exchange data with "right-hand" neighbor
	 if (first+npoints-1 .ne. tpoints) then
            call mpi_isend(values(npoints), 1, MPI_DOUBLE_PRECISION,
     .       right, E_LtoR, MPI_COMM_WORLD, request, ierr)
            call mpi_wait(request, status, ierr)
            call mpi_recv(values(npoints+1), 1, MPI_DOUBLE_PRECISION,
     .       right, E_RtoL, MPI_COMM_WORLD, status, ierr)
         end if

C        Update points along line
         do j = 1, npoints
C           Global endpoints
            if ((first+j-1 .eq. 1).or.(first+j-1 .eq. tpoints))then
               newval(j) = 0.0
	    else
               call do_math(j)
            end if
         end do

         do j = 1, npoints
            oldval(j) = values(j)
	    values(j) = newval(j)
         end do

      end do

      end

C  ------------------------------------------------------------------------
C     Receive results from workers and print
C  ------------------------------------------------------------------------ 

      subroutine out_master

      implicit none
      include 'parameters.h'
      include 'mpif.h'
      integer ierr, status(MPI_STATUS_SIZE), request
      integer taskid, nproc
      common/config/taskid, nproc

      integer tpoints, nsteps
      common/inputs/tpoints, nsteps

      integer npoints, first
      common/decomp/npoints, first

      real*8 values(0:1001), oldval(0:1001), newval(0:1001)
      common/data/values, oldval, newval

      integer i,  start, npts, buffer(2), tpts
      real*8 results(1000) 
 
C     Store worker's results in results array
      do i = 1, nproc - 1
C        Receive number of points and first point
         call mpi_recv(buffer, 2, MPI_INTEGER, MPI_ANY_SOURCE,
     .       E_OUT1, MPI_COMM_WORLD, status, ierr)
         start = buffer(1)
         npts = buffer(2)

C        Receive results
         call mpi_recv(results(start), npts, MPI_DOUBLE_PRECISION,
     .       MPI_ANY_SOURCE, E_OUT2,  MPI_COMM_WORLD, status, ierr)
      end do

C     Store master's results in results array
      do i = first, first+npoints-1
         results(i) = values(i)
      end do

      if (tpoints .lt. 10) then
         tpts = tpoints
      else 
         tpts = 10
      end if
      write (*,200) tpts, (results(i), i = 1, tpts)
 200  format('first ', I5, ' points (for validation):'/
     &       10(f4.2, '  '))
      end

C  -------------------------------------------------------------------------
C     Send the updated values to the master
C  -------------------------------------------------------------------------
 
      subroutine out_workers

      implicit none
      include 'mpif.h'
      integer ierr, status(MPI_STATUS_SIZE), request
      include 'parameters.h'

      integer npoints, first
      common/decomp/npoints, first

      real*8 values(0:1001), oldval(0:1001), newval(0:1001)
      common/data/values, oldval, newval

      integer buffer(2)


C     Send first point and number of points handled to master
      buffer(1) = first
      buffer(2) = npoints
      call mpi_isend(buffer, 2, MPI_INTEGER, 0, E_OUT1,
     .        MPI_COMM_WORLD, request, ierr)
      call mpi_wait(request, status, ierr) 

C     Send results to master
      call mpi_isend(values(1), npoints, MPI_DOUBLE_PRECISION,
     .      0, E_OUT2, MPI_COMM_WORLD, request, ierr)
      call mpi_wait(request, status, ierr)     
      end
