! Implement the game of Life, as invented by J.H. Conway.  Note that
! it uses a fair number of programming techniques that may be unfamiliar
! to some people; many of them are for portability and extensibility.
! For example, the topology management is localised in Parallel and
! executed only in the root process.

! This uses process 0 as the root process, for reading from stdin;
! that usually works, but there is a later lecture on I/O.  

! Look for UseMPI to see the MPI extensions to the serial code.
! There are also a lot of auxiliary procedures used for parallel
! logic.

PROGRAM Life
    USE mpi
    IMPLICIT NONE
    INTEGER, PARAMETER :: root = 0, unit = 42, boundary = 3
    INTEGER, PARAMETER ::  byte = SELECTED_INT_KIND(2)
    LOGICAL, PARAMETER :: UseMPI = .True.
    INTEGER :: nprocs, nprocs_2, myrank, steps, bytetype, size,     &
        count, rate, error, n
    REAL(KIND=KIND(0.0D0)) :: stamp
    CHARACTER(LEN=100) :: filename
    CHARACTER(LEN=20) :: format
    INTEGER(KIND=byte), ALLOCATABLE :: total(:,:)
    CHARACTER, PARAMETER :: chars(0:2) = (/ '.', '*', '?' /)
    CHARACTER, ALLOCATABLE :: line(:)

! Initialise MPI and find out our environment.  Also create an MPI
! derived datatype for INTEGER(KIND=byte) - see the lecture on
! composite types and language standards for details.

    IF (UseMPI) THEN
        CALL MPI_Init(error)
        CALL MPI_Comm_size(MPI_COMM_WORLD,nprocs,error)
        CALL MPI_Comm_rank(MPI_COMM_WORLD,myrank,error)
        nprocs_2 = INT(SQRT(nprocs+0.0D0))
        IF (nprocs_2**2 /= nprocs) THEN
            PRINT *, 'Invalid number of processors'
            CALL MPI_Abort(MPI_COMM_WORLD,1,error)
        END IF
        CALL MPI_Type_create_f90_integer(1,bytetype,error)
        stamp = MPI_Wtime()
    ELSE
        nprocs = 1
        nprocs_2 = 1
        myrank = root
        CALL SYSTEM_CLOCK(count,rate)
        stamp = count/(rate+0.0D0)
    END IF

! Read in the parameters.

    IF (myrank == root) THEN
        READ '(A)', filename
        READ *, steps
        IF (steps < 1 .OR. steps > 1000000) THEN
            PRINT *, 'Invalid number of steps'
            IF (UseMPI) THEN
                CALL MPI_Abort(MPI_COMM_WORLD,1,error)
            ELSE
                STOP 1
            END IF
        END IF

! Read in the size and allocate the grid, including a boundary in
! parallel mode, to simplify later code.  Note that the code assumes
! that a cell is at least as wide as the boundary, for parallel mode,
! though the boundary is fixed at 1 for serial mode.

        OPEN (unit,file=filename,STATUS='OLD',ACTION='READ')
        READ (unit,*) size
        IF (size < 1 .OR. size < nprocs_2*boundary .OR.    &
                size > 1000 .OR. MOD(size,nprocs_2) /= 0) THEN
            PRINT *, 'Invalid grid size'
            IF (UseMPI) THEN
                CALL MPI_Abort(MPI_COMM_WORLD,1,error)
            ELSE
                STOP 1
            END IF
        END IF
        IF (UseMPI .AND. nprocs > 1) THEN
            ALLOCATE (total(1-boundary:size+boundary,    &
                    1-boundary:size+boundary))
        ELSE
            ALLOCATE (total(0:size+1,0:size+1))
        END IF
        total = 0

! Read in the data, and display the initial grid.  This is not very
! efficient in memory use, but so what?  For real programs, this
! often needs to be handled more subtly, which can be tricky.  Also,
! a loop may be more efficient than WHERE, but the latter is clearer.

        ALLOCATE (line(size))
        WRITE (format,'("(",I0,"A1)")') size
        DO n = 1,size
            READ (unit,format) line
            WHERE (line == '.')
                total(1:size,n) = 0
            ELSEWHERE (line == '*')
                total(1:size,n) = 1
            ELSEWHERE
                total(1:size,n) = 2
            END WHERE
            PRINT format, chars(total(1:size,n))
        END DO
        PRINT *
        DEALLOCATE(line)

! Check the input grid for validity.

        IF (MAXVAL(total) >= 2) THEN
            PRINT *, 'Invalid input grid'
            IF (UseMPI) THEN
                CALL MPI_Abort(MPI_COMM_WORLD,1,error)
            ELSE
                STOP 1
            END IF
        END IF
    END IF

