This was originally asked on Stackoverflow but instead of migrating it was suggested I repost here and delete that one. I've tried to clarify based on comments there.
I'm new to MPI, but I am trying to use it to speed up some code I have. A minimal version of that code is below:
program main !! takes a list, then for each element randomly generates an index and adds !! the element to that location !! while this program is useless the basic features are the same as a Monte !! Carlo program I am writing integer, parameter :: N=5 integer, parameter :: niter=10 integer :: arr(N), arr_new(N) ! will want dp real ! dummy integer :: i,j,step do i=1,N arr(i) = i enddo do step=1,niter arr_new = 0 ! initialise to zero do i=1,N j = randint_exc(1,N,i) arr_new(j) = arr_new(j) + i enddo arr = arr + arr_new print*, "newarr", arr_new print*, "uptarr", arr enddo contains function randint_exc(a, b, exclude) result(retval) !! get random integer between a and b, but exclude arg exclude implicit none integer, parameter :: dp = kind(1.d0) integer, intent(in) :: a, b integer, intent(in) :: exclude integer :: retval real(dp) :: u call random_number(u) retval = a + floor((b-a)*u) ! randint between a and b-1 if (retval >= exclude) then retval = retval + 1 endif end function randint_exc end program main
(FWIW I am parallelising my own implementation of FCIQMC, just for fun; I know there are good programs for that out there. I thought I'd boil it down here so that you don't have to worry about details)
Basically, I have some array of values (I know its initial values), and for each element of this array, I want to randomly select another element in the array, and add the current element to it. I then do this for some fixed number of iterations. As you can see, the way I do this is to initialise a new array to zero and add values to that, then add that new array to the original. Rinse and repeat.
My attempt at parallelising it with MPI was for each process to generate its own arrays, but I am stuck at the part where it may generate elements outside its own block. I guess I would have to check which rank process index
j belongs to, then send the index along with the value to the process and also receive (an arbitrary count). I have struggled with using
MPI_Recv for this (my attempts haven't even compiled). How would I do this, and is there some more elegant/easy way of doing it? (Also about the blocking into sections; is there a built-in MPI function?) Here is my attempt with a
! TODO ??? comment where I got stuck... Otherwise, is there a way for all MPI processes to share memory to the new array, so that I can send to any arbitrary index any time?
program main use mpi implicit none ! MPI variables integer :: ierr, nproc, rank integer :: N=5 integer, parameter :: niter=10 ! variables introduced because I'm trying to move to MPI integer :: Nlocal, r integer, allocatable :: arrlocal(:), arr_newlocal(:) ! dummy integer :: i,j,step call MPI_Init(ierr) call MPI_Comm_size(MPI_COMM_WORLD, nproc, ierr) call MPI_Comm_rank(MPI_COMM_WORLD, rank, ierr) Nlocal = N/nproc if (rank == nproc-1) then ! add remaining elements to last processor's list r = modulo(N,nproc) else r = 0 endif allocate(arrlocal(Nlocal+r), arr_newlocal(Nlocal+r)) do i=1,Nlocal+r arrlocal(i) = Nlocal*rank+i enddo print*, rank, Nlocal, "array", arrlocal do step=1,niter ! NOTE you can only start the next step when all the other processes are ! done (I think), since it will depend on the new full array ! so force all the processes to reach this point call MPI_Barrier(MPI_COMM_WORLD, ierr) arr_newlocal = 0 ! initialise to zero do i=1,N ! this is the part I am most confused about parallelising j = randint_exc(1,N,i) ! NOTE N, *not* Nlocal ! TODO ??? ! j might be outside the scope of this process ! arr_newlocal(j) = arr_newlocal(j) + i enddo arrlocal = arrlocal + arr_newlocal ! print*, step, rank, "newarr", arr_newlocal ! print*, step, rank, "uptarr", arrlocal enddo call MPI_Finalize(ierr) contains function randint_exc(a, b, exclude) result(retval) !! get random integer between a and b, but exclude arg exclude implicit none integer, parameter :: dp = kind(1.d0) integer, intent(in) :: a, b integer, intent(in) :: exclude integer :: retval real(dp) :: u call random_number(u) retval = a + floor((b-a)*u) ! randint between a and b-1 if (retval >= exclude) then retval = retval + 1 endif end function randint_exc end program main
Some of the comments made it sound like this is a lot harder to parallelise than I realised, so I'd be happy to accept an answer that's basically just a reference with "use X pattern and Y method." I am interested in the intermediate steps as well and I want to go to very large arrays (so it's memory-intensive), and this is more just a programming exercise than anything else so I'd like to parallelise the algorithm itself rather than solve the same algorithm many times in parallel and get stats from there.