PROGRAM life
      IMPLICIT NONE
      INCLUDE "mpif.h"
!
!     MAXLOOP iterations of an NxN life board
!
      INTEGER :: N, MAXLOOP
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: step
      INTEGER loop, i, j, cr, nx, r, c, nb, ls
      CHARACTER*(16) name
	  REAL(8) :: wall1, wall2, wall
!    
!     MPI variables
	  INTEGER :: my_rank, numprocs, ierr, tag
	  INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status
	  CHARACTER(MPI_MAX_PROCESSOR_NAME) :: server
!  
	  INTEGER :: p, my_start, my_end, my_out_col0, my_out_col1, my_cols, loc_out_col0, loc_out_col1
	  INTEGER :: my_inn_col0, my_inn_col1, loc_inn_col0, loc_inn_col1, istart, iend
!
!     MPI initialization
!
	  call MPI_INIT( ierr )
	  call MPI_COMM_RANK( MPI_COMM_WORLD, my_rank, ierr )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
	  call MPI_Get_processor_name(server,ls,ierr)
	  WRITE(*,*) " Activated process ",my_rank," of ",numprocs, &
	  &          " on server ",TRIM(server)
!
!     Initialization
!
	  IF ( my_rank == 0 ) THEN
	  	  WRITE(*,*) " N, MaxLoop ?"
	  	  READ*,N,MAXLOOP
	  	  ! Allocate full step board
	  	  ALLOCATE(step(N,N,2), STAT=ierr)
	  	  IF ( ierr /= 0 ) THEN
	  	  	  WRITE(*,*) " Error allocating step(N,N,2)"
	  	  	  ierr = -1
	  	  	  CALL Terminate(my_rank,ierr)
	  	  	  STOP
	  	  END IF
!         Initialize step
          step = 0
!         Ten active cells in a line
          step(N/2,((N/2)-5):((N/2)+4),1) = 1
