! parallelpascalstriangle.f90

! © Copyright 2003-4 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 pascallib
! MPI library 
      include 'mpif.h'
      contains 

! Routines performing elemental work. 
      subroutine addright(inputArray, inputSize)
      implicit none
      integer, dimension(:) :: inputArray
      integer :: inputSize, i 
      if (size(inputArray).gt.1) then
		do i = inputSize, 1, -1 
			inputArray(i+1)=inputArray(i+1)+inputArray(i)
        end do
      end if
      end subroutine
    
      subroutine addleft(inputArray, inputSize)
      implicit none
      integer, dimension(:) :: inputArray
      integer :: inputSize, i 

      if ((size(inputArray).gt.1)) then
        do i=1, inputSize
			inputArray(i)=inputArray(i)+inputArray(i+1)
        end do
      end if
      end subroutine

      subroutine iterationloop(theArray, arraySize, lastIteration, verbose, idproc, nproc)
      implicit none
      integer, dimension(:) :: theArray
      integer :: arraySize, lastIteration, verbose, index, i
      integer :: idproc, nproc
      integer :: leftIDProc, rightIDProc
      integer :: request, ierr
      integer, dimension(MPI_STATUS_SIZE) :: status 
	  
      if (arraySize.gt.0) then
        leftIDProc=idproc-1; rightIDProc=idproc+1
		
! get to know my (processor's) neighbors 
        if (leftIDProc.lt.0) leftIDProc=nproc-1 
        if (rightIDProc.ge.nproc) rightIDProc=0
		
! loop over indicies 
		do index = 0, lastIteration-1
! call the left or right routine for each index, reusing array 
          if (iand(1,index)/=0) then
				call addleft(theArray, arraySize)
				
				call MPI_IRECV(theArray(arraySize), 1, MPI_INTEGER, rightIDProc, index, MPI_COMM_WORLD, request, ierr)
				call MPI_SEND(theArray(1), 1, MPI_INTEGER, leftIDProc, index, MPI_COMM_WORLD, ierr) 
				call MPI_WAIT(request, status, ierr)
				
          else 
				call addright(theArray, arraySize)
				
				call MPI_IRECV(theArray(1), 1, MPI_INTEGER, leftIDProc, index, MPI_COMM_WORLD, request, ierr)
				call MPI_SEND(theArray(arraySize), 1, MPI_INTEGER, rightIDProc, index, MPI_COMM_WORLD, ierr)
				call MPI_WAIT(request, status, ierr)
          end if
			
          if (verbose/=0) then
            do i=1,arraySize
              if (iand(1,verbose)/=0) then
                write(*, '(I6)', advance="no") theArray(i)
              else 
                if (iand(1,theArray(i))/=0) then
                  write(*, '(A)', advance="no") '*'
                else
                  write(*, '(A)', advance="no") ' '
                end if
              end if
            end do
            print *," " 
          end if
          
        end do
		
      end if
	
      end subroutine
	
      subroutine computeloop(theArray, binomialExponent, idproc, nproc)
      implicit none
      integer, dimension(:) :: theArray
      integer :: binomialExponent
! local copies of the data to be output 
      integer :: arraySize, i
      integer :: idproc, nproc, ierr
! allocate an array to hold all the results 
      integer, dimension(((binomialExponent+nproc)/nproc)+1) :: myArray
      
      arraySize=((binomialExponent+nproc)/nproc)+1
	
      if (arraySize.gt.1) then
		
	! clear the array 
        do i=1, arraySize
          myArray(i)=0
        end do
		
	! seed the array 
		if (idproc.eq.((binomialExponent/2)/(arraySize-1))) then
			myArray((binomialExponent/2)-idproc*(arraySize-1)+1)=1
		end if
		
	! propagate 
        call iterationloop(myArray, arraySize, binomialExponent, 2, idproc, nproc)
		
	! proc 0's has a bigger array so we can gather all the data at the end 
		
	! gathers the data from the other arrays ... 
		call MPI_GATHER(myArray, arraySize-1, MPI_INTEGER, theArray, arraySize-1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
	! ... to proc 0
		
      end if
	
! return the array 

      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 pascal
      use pascallib
      implicit none
! BinomialExponent specifies how far to propagate 
      integer, parameter :: BinomialExponent = 30
! main copies of the sum and the array 
      integer, dimension(BinomialExponent*2) :: theArray
      integer :: ios, idproc, nproc, ierr
	
! initialize parallel processing 
      ierr = ppinit(idproc, nproc)
	
      if (ierr.ne.0) stop ! stop right there if there's a problem 
	
      print *, "I'm processor #", idproc, " in a ", nproc, "-processor cluster." 
	
      print *,"Beginning computation..." 
	
      call computeloop(theArray, BinomialExponent, idproc, nproc)
	
      if (idproc.eq.0) then! only processor 0 
		
! save the array into a data file 
        open(unit = 2, file = "output", status = "replace", iostat = ios, form = "unformatted") 
        if (ios .eq. 0) then
          print *, "writing array...\n"
	
          write(unit = 2) theArray
			
          endfile 2 
        end if
		
      end if
		
      write( *, '(I16)', advance="no") theArray(1:BinomialExponent)
      print *, " "
		
! clean up after myself 
	
      call ppexit()
	
      stop
      end program