Affecting parts of an array outside MPI process efficiently

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_Send and 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.

The basic model of MPI is "two-sided communication": you have a sender who knows where to send, and a receiver who knows from where to expect something. In your description that is not the case: the sender sends to a randomly generated receiver. You could do this with one-sided communication in MPI which will be a bit of a learning curve. There the sender (more correctly: "origin") can pick any process as receiver ("target") and put data there. So you'd use MPI_Put instead of MPI_Send.

Other possibility: the processes do a wildcard receive. This has the problem that they don't know how many times they have to issue such. But you can use a non-blocking barrier to solve that. If you can pull that off you can pat yourself on the back for being totally state-of-the-art.

Or you could use a totally different paradigm, that makes your distributed memory look like shared. Unfortunately those are generally far from efficient.

EDIT another answer suggests Isend/Irecv. I'm not sure that that is going to work: the recipient is dependent on that random number. If only one point sends you're in big trouble because everyone has to listen for that message but only one actually gets it. If everyone sends you can be a little more clever by accumulating the sends and receives. This can work in principle if your application has "supersteps"; it doesn't work if the sends are also random in time. But statistically it's still possible that one process gets zero data, so it will have an outstanding Irecv that is never satisfied. (Another problem is that you don't know how much buffer space to create.) So a better solution is to use a Reduce_scatter to find out how much data you're going to get and then receiving that.

• With one-sided communication, the receiver still has to know that data was placed in a buffer or nothing useful will ever happen to the data so sent :-) Jun 14 at 16:08
• Indeed, synchronization is a big headache, but at least the receiver only needs to know that something may be coming, not from where. Jun 14 at 16:39
• PS I'm interested in how your code is a model for FCIQMC calculations. Send me a pointer? Jun 14 at 21:12
• Sorry I'm rather busy this week so I'll try out these answers a bit later; thank you. As for your question, basically the spawning step of FCIQMC is the randomised part with pgen=1/(N-1) (i.e. uniform spawning probability to all other determinants, which is inefficient but works), death/cloning step is missing from this code (it is just diagonal so can be done on each processor easily), and the annihilation step is arr = arr + arr_new. There are certainly details missing throughout but I think if I can parallelise this, I can parallelise my FCIQMC code.
– tmph
Jun 15 at 14:47
• Well, if you Put you have to make sure not everyone puts in the same location, so you need a "stack pointer" or something like that on each process. And you update that atomically. It's all a big mess..... But otherwise: yes, this can deal with arbitary amounts of data arriving, including none. Aside: If you Put there is no Recv neeeded, so I'm not sure why you mention it. Jun 15 at 14:57

One sided communication is certainly one route. Another mechanism you could use to tackle this is use non-blocking two-sided communication, things like MPI_ISend and MPI_IRecv. If you fix some MPI process, you could set it up to have two main tasks:

1. Process requests from other processes asking for array data local to this process
2. When this process gets a response for one of its requests for data, it updates the corresponding local array element

You would probably want to make your array effectively distributed and use some global ID $$g_{\text{id}}$$ that indexes elements across the whole distributed array. One thing you could do for the global ID is construct it such that you can in constant time use it to compute the rank $$r$$ the element lives on and compute the local ID $$\ell_{\text{id}}$$ of where this element is located in the array chunk living on the rank $$r$$ process. For example, you might use something like

$$g_{\text{id}} = r N + \ell_{\text{id}}$$

where, say, $$N$$ is the number of elements on any rank (assuming the same across all ranks). Then clearly given $$g_{\text{id}}$$ we can compute $$r$$ and $$\ell_{\text{id}}$$. When you randomly choose an element, you can do so by computing a random global ID $$g_{\text{id}}$$ between $$0$$ and your largest global index, which I will assume you get at initialization of your code, and use this to request the data associated with element $$g_{\text{id}}$$ from the parent rank.

As long as each process has the mechanism in place to respond to requests and handle when it gets its response, this seems like it would be relatively straight forward to get done. Note that due to the randomness, you will want to use something like MPI_Iprobe to help you figure out if you have any messages in the message queue that you should process (be it requests or responses).

• I'm adding a discussion if Isend/Irecv to my answer. Jun 15 at 14:32
• @VictorEijkhout One comment I have on what you said is a process does not need to post an irecv in advance for the unlikely chance a given process is the randomly chosen one. Instead, the given process can use iprobe to check if some process sent them a request message and then they can setup an irecv to receive it. Obviously this is not going to be as efficient as having an irecv already posted, but it simplifies the approach. One can still post an immediate irecv corresponding to a request being sent out, though, to at least benefit from the lower latency there. Thoughts? Jun 15 at 14:52
• The problem with Iprobe is that if it returns false does that mean that no message is coming, or is it only late? In other words, when do you stop probing? Jun 15 at 14:57
• @VictorEijkhout I would actually expect the probing to be a consistent part of, say, the main loop for the algorithm. So in that sense, you may not ever really stop probing. So for a single iteration of this main loop, you (1) try to probe and process anything you get from that (if it returns true) and then (2) you do an MPI_Test to see if the request objects from your request irecv calls are satisfied. The moment you get a response for some data request, you process that accordingly and then maybe send out new request(s) if applicable. Jun 15 at 15:06
• I should add that I see the probing using MPI_ANY_SOURCE and maybe MPI_ANY_TAG in this loop. Jun 15 at 15:12