Affecting parts of an array outside MPI process efficiently

Question

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.

Answer 1

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.

Answer 2

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).