! Now run the game, either in parallel or serial.  We need to pass
! down most of the variables, so we rely on host association.  We
! need to keep the boundary empty for repeatable results.

    IF (UseMPI .AND. nprocs > 1) THEN
        CALL Parallel
    ELSE
        DO n = 1,steps
            CALL Iterate(total)
            total(:0,:) = 0
            total(:,:0) = 0
            total(size+1:,:) = 0
            total(:,size+1:) = 0
        END DO
    END IF

! Print the result and time taken.

    IF (myrank == root) THEN
        DO n = 1,size
            PRINT format, chars(total(1:size,n))
        END DO
        PRINT *
        IF (UseMPI) THEN
            PRINT '(A,F0.1," seconds")', 'Time taken was about ',    &
                MPI_Wtime()-stamp
        ELSE
            CALL SYSTEM_CLOCK(count,rate)
            PRINT '(A,F0.1," seconds")', 'Time taken was about ',    &
                count/(rate+0.0D0)-stamp
        END IF
        DEALLOCATE(total)
    END IF

! Clean up and finish.

    IF (UseMPI) THEN
        CALL MPI_Finalize(error)
    END IF

CONTAINS



! This code contains all of the parallel control logic, and most of
! the MPI use.  It's out out of line purely for clarity.

    SUBROUTINE Parallel
        INTEGER :: buffer(2), size_2, extent_2, here(2) = 0,    &
            error, i, j, k
        INTEGER(KIND=byte), ALLOCATABLE :: grid(:,:),  bundle(:,:,:)
        INTEGER, ALLOCATABLE :: topology(:,:)

! Distribute the parameters and allocate the local grid.

       IF (myrank == root) THEN
            buffer(1) = steps
            buffer(2) = size
        END IF
        CALL MPI_Bcast(buffer,2,MPI_INTEGER,root,MPI_COMM_WORLD,error)
        steps = buffer(1)
        size = buffer(2)
        size_2 = size/nprocs_2
        extent_2 = size_2+2*boundary
        ALLOCATE (topology(0:nprocs_2-1,0:nprocs_2-1),    &
            grid(1-boundary:size_2+boundary,    &
                1-boundary:size_2+boundary))

! Distributing the data is non-trivial, as it needs reordering, and
! boundaries need considering.  You do not have to use modern Fortran's
! array sections, but they are useful.  Note why I added a boundary
! around the overall grid!  Again, it's not very efficient in memory
! use.

        IF (myrank == root) THEN
            ALLOCATE (bundle(1-boundary:size_2+boundary,    &
                    1-boundary:size_2+boundary,nprocs))
            k = 0
            DO j = 0,nprocs_2-1
                DO i = 0,nprocs_2-1
                    topology(i,j) = k
                    IF (k == myrank) here = (/ i, j /)
                    k = k+1
                    bundle(:,:,k) =    &
                        total(i*size_2+1-boundary:    &
                                i*size_2+size_2+boundary,    &
                            j*size_2+1-boundary:    &
                                j*size_2+size_2+boundary)
                END DO
            END DO
        ELSE
            ALLOCATE (bundle(1,1,1))
        END IF

! Now do the actual distribution, and work out where we are.

        CALL MPI_Bcast(topology,nprocs,MPI_INTEGER,root,    &
            MPI_COMM_WORLD,error)
        CALL MPI_Scatter(bundle,extent_2**2,bytetype,    &
            grid,extent_2**2,bytetype,root,MPI_COMM_WORLD,error)
        DEALLOCATE(bundle)
        DO j = 0,nprocs_2-1
            DO i = 0,nprocs_2-1
                IF (topology(i,j) == myrank) here = (/ i, j /)
            END DO
        END DO

! Run the game for the appropriate number of steps, synchronising
! appropriately, and remembering to clean up the true boundary.

        DO k = 1,steps
            CALL Iterate(grid)
            IF (here(1) <= 0) grid(:0,:) = 0
            IF (here(1) >= nprocs_2-1) grid(size_2+1:,:) = 0
            IF (here(2) <= 0) grid(:,:0) = 0
            IF (here(2) >= nprocs_2-1) grid(:,size_2+1:) = 0
            IF (MOD(k,boundary) == 0 .AND. k < steps)    &
                CALL Synchronise(grid,size_2,topology,here)
        END DO

