      PROGRAM ASSIGNMENT4
      INCLUDE 'mpif.h'
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! For each body in the array x, this program calculates a correlation 
! coefficient between that body and the bodies in array y.  The program
! is designed to run with Procs MPI processes, which can be altered 
! in the parameter below.  The number of bodies in the x and y arrays
! can also be altered below.  However, if bodies is set larger than
! 10,000, the provided data file will be insufficient to provide data.
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
! Setting the number of bodies and the number of MPI processes. 
! Adjust bodies to the number of bodies in the x and y arrays.
! bodies cannot be set to more than 10,000(Will run out of data).
! Adjust procs to the number of processors being used.
      Integer   	bodies, procs, n
      Parameter         (bodies=10000)
      Parameter         (procs=4)
      Parameter         (n=bodies/procs+1)

! Declaring other variables.
      Integer           status(MPI_STATUS_SIZE),left,right
      Integer           i, j, k, rank, size, ierr, remain
      Integer           iread, xwork, ywork, request
      Real              xx,yy
      Real*8            correlation(n), x(n), y(n), y2(n)
      Real*8            result, sum


! Initializing MPI and obtaining the process no. and no. of processes.
      call MPI_INIT(ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
      if (size .ne. procs) then
         if (rank .eq. 0) then
            print*,"The number of processors must match the"
            print*,"number of processors specified in the program."
            print*,"Aborting the program."
         endif
         call MPI_FINALIZE(ierr)
         stop
      endif
    
! Determining the remainder and assigning the remaining
! elements to a processor.
      xwork=bodies/size
      iread=xwork*rank
      remain=bodies-(xwork*size)
      do i=1,remain
         if (rank .eq. i-1) xwork=xwork+1
         if (rank .ge. i) iread=iread+1
      enddo
      ywork=xwork

! Initializing the x and y arrays
         x=0.0
         y=0.0

! Initializing the correlation array
         correlation=0.0

! Determining the left and right neighbors
      left=mod(rank-1+procs,procs)
      right=mod(rank+1+procs,procs)

! Reading in the x and y arrays.
      open(unit=11,file='data',status='old', &
          access='direct',recl=8)
      do i=1,xwork
         read(11, rec=i+iread) xx,yy
         x(i) = xx
         y(i) = yy
      enddo 

! Looping over the number of processors
      do k=1,size

! Sending the y array to the right neighbor.
         call MPI_Isend(y,n,MPI_DOUBLE_PRECISION,right, &
                       1,MPI_COMM_WORLD,request,ierr)

! Calculating the correlation for the currently available ys
         do i=1,xwork
            temp=x(i)*x(i)
            do j=1,ywork
               correlation(i)=correlation(i)+temp+y(j)*y(j)
            enddo
         enddo

! Waiting for the send to complete before overwriting
         call MPI_RECV(y2,n,MPI_DOUBLE_PRECISION,left, &
                     1,MPI_COMM_WORLD,status,ierr)

! Determining if the y's have a remainder and setting the
! amount of work (ywork) appropriately.
         if (y2(n) .eq. 0.0) then
            ywork=n-1
         else 
            ywork=n
         endif

! Replacing my set of ys with the new set after my set
! has been received.
         call MPI_Wait(request, status, ierr)
            y=y2
      enddo
! Summing my correlation array
      sum=0.0
      do i=1,xwork
         sum=sum+correlation(i)
      enddo

! Combining the sums of the correlation arrays from ea. proc.
      call MPI_Reduce(sum,result,1,MPI_DOUBLE_PRECISION, &
                    MPI_SUM,0,MPI_COMM_WORLD,ierr)

      call MPI_FINALIZE(ierr)

! Printing the final answer
      if (rank .eq. 0) then
         print*,"The correlation is",result
         print*,"Time used is",t2-t1
      endif
      end

