! parallellife.f

! © Copyright 2005 Dauger Research, Inc. 
    
! Our lawyers made us say this: 
! DISCLAIMER: We provide the following on an "AS IS" basis. Use it at your own risk.  

      module lifelib
! MPI library 
      include 'mpif.h'
      contains 

      subroutine fractionalsleep(sleepRequest)
! like sleep(), except that it accepts fractions of a second 
      integer, dimension(2) :: request, returned
      real :: sleepRequest 
      if (sleepRequest>0) then
         ! <time.h>
         request(2)=int(1e9*modulo(sleepRequest, 1.0))
         request(1)=int(sleepRequest) 
         call nanosleep(request, returned)
      end if
      end subroutine


      subroutine propogatelife(out, in) ! assumed to be equal in size 
! Propogates life according to rules by J. Conway in 1970 
      implicit none
      integer, dimension(0:,0:) :: out, in
      
      integer :: row, column, neighborCount 
      
      if ((size(out).ge.9).and.(size(in).ge.9)) then
        do row=1, ubound(in,2)-1
			
            do column=1, ubound(in,1)-1
              neighborCount = 0
              
              if (in(column-1,row).gt.0) neighborCount = neighborCount + 1
              if (in(column+1,row).gt.0) neighborCount = neighborCount + 1
              
              if (in(column-1,row+1).gt.0) neighborCount = neighborCount + 1
              if (in(column,row+1).gt.0) neighborCount = neighborCount + 1
              if (in(column+1,row+1).gt.0) neighborCount = neighborCount + 1
              
              if (in(column-1,row-1).gt.0) neighborCount = neighborCount + 1
              if (in(column,row-1).gt.0) neighborCount = neighborCount + 1
              if (in(column+1,row-1).gt.0) neighborCount = neighborCount + 1
				
				select case (neighborCount) 
					case default 
						out(column,row) = 0 ! die 
					case (0:1) 
						! not enough 
						out(column,row) = 0 ! die 
					case (4:9)
						! too much 
						out(column,row) = 0 ! die 
					case (2)
                      if (in(column,row).eq.0) then
                        out(column,row) = 0 ! stay dead 
                      else
                        out(column,row) = in(column,row) + 1  ! live 
                        if (out(column,row).gt.1000000) out(column,row)=1000000 
                      end if
					case (3) ! just right 
                      out(column,row) = in(column,row) + 1  ! live 
                      if (out(column,row).gt.1000000) out(column,row)=1000000 
				end select
				
			end do
			
        end do
		
      end if
	
      end subroutine

      subroutine maintainboundaryconditions(in, idproc, nproc) 
      ! Maintains periodic boundary conditions 
      implicit none
      integer, dimension(0:,0:) :: in
      integer :: idproc, nproc
      
      integer :: row, column
      integer :: i, ierr, leftReq, rightReq, rightIDProc, leftIDProc
      integer, dimension(MPI_STATUS_SIZE) :: status 
      
      if (size(in).ge.9) then
        
        do row=1, ubound(in,2)-1
		! copy ends of rows 
			in(0,row) = in(ubound(in,1)-1, row)
			in(ubound(in,1),row) = in(1, row)
        end do
		
        if (nproc.gt.0) then

            rightIDProc=idproc+1 
            leftIDProc=idproc-1 
			
			if (rightIDProc>=nproc) rightIDProc=0 
            if (leftIDProc<0) leftIDProc=nproc-1 
			
            do i=1,2
              if (iand(ieor(i,idproc),1).ne.0) then
				call MPI_IRECV(in(0,0), ubound(in,1)+1, MPI_INTEGER, leftIDProc, leftIDProc, MPI_COMM_WORLD, leftReq, ierr)
				call MPI_SEND(in(0,1), ubound(in,1)+1, MPI_INTEGER, leftIDProc, idproc, MPI_COMM_WORLD, ierr) 
				call MPI_WAIT(leftReq, status, ierr)
              else 
				call MPI_IRECV(in(0,ubound(in,2)), ubound(in,1)+1, MPI_INTEGER, rightIDProc, rightIDProc, MPI_COMM_WORLD, rightReq, ierr)
				call MPI_SEND(in(0,ubound(in,2)-1), ubound(in,1)+1, MPI_INTEGER, rightIDProc, idproc, MPI_COMM_WORLD, ierr)
				call MPI_WAIT(rightReq, status, ierr)
              end if
            end do
			
		else
		! copy first and last rows 
          do column=0, ubound(in,1)
            in(column, 0) = in(column, ubound(in,2)-1)
            in(column, ubound(in,2)) = in(column, 1)
  		  end do
		end if
		
      end if
      
      end subroutine


      integer function myrandom() 
      implicit none
      real :: r
      
      r = rand()*2000000.
      
      if (r<0) r = - r
      
      myrandom = int(r)
      
      end function 


      integer function printresult(in) 
      implicit none
      integer, dimension(0:,0:) :: in
      
      integer :: row, column, c
      character(len=ubound(in,2)) :: line

      printresult = 0
      if (size(in).gt.9) then ! print the array 
        
        do column = 1, ubound(in,1)-1
          
          line = ""
          
          do row = 1, ubound(in,2)-1
            
            select case (in(column, row))
              case (0) 
                c = 32  
              case (2:10) 
                printresult = printresult + 2 
                c = in(column, row) + 48 -1 
              case (1)
                printresult = printresult + 1
                c = in(column, row) + 48 -1 
              case default
                c = 42
            end select 
            
            line(row:row) = char(c)
            
          end do
          line(ubound(in,2):ubound(in,2)) = char(0)
          
          print *, line 
        end do
		
      end if

      printresult = printresult/2
      end function


      subroutine addmaterial(in)