! Collect the data together again - this is just the converse of how
! we distributed it, though we don't need (or want) the boundaries.
! Yes, it forces an array copy.  Then clean up.

        IF (myrank == root) THEN
            ALLOCATE (bundle(size_2,size_2,nprocs))
        ELSE
            ALLOCATE (bundle(1,1,1))
        END IF
        CALL MPI_Gather(    &
            grid(1:size_2,1:size_2),size_2*size_2,bytetype,    &
            bundle,size_2*size_2,bytetype,    &
            root,MPI_COMM_WORLD,error)
        IF (myrank == root) THEN
            k = 0
            DO j = 0,nprocs_2-1
                DO i = 0,nprocs_2-1
                    k = k+1
                    total(i*size_2+1:i*size_2+size_2,    &
                            j*size_2+1:j*size_2+size_2) =    &
                        bundle(:,:,k)
                END DO
            END DO
        END IF
        DEALLOCATE(topology,grid,bundle)

    END SUBROUTINE Parallel



! This synchronises the grids.  It assumes that the grid and topology
! are valid, both in shape and content. It transfers the boundaries in
! grid order, rather than proximity order.
!
! The communication arrays are rather larger than is needed, and are
! still allocated at termination, but that doesn't matter.  If it does,
! they could be made smaller (with some loss of clarity), or allocated
! and deallocated in the main program or Parallel.
!
! Note that it is both long and fairly complicated, and that is
! reflected in the time taken to debug it.  As the lecture says,
! actually calling MPI is the easy part!

    SUBROUTINE Synchronise (grid, size_2, topology, here)
        INTEGER(KIND=byte) :: grid(1-boundary:,1-boundary:)
        INTEGER, INTENT(IN) :: size_2, topology(0:,0:), here(2)
        INTEGER :: requests(16), statuses(MPI_STATUS_SIZE,16), error, i
        INTEGER, PARAMETER :: tag = 123
        INTEGER, SAVE :: edge = 0, corner = 0, procs(8) = (/ (-1, i = 1,8) /)
        INTEGER(KIND=byte) :: send_1(boundary,size_2,2),    &
            send_2(size_2,boundary,2), send_3(boundary,boundary,4),    &
            receive_1(boundary,size_2,2),    &
            receive_2(size_2,boundary,2),    &
            receive_3(boundary,boundary,4)

! Initialise the processor numbers for the edges and corners.

        IF (procs(1) < 0) THEN
            IF (here(1) > 0) procs(1) = topology(here(1)-1,here(2))
            IF (here(1) < UBOUND(topology,1))    &
                procs(2) = topology(here(1)+1,here(2))
            IF (here(2) > 0) procs(3) = topology(here(1),here(2)-1)
            IF (here(2) < UBOUND(topology,2))    &
                procs(4) = topology(here(1),here(2)+1)
            IF (here(1) > 0 .AND. here(2) > 0)    &
                procs(5) = topology(here(1)-1,here(2)-1)
            IF (here(1) < UBOUND(topology,1) .AND.    &
                    here(2) < UBOUND(topology,2))    &
                procs(6) = topology(here(1)+1,here(2)+1)
            IF (here(1) > 0 .AND. here(2) < UBOUND(topology,2))    &
                procs(7) = topology(here(1)-1,here(2)+1)
            IF (here(1) < UBOUND(topology,2) .AND. here(2) > 0)    &
                procs(8) = topology(here(1)+1,here(2)-1)
            edge = boundary*size_2
            corner = boundary*boundary
        END IF

! Initialise for transferring.  Setting the elements of the buffers
! to an invalid value isn't needed, but is useful when debugging.
! Initialising the requests IS critical.

        send_1 = -1
        send_2 = -1
        send_3 = -1
        receive_1 = -1
        receive_2 = -1
        receive_3 = -1
        requests = MPI_REQUEST_NULL
 
! Send the edges.  Note that we need to copy the data, because we need
! to access non-contiguous sections and are using non-blocking transfers.

        IF (procs(1) >= 0) THEN
            send_1(:,:,1) = grid(1:boundary,1:size_2)
            CALL MPI_Isend(send_1(:,:,1),edge,bytetype,    &
                procs(1),tag,MPI_COMM_WORLD,requests(1),error)
        END IF
        IF (procs(2) >= 0) THEN
             send_1(:,:,2) = grid(size_2-boundary+1:size_2,1:size_2)
             CALL MPI_Isend(send_1(:,:,2),edge,bytetype,    &
                procs(2),tag,MPI_COMM_WORLD,requests(2),error)
        END IF
        IF (procs(3) >= 0) THEN
            send_2(:,:,1) = grid(1:size_2,1:boundary)
            CALL MPI_Isend(send_2(:,:,1),edge,bytetype,    &
                procs(3),tag,MPI_COMM_WORLD,requests(3),error)
        END IF
        IF (procs(4) >= 0) THEN
            send_2(:,:,2) = grid(1:size_2,size_2-boundary+1:size_2)
            CALL MPI_Isend(send_2(:,:,2),edge,bytetype,    &
                procs(4),tag,MPI_COMM_WORLD,requests(4),error)
        END IF

