The integration of MPI transfers within task parallelism is done in a very natural way by the means of asynchronous interactions between the application and StarPU.  This is implemented in a separate libstarpumpi library which basically provides "StarPU" equivalents of MPI_* functions, where void * buffers are replaced with starpu_data_handle_ts, and all GPU-RAM-NIC transfers are handled efficiently by StarPU-MPI.  The user has to use the usual mpirun command of the MPI implementation to start StarPU on the different MPI nodes.

An MPI Insert Task function provides an even more seamless transition to a distributed application, by automatically issuing all required data transfers according to the task graph and an application-provided distribution.

Simple Example
==============

The flags required to compile or link against the MPI layer are then accessible with the following commands:

     $ pkg-config --cflags starpumpi-1.0  # options for the compiler
     $ pkg-config --libs starpumpi-1.0    # options for the linker

Also pass the --static option if the application is to be linked statically.

     void increment_token(void)
     {
         struct starpu_task *task = starpu_task_create();

         task->cl = &increment_cl;
         task->handles[0] = token_handle;

         starpu_task_submit(task);
     }

     int main(int argc, char **argv)
     {
         int rank, size;

         starpu_init(NULL);
         starpu_mpi_initialize_extended(&rank, &size);

         starpu_vector_data_register(&token_handle, 0, (uintptr_t)&token, 1, sizeof(unsigned));

         unsigned nloops = NITER;
         unsigned loop;

         unsigned last_loop = nloops - 1;
         unsigned last_rank = size - 1;

         for (loop = 0; loop < nloops; loop++) {
             int tag = loop*size + rank;

             if (loop == 0 && rank == 0)
             {
                 token = 0;
                 fprintf(stdout, "Start with token value %d\n", token);
             }
             else
             {
                 starpu_mpi_irecv_detached(token_handle, (rank+size-1)%size, tag,
                         MPI_COMM_WORLD, NULL, NULL);
             }

             increment_token();

             if (loop == last_loop && rank == last_rank)
             {
                 starpu_data_acquire(token_handle, STARPU_R);
                 fprintf(stdout, "Finished: token value %d\n", token);
                 starpu_data_release(token_handle);
             }
             else
             {
                 starpu_mpi_isend_detached(token_handle, (rank+1)%size, tag+1,
                         MPI_COMM_WORLD, NULL, NULL);
             }
         }

         starpu_task_wait_for_all();

         starpu_mpi_shutdown();
         starpu_shutdown();

         if (rank == last_rank)
         {
             fprintf(stderr, "[%d] token = %d == %d * %d ?\n", rank, token, nloops, size);
             STARPU_ASSERT(token == nloops*size);
         }

Exchanging User Defined Data Interface
======================================

New data interfaces defined as explained in An example of data interface can also be used within StarPU-MPI and exchanged between nodes. Two functions needs to be defined through the type struct starpu_data_interface_ops (see Data Interface API). The pack function takes a handle and returns a contiguous memory buffer along with its size where data to be conveyed to another node should be copied. The reversed operation is implemented in the unpack function which takes a contiguous memory buffer and recreates the data handle.

     static int complex_pack_data(starpu_data_handle_t handle, unsigned node, void **ptr, size_t *count)
     {
       STARPU_ASSERT(starpu_data_test_if_allocated_on_node(handle, node));

       struct starpu_complex_interface *complex_interface =
         (struct starpu_complex_interface *) starpu_data_get_interface_on_node(handle, node);

       *count = complex_get_size(handle);
       *ptr = malloc(*count);
       memcpy(*ptr, complex_interface->real, complex_interface->nx*sizeof(double));
       memcpy(*ptr+complex_interface->nx*sizeof(double), complex_interface->imaginary,
              complex_interface->nx*sizeof(double));

       return 0;
     }

     static int complex_unpack_data(starpu_data_handle_t handle, unsigned node, void *ptr, size_t count)
     {
       STARPU_ASSERT(starpu_data_test_if_allocated_on_node(handle, node));

       struct starpu_complex_interface *complex_interface =
         (struct starpu_complex_interface *)	starpu_data_get_interface_on_node(handle, node);

       memcpy(complex_interface->real, ptr, complex_interface->nx*sizeof(double));
       memcpy(complex_interface->imaginary, ptr+complex_interface->nx*sizeof(double),
              complex_interface->nx*sizeof(double));

       return 0;
     }

     static struct starpu_data_interface_ops interface_complex_ops =
     {
       ...
       .pack_data = complex_pack_data,
       .unpack_data = complex_unpack_data
     };

MPI Insert Task Utility
=======================

To save the programmer from having to explicit all communications, StarPU provides an "MPI Insert Task Utility". The principe is that the application decides a distribution of the data over the MPI nodes by allocating it and notifying StarPU of that decision, i.e. tell StarPU which MPI node "owns" which data. It also decides, for each handle, an MPI tag which will be used to exchange the content of the handle. All MPI nodes then process the whole task graph, and StarPU automatically determines which node actually execute which task, and trigger the required MPI transfers.

 -- Function: int starpu_data_set_tag (starpu_data_handle_t HANDLE, int TAG)
     Tell StarPU-MPI which MPI tag to use when exchanging the data.

 -- Function: int starpu_data_get_tag (starpu_data_handle_t HANDLE)
     Returns the MPI tag to be used when exchanging the data.

 -- Function: int starpu_data_set_rank (starpu_data_handle_t HANDLE, int RANK)
     Tell StarPU-MPI which MPI node "owns" a given data, that is, the node which will always keep an up-to-date value, and will by default execute tasks which write to it.

 -- Function: int starpu_data_get_rank (starpu_data_handle_t HANDLE)
     Returns the last value set by starpu_data_set_rank.

 -- Macro: STARPU_EXECUTE_ON_NODE
     this macro is used when calling starpu_mpi_insert_task, and must be followed by a integer value which specified the node on which to execute the codelet.

 -- Macro: STARPU_EXECUTE_ON_DATA
     this macro is used when calling starpu_mpi_insert_task, and must be followed by a data handle to specify that the node owning the given data will execute the codelet.

 -- Function: int starpu_mpi_insert_task (MPI_Comm COMM, struct starpu_codelet *CODELET, ...)
     Create and submit a task corresponding to CODELET with the following arguments.  The argument list must be zero-terminated.

     The arguments following the codelets are the same types as for the function starpu_insert_task defined in Insert Task Utility. The extra argument STARPU_EXECUTE_ON_NODE followed by an integer allows to specify the MPI node to execute the codelet. It is also possible to specify that the node owning a specific data will execute the codelet, by using STARPU_EXECUTE_ON_DATA followed by a data handle.

     The internal algorithm is as follows:

      1. Find out which MPI node is going to execute the codelet.

           1. If there is only one node owning data in W mode, it will be selected;

           2. If there is several nodes owning data in W node, the one selected will be the one having the least data in R mode so as to minimize the amount of data to be transfered;

           3. The argument STARPU_EXECUTE_ON_NODE followed by an integer can be used to specify the node;

           4. The argument STARPU_EXECUTE_ON_DATA followed by a data handle can be used to specify that the node owing the given data will execute the codelet.

      2. Send and receive data as requested. Nodes owning data which need to be read by the task are sending them to the MPI node which will execute it. The latter receives them.

      3. Execute the codelet. This is done by the MPI node selected in the 1st step of the algorithm.

      4. If several MPI nodes own data to be written to, send written data back to their owners.

     The algorithm also includes a communication cache mechanism that allows not to send data twice to the same MPI node, unless the data has been modified. The cache can be disabled (see STARPU_MPI_CACHE).

 -- Function: void starpu_mpi_get_data_on_node (MPI_Comm COMM, starpu_data_handle_t DATA_HANDLE, int NODE)
     Transfer data DATA_HANDLE to MPI node NODE, sending it from its owner if needed. At least the target node and the owner have to call the function.

Here an stencil example showing how to use starpu_mpi_insert_task. One first needs to define a distribution function which specifies the locality of the data. Note that that distribution information needs to be given to StarPU by calling starpu_data_set_rank.

     /* Returns the MPI node number where data is */
     int my_distrib(int x, int y, int nb_nodes) {
       /* Block distrib */
       return ((int)(x / sqrt(nb_nodes) + (y / sqrt(nb_nodes)) * sqrt(nb_nodes))) % nb_nodes;

       // /* Other examples useful for other kinds of computations */
       // /* / distrib */
       // return (x+y) % nb_nodes;

       // /* Block cyclic distrib */
       // unsigned side = sqrt(nb_nodes);
       // return x % side + (y % side) * size;
     }

Now the data can be registered within StarPU. Data which are not owned but will be needed for computations can be registered through the lazy allocation mechanism, i.e. with a home_node set to -1. StarPU will automatically allocate the memory when it is used for the first time.

One can note an optimization here (the else if test): we only register data which will be needed by the tasks that we will execute.

         unsigned matrix[X][Y];
         starpu_data_handle_t data_handles[X][Y];

         for(x = 0; x < X; x++) {
             for (y = 0; y < Y; y++) {
                 int mpi_rank = my_distrib(x, y, size);
                  if (mpi_rank == my_rank)
                     /* Owning data */
                     starpu_variable_data_register(&data_handles[x][y], 0,
                                                   (uintptr_t)&(matrix[x][y]), sizeof(unsigned));
                 else if (my_rank == my_distrib(x+1, y, size) || my_rank == my_distrib(x-1, y, size)
                       || my_rank == my_distrib(x, y+1, size) || my_rank == my_distrib(x, y-1, size))
                     /* I don't own that index, but will need it for my computations */
                     starpu_variable_data_register(&data_handles[x][y], -1,
                                                   (uintptr_t)NULL, sizeof(unsigned));
                 else
                     /* I know it's useless to allocate anything for this */
                     data_handles[x][y] = NULL;
                 if (data_handles[x][y])
                     starpu_data_set_rank(data_handles[x][y], mpi_rank);
             }
         }

Now starpu_mpi_insert_task() can be called for the different steps of the application.

         for(loop=0 ; loop<niter; loop++)
             for (x = 1; x < X-1; x++)
                 for (y = 1; y < Y-1; y++)
                     starpu_mpi_insert_task(MPI_COMM_WORLD, &stencil5_cl,
                                            STARPU_RW, data_handles[x][y],
                                            STARPU_R, data_handles[x-1][y],
                                            STARPU_R, data_handles[x+1][y],
                                            STARPU_R, data_handles[x][y-1],
                                            STARPU_R, data_handles[x][y+1],
                                            0);
         starpu_task_wait_for_all();

I.e. all MPI nodes process the whole task graph, but as mentioned above, for each task, only the MPI node which owns the data being written to (here, data_handles[x][y]) will actually run the task. The other MPI nodes will automatically send the required data.

This can be a concern with a growing number of nodes. To avoid this, the application can prune the task for loops according to the data distribution, so as to only submit tasks on nodes which have to care about them (either to execute them, or to send the required data).

MPI Collective Operations
=========================

 -- Function: int starpu_mpi_scatter_detached (starpu_data_handle_t *DATA_HANDLES, int COUNT, int ROOT, MPI_Comm COMM, void (*SCALLBACK)(void *), void *SARG, void (*RCALLBACK)(void *), void *RARG)
     Scatter data among processes of the communicator based on the ownership of the data. For each data of the array DATA_HANDLES, the process ROOT sends the data to the process owning this data. Processes receiving data must have valid data handles to receive them. On completion of the collective communication, the SCALLBACK function is called with the argument SARG on the process ROOT, the RCALLBACK function is called with the argument RARG on any other process.

 -- Function: int starpu_mpi_gather_detached (starpu_data_handle_t *DATA_HANDLES, int COUNT, int ROOT, MPI_Comm COMM, void (*SCALLBACK)(void *), void *SARG, void (*RCALLBACK)(void *), void *RARG)
     Gather data from the different processes of the communicator onto the process ROOT. Each process owning data handle in the array DATA_HANDLES will send them to the process ROOT. The process ROOT must have valid data handles to receive the data. On completion of the collective communication, the RCALLBACK function is called with the argument RARG on the process ROOT, the SCALLBACK function is called with the argument SARG on any other process.

     if (rank == root)
     {
         /* Allocate the vector */
         vector = malloc(nblocks * sizeof(float *));
         for(x=0 ; x<nblocks ; x++)
         {
             starpu_malloc((void **)&vector[x], block_size*sizeof(float));
         }
     }

     /* Allocate data handles and register data to StarPU */
     data_handles = malloc(nblocks*sizeof(starpu_data_handle_t *));
     for(x = 0; x < nblocks ;  x++)
     {
         int mpi_rank = my_distrib(x, nodes);
         if (rank == root) {
             starpu_vector_data_register(&data_handles[x], 0, (uintptr_t)vector[x],
                                         blocks_size, sizeof(float));
         }
         else if ((mpi_rank == rank) || ((rank == mpi_rank+1 || rank == mpi_rank-1))) {
             /* I own that index, or i will need it for my computations */
             starpu_vector_data_register(&data_handles[x], -1, (uintptr_t)NULL,
                                        block_size, sizeof(float));
         }
         else {
             /* I know it's useless to allocate anything for this */
             data_handles[x] = NULL;
         }
         if (data_handles[x]) {
             starpu_data_set_rank(data_handles[x], mpi_rank);
         }
     }

     /* Scatter the matrix among the nodes */
     starpu_mpi_scatter_detached(data_handles, nblocks, root, MPI_COMM_WORLD);

     /* Calculation */
     for(x = 0; x < nblocks ;  x++) {
         if (data_handles[x]) {
             int owner = starpu_data_get_rank(data_handles[x]);
             if (owner == rank) {
                 starpu_insert_task(&cl, STARPU_RW, data_handles[x], 0);
             }
         }
     }

     /* Gather the matrix on main node */
     starpu_mpi_gather_detached(data_handles, nblocks, 0, MPI_COMM_WORLD);