!
!         Print starting config to file Life0000.pgm
!
	      loop = 0
	      WRITE(name,"(A4,I4.4)") "Life",loop
          CALL INTDATA2MPGM(N,N,step(:,:,1),name,4)
      END IF
      !  Send initial data to processes
      CALL MPI_BCAST(N,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      CALL MPI_BCAST(MAXLOOP,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      IF ( my_rank == 0 ) THEN
         CALL ProcessBand(N,my_start,my_end,numprocs,my_rank)    !  Global columns
	  my_out_col0 = MAX(my_start-1,1)   ! First allocated column in global numeration
          my_out_col1 = MIN(my_end+1,N)     ! Last allocated column in global numeration
          my_inn_col0 = my_out_col0 + 1
          my_inn_col1 = my_out_col1 - 1
          my_cols = (my_out_col1-my_out_col0+1) ! Total allocated columns
         DO p = 1, numprocs-1
         !  Send process data
            CALL ProcessBand(N,istart,iend,numprocs,p)    !  Global columns
            tag = 100   ! istart
            call MPI_SEND(istart,1,MPI_INTEGER, &
              &     p,tag,MPI_COMM_WORLD,ierr)      
            IF ( ierr /= 0 ) THEN
	  	  	   WRITE(*,*) " Error MPI_SEND(istart,  from 0 to ",p
	  	  	   ierr = -1
	  	  	   CALL Terminate(my_rank,ierr)
	  	  	   STOP
	     END IF

            tag = 110   ! iend
            call MPI_SEND(iend,1,MPI_INTEGER, &
              &     p,tag,MPI_COMM_WORLD,ierr)      
            IF ( ierr /= 0 ) THEN
	  	  	   WRITE(*,*) " Error MPI_SEND(iend,  from 0 to ",p
	  	  	   ierr = -1
	  	  	   CALL Terminate(my_rank,ierr)
	  	  	   STOP
	    END IF

	  	    tag = 120   ! step board
            call MPI_SEND(step(1,istart,1),(iend-istart+1)*N,MPI_INTEGER, &
              &     p,tag,MPI_COMM_WORLD,ierr)      
            IF ( ierr /= 0 ) THEN
	  	  	   WRITE(*,*) " Error MPI_SEND(step(1,istart,1) from 0 to ",p
	  	  	   ierr = -1
	  	  	   CALL Terminate(my_rank,ierr)
	  	  	   STOP
	    END IF
         END DO
      ELSE   ! my_rank > 0   
      	  tag = 100   ! istart
          CALL MPI_RECV(my_start,1, &
          &             MPI_INTEGER, 0, tag, MPI_COMM_WORLD, status, ierr)   
          IF ( ierr /= 0 ) THEN
	  	  	   WRITE(*,*) " Error MPI_RECV(istart,1, from 0 in process ",my_rank
	  	  	   ierr = -1
	  	  	   CALL Terminate(my_rank,ierr)
	  	  	   STOP
	  END IF

	  	  tag = 110   ! iend
          CALL MPI_RECV(my_end,1, &
          &             MPI_INTEGER, 0, tag, MPI_COMM_WORLD, status, ierr)     	  
          IF ( ierr /= 0 ) THEN
	  	  	   WRITE(*,*) " Error MPI_RECV(iend,1, from 0 in process ",my_rank
	  	  	   ierr = -1
	  	  	   CALL Terminate(my_rank,ierr)
	  	  	   STOP
	  END IF

	  my_out_col0 = MAX(my_start-1,1)   ! First allocated column in global numeration
          my_out_col1 = MIN(my_end+1,N)     ! Last allocated column in global numeration
          my_inn_col0 = my_out_col0 + 1
          my_inn_col1 = my_out_col1 - 1
          my_cols = (my_out_col1-my_out_col0+1) ! Total allocated columns
      	  ALLOCATE(step(N,my_cols,2), STAT=ierr)
	  	  IF ( ierr /= 0 ) THEN
	  	  	  WRITE(*,*) " Error allocating step(N,my_cols,2) in process ",my_rank
	  	  	  ierr = -1
	  	  	  CALL Terminate(my_rank,ierr)
	  	  	  STOP
	  	  END IF
	  	  step = 0   ! board initialization
	  	  tag = 120   ! step board
          CALL MPI_RECV( step(1,(my_start-my_out_col0+1),1),(my_end-my_start+1)*N, &
          &             MPI_INTEGER, 0, tag, MPI_COMM_WORLD, status, ierr)
          IF ( ierr /= 0 ) THEN
	  	  	   WRITE(*,*) " Error MPI_RECV( step(1,(my_start-my_out_col0+1),1) &
	  	  	   & from 0 in process ",my_rank
	  	  	   ierr = -1
	  	  	   CALL Terminate(my_rank,ierr)
	  	  	   STOP
	  	    END IF
      END IF
      loc_out_col0 = 1   ! First allocated column in local numeration
      loc_inn_col0 = 2   ! First column to be computed in local numeration
      loc_out_col1 = loc_out_col0 + my_out_col1 - my_out_col0
      loc_inn_col1 = loc_inn_col0 + my_inn_col1 - my_inn_col0   ! Last column to be computed in local numeration
!
!     Perform MAXLOOP updates
!  

      DO loop = 1, MAXLOOP
         cr = 2-mod(loop,2)
         nx = 2-mod(loop+1,2)
         ! Send border values
         IF ( my_end < N ) then   
 	  	    tag = 310   ! Send my rightmost column to left border of receiver
            CALL MPI_SEND( step(1,loc_inn_col1,cr),N, &
            &             MPI_INTEGER, my_rank+1, tag, MPI_COMM_WORLD, ierr)
            IF ( ierr /= 0 ) THEN
	  	  	     WRITE(*,*) " Error MPI_SEND( step(1,loc_inn_col1+1,1),N, &
	  	  	        & in process ",my_rank
	  	  	     ierr = -1
	  	  	     CALL Terminate(my_rank,ierr)
	  	  	     STOP
	    END IF
         END IF  
         IF ( my_start > 1 ) then   
 	  	    tag = 320   ! Send my leftmost column to right border of receiver
            CALL MPI_SEND( step(1,loc_inn_col0,cr),N, &
            &             MPI_INTEGER, my_rank-1, tag, MPI_COMM_WORLD, ierr)
            IF ( ierr /= 0 ) THEN
	  	  	     WRITE(*,*) " Error MPI_SEND( step(1,loc_inn_col0-1,1),N, &
	  	  	        &  in process ",my_rank
	  	  	     ierr = -1
	  	  	     CALL Terminate(my_rank,ierr)
	  	  	     STOP
	    END IF
         END IF
         ! Acquire border values
         IF ( my_start > 1 ) then   
 	  	    tag = 310   ! Acquire left border
            CALL MPI_RECV( step(1,loc_out_col0,cr),N, &
            &             MPI_INTEGER, my_rank-1, tag, MPI_COMM_WORLD, status, ierr)
            IF ( ierr /= 0 ) THEN
	  	  	     WRITE(*,*) " Error MPI_RECV( step(1,loc_out_col0,cr),N, &
	  	  	        & in process ",my_rank
	  	  	     ierr = -1
	  	  	     CALL Terminate(my_rank,ierr)
	  	  	     STOP
	    END IF
         END IF  
         IF ( my_end < N ) then   
 	  	    tag = 320   ! Acquire right border
            CALL MPI_RECV( step(1,loc_out_col1,cr),N, &
            &             MPI_INTEGER, my_rank+1, tag, MPI_COMM_WORLD, status, ierr)
            IF ( ierr /= 0 ) THEN
	  	  	     WRITE(*,*) " Error MPI_RECV( step(1,loc_out_col1,1),N, &
	  	  	        & in process ",my_rank
	  	  	     ierr = -1
	  	  	     CALL Terminate(my_rank,ierr)
	  	  	     STOP
	    END IF
         END IF
         DO r = 2, N-1
            DO c = loc_inn_col0, loc_inn_col1   ! local columns
               step(r,c,nx) = step(r,c,cr)   ! Inizialize cell
               !  Per each cell
               nb = 0
               DO i = r-1,r+1
                  DO j = c-1,c+1
                      ! Count active neighbours
                  	  IF ( .NOT. (i == r .AND. j == c) ) THEN
                  	  	  nb = nb + step(i,j,cr)
                  	  END IF
                  END DO
               END DO
               IF (nb.lt.2 .or. nb.gt.3) THEN
                  step(r,c,nx) = 0
               ELSE IF (nb.eq.3) THEN
                  step(r,c,nx) = 1
               END IF
            END DO
         END DO
!
!        Send data to root process
!        
	 tag = 200
	 IF ( my_rank == 0 ) THEN
	    DO p = 1, numprocs-1
               CALL ProcessBand(N,istart,iend,numprocs,p)
               call MPI_RECV(step(1,istart,nx),(iend-istart+1)*N,MPI_INTEGER, &
                 &     p,tag,MPI_COMM_WORLD,status,ierr)      
               IF ( ierr /= 0 ) THEN
	  	  WRITE(*,*) " Error MPI_RECV(step(1,istart,nx) from ",p
	  	  ierr = -1
	  	  CALL Terminate(my_rank,ierr)
	  	  STOP
	        END IF
            END DO
	 ELSE   ! Send data to root 
            call MPI_SEND(step(1,(my_start-my_out_col0+1),nx),(my_end-my_start+1)*N, MPI_INTEGER, &
              &     0,tag,MPI_COMM_WORLD,ierr)      
            IF ( ierr /= 0 ) THEN
	  	WRITE(*,*) " Error MPI_SEND(step(1,(my_start-my_out_col0+1),1),(my_end-my_start+1)  &
	  	& from process ",my_rank
	  	ierr = -1
	  	CALL Terminate(my_rank,ierr)
	  	STOP
	    END IF
	 END IF	    
         IF ( my_rank == 0 ) THEN
!
!           Write out computed state of board
!
       	    WRITE(name,"(A4,I4.4)") "Life",loop
            CALL INTDATA2MPGM(N,N,step(:,:,nx),name,4)
         END IF
      END DO

      CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
      ierr = 0
      CALL Terminate(my_rank,ierr)

   END PROGRAM life

      
subroutine IntData2Mpgm(s1,s2,idata,name,Mag)
!  Simple subroutine to dump integer data in a PGM format
   implicit none
   INTEGER, INTENT(IN) :: s1, s2, Mag
   INTEGER, DIMENSION(s1,s2), intent(in) :: idata
   character(*), intent(in) :: name
!
   integer :: i, j, ouni, r, m
   integer :: vp, vs
   real :: rmin, rmax
   character(80) :: fname
!
!  Write on unit 700 with  PGM format
   ouni = 700
   fname = TRIM(name)//".pgm"
   open(ouni,file=trim(fname),status="replace",iostat=vs)
   if ( vs /= 0 ) then
      WRITE(*,*) "!!!! Error write access to file ",trim(fname)
   end if
!  Magic code
   write(ouni,"(a)") "P2"
!  Dimensions
   write(ouni,*) s1*Mag, s2*Mag
!  Maximum value
   write(ouni,"(a)") "255"
!  Values from 0 to 255
   rmin = MINVAL(idata); rmax = MAXVAL(idata)
   vs = 0
   do i = 1, s1
      DO r = 1, Mag
      do j = 1, s2
         vp = INT ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) )
         DO m = 1, Mag
            write(ouni,"(i4)",advance="no") vp
            vs = vs + 1
         END DO
         if (vs >= 10 ) then
            write(ouni,"(a)") " "
            vs = 0
         end if
      end do
      write(ouni,"(a)") " "
      vs = 0
      END DO
   end do
   close(ouni)
   return
end subroutine IntData2Mpgm


SUBROUTINE ProcessBand(length,istart,iend,intervals,n)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: length,intervals,n
      INTEGER, INTENT(OUT) :: istart,iend
!
      INTEGER :: delta
!
      delta = ( length + (intervals - 1) ) / intervals
      istart = delta * n + 1; iend = delta * ( n + 1)
      if ( iend > length ) iend = length
!  
      RETURN
END SUBROUTINE ProcessBand
      

SUBROUTINE Terminate(p,rc)
    IMPLICIT NONE
    INCLUDE "mpif.h"
!
    INTEGER, INTENT(IN) :: P
    INTEGER, INTENT(IN OUT) :: rc
    INTEGER :: ierr
 
    WRITE(*,*) "Process ",p," terminates"
    IF ( rc < 0 ) THEN
       CALL MPI_Abort(MPI_COMM_WORLD, rc, ierr)
    ELSE
       call MPI_FINALIZE(ierr)
    END IF
 
    STOP
END SUBROUTINE Terminate