! Send the corners.

        IF (procs(5) >= 0) THEN
            send_3(:,:,1) = grid(1:boundary,1:boundary)
            CALL MPI_Isend(send_3(:,:,1),corner,bytetype,    &
                procs(5),tag,MPI_COMM_WORLD,requests(5),error)
        END IF
        IF (procs(6) >= 0) THEN
            send_3(:,:,2) = grid(size_2-boundary+1:size_2,    &
                    size_2-boundary+1:size_2)
            CALL MPI_Isend(send_3(:,:,2),corner,bytetype,    &
                procs(6),tag,MPI_COMM_WORLD,requests(6),error)
        END IF
        IF (procs(7) >= 0) THEN
            send_3(:,:,3) = grid(1:boundary,size_2-boundary+1:size_2)
            CALL MPI_Isend(send_3(:,:,3),corner,bytetype,    &
                procs(7),tag,MPI_COMM_WORLD,requests(7),error)
        END IF
        IF (procs(8) >= 0) THEN
            send_3(:,:,4) = grid(size_2-boundary+1:size_2,1:boundary)
            CALL MPI_Isend(send_3(:,:,4),corner,bytetype,    &
                procs(8),tag,MPI_COMM_WORLD,requests(8),error)
        END IF

! Receive the data that has been sent.

        DO i = 1,2
            IF (procs(i) >= 0)    &
                CALL MPI_Irecv(receive_1(:,:,i),edge,bytetype,    &
                    procs(i),tag,MPI_COMM_WORLD,requests(i+8),error)
        END DO
        DO i = 3,4
            IF (procs(i) >= 0)    &
                CALL MPI_Irecv(receive_2(:,:,i-2),edge,bytetype,    &
                    procs(i),tag,MPI_COMM_WORLD,requests(i+8),error)
        END DO
        DO i = 5,8
            IF (procs(i) >= 0)    &
                CALL MPI_Irecv(receive_3(:,:,i-4),corner,bytetype,    &
                    procs(i),tag,MPI_COMM_WORLD,requests(i+8),error)
        END DO

! Wait for the transfers to finish and then extract the data.

        CALL MPI_Waitall(16,requests,statuses,error)
        IF (procs(1) >= 0) grid(:0,1:size_2) = receive_1(:,:,1)
        IF (procs(2) >= 0) grid(size_2+1:,1:size_2) = receive_1(:,:,2)
        IF (procs(3) >= 0) grid(1:size_2,:0) = receive_2(:,:,1)
        IF (procs(4) >= 0) grid(1:size_2,size_2+1:) = receive_2(:,:,2)
        IF (procs(5) >= 0) grid(:0,:0) = receive_3(:,:,1)
        IF (procs(6) >= 0) grid(size_2+1:,size_2+1:) = receive_3(:,:,2)
        IF (procs(7) >= 0) grid(:0,size_2+1:) = receive_3(:,:,3)
        IF (procs(8) >= 0) grid(size_2+1:,:0) = receive_3(:,:,4)

    END SUBROUTINE Synchronise 



! This performs one step in the game on the local grid.  It assumes
! that the grid is valid, both in shape and content.  It would be
! possible to do this more subtly using WHERE, but there is little
! point.  One of the higher bits is used to store what the value will
! be next step (note that the neighbour count can reach 8).  Using a
! second grid would be needed when using a bitmap and provide the
! compiler more scope for autoparallelisation, but that's an issue for
! an OpenMP course, not a MPI one!

    SUBROUTINE Iterate (grid)
        INTEGER(KIND=byte) :: grid(:,:)
        INTEGER :: i, j, k, n

        n = UBOUND(grid,1)
        DO j = 2,n-1
            DO i = 2,n-1
                k = grid(i-1,j-1)+grid(i-1,j+1)+grid(i+1,j-1)+    &
                        grid(i+1,j+1)+grid(i-1,j)+grid(i+1,j)+    &
                        grid(i,j-1)+grid(i,j+1)
                SELECT CASE (IAND(k,15))
CASE (2);             grid(i,j) = 17*grid(i,j)
CASE (3);             grid(i,j) = grid(i,j)+16
                END SELECT
            END DO
        END DO
        grid = ISHFT(grid,-4)
    END SUBROUTINE Iterate

END PROGRAM Life