! when called add material to the life array 
      implicit none
      integer, dimension(0:,0:) :: in
      
      integer, dimension(9) :: genArray=(/z'0247', z'07d9', z'09be', z'0fb9', z'17d9', z'1f4e', z'1fdc', z'27df', z'8957'/)
      integer :: rowstart, colstart, s, i, j
      
      if (size(in).ge.9) then
        rowstart=modulo(myrandom(),(ubound(in,2)-7)) 
        colstart=modulo(myrandom(),(ubound(in,1)-7)) 
		if (rowstart<0) rowstart=-rowstart
		if (colstart<0) colstart=-colstart 
		rowstart = rowstart + 1 
		colstart = colstart + 1 
		
        if (iand(myrandom(),3).ne.0) then 
          s = 1
        else
          s=modulo(myrandom(),size(genArray))+1
        end if
        if (s<1) s=1
		
		print *, s
		
		s=genArray(s)
		
		print "('adding material ',z4)", s
		call fractionalsleep(1.0)
		
		if (iand(myrandom(),1).ne.0) then ! flip 
           s=ishft(iand(s,z'0f000'),-12)+ &
                ishft(iand(s,z'0f00'),-4)+ &
                ishft(iand(s,z'0f0'),4)+ &
                ishft(iand(s,z'0f'),12)
		end if
		
		if (iand(myrandom(),1).ne.0) then ! flip  
             s=ishft(iand(s,z'08888'),-3)+ &
                ishft(iand(s,z'04444'),-1)+ &
                ishft(iand(s,z'02222'),1)+ &
                ishft(iand(s,z'01111'),3)
		end if
		
        do j=1, 4
          do i=1, 4
            in(i+colstart,rowstart+j) = in(i+colstart,rowstart+j) + modulo(s,2)*4 
            s=s/2 
          end do
        end do
		
      end if
      
      end subroutine 


      subroutine performcalculation()
      implicit none
      
      integer, parameter :: RowDimension=78, ColumnDimension=22
! allocate 
      integer, dimension(0:ColumnDimension+1,0:RowDimension+1), target :: arrayA, arrayB
      
      integer :: frameCount, countdown, row, column, activecellcount
      integer, dimension(:,:), pointer :: lastArray, nextArray
      integer :: nproc, idproc, ierr
      
      countdown = 8
		
! set initial conditions  
      do row = 1, RowDimension
        do column = 1, ColumnDimension
          if (modulo(myrandom()/4,2).eq.0) then
            arrayA(column, row) = 1
          else
            arrayA(column, row) = 0
          end if
        end do
      end do
      
      frameCount = printresult(arrayA) 
		
! initialize parallel processing 
      ierr = ppinit(idproc, nproc)
	
      if (ierr.ne.0) then ! stop right there if there's a problem 
	    print *, "MPI_INIT has returned error ",ierr,".  Assuming single-processor... " 
        nproc = 0
        idproc = 0
      end if
	
      do frameCount = 1, 5000
        if (modulo(frameCount,2).eq.0) then
          lastArray => arrayB
          nextArray => arrayA
        else
          lastArray => arrayA
          nextArray => arrayB
        end if

        call maintainboundaryconditions(lastArray, idproc, nproc)
        call propogatelife(nextArray, lastArray)
			
        call system("clear")
			
        print *, "Frame ",frameCount," idproc = ", idproc 
        activecellcount = printresult(nextArray) 
			
        call fractionalsleep(0.125)
			
        if (activecellcount*100<(ColumnDimension*RowDimension)) then
          if (countdown.eq.0) then
            call addmaterial(nextArray)
            countdown = 16
          else 
            countdown = countdown - 1 
          end if
        else 
          countdown = 32 
        end if
		
      end do
	
      if (nproc.gt.0) call ppexit() 
     
      end subroutine
      
      function ppinit(idproc, nproc)
      implicit none
! this subroutine initializes parallel processing
      integer :: idproc ! = processor id
      integer :: nproc ! = number of real or virtual processors obtained 
      integer :: ierr, ppinit
      
      nproc=0 

! initialize the MPI execution environment 

      call MPI_INIT(ierr)
      if (ierr.eq.0) then

	! determine the rank of the calling process in the communicator 
        call MPI_COMM_RANK(MPI_COMM_WORLD, idproc, ierr)

	! determine the size of the group associated with a communicator 
        call MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr)

      end if
	
      ppinit = ierr
      
      end function

      subroutine ppexit()
      implicit none
! this subroutine terminates parallel processing 
      integer :: ierr

! terminate MPI execution environment 
      call MPI_Finalize(ierr)
	
      end subroutine
    	
      end module
      
      program life
      use lifelib
      implicit none

      integer, dimension(4) :: throwoutcount
      integer :: i
	
      call itime_(throwoutcount)
      throwoutcount(1) = iand(throwoutcount(1)+throwoutcount(2)+throwoutcount(3),z'03ff') 
      print *, "Initializing random number generator by throwing out ",throwoutcount(1)," random numbers..."
      do i=1,throwoutcount(1)
       throwoutcount(2)=myrandom()
      end do
      
      print *, "Beginning calculation... (first random is ",myrandom(),")"
	
      call performcalculation

! clean up after myself 
	
      call ppexit()
	
      stop
      end program