Is there an MPI All Gather operation for matrices?

Question

I have a distributed matrix, in block column format. I know that I can reshape the matrix into one long vector and use an all_gatherv operation. I just wanted to avoid the trouble of having to reshape the matrix in my code. So, I was wondering if there's an mpi all gather operation so that in the end, every processor has an exact copy of the full matrix.

Answer 1

If this is a dense matrix, it's pretty straightforward; you use MPI_Type_create_subarray or something similar (you could build it yourself out of MPI_Type_vector or whatever; MPI_Type_struct() is the just about the most general possible option) to define a single column as a datatype. (In Fortran, you wouldn't even have to do that, you'd just send (nrows) values at a time, but you'd have the same issue if you wanted to send rows instead.) Then you need to resize the datatype so that they "start" and "end" at the right place, and you're ready to start allgatherv'ing in units of columns:

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

void printMatrix (float **m, int rows, int cols)
{
    for (int i = 0; i < rows; ++i) {
        printf ("%3d: ", i);
        for (int j = 0; j < cols; ++j)
            printf ("%2.0f ", m[i][j]);
        printf ("\n");
    }
}

float **allocMat (int rows, int cols)
{
    float  *data   = (float  *) malloc (rows * cols * sizeof(float));
    float **matrix = (float **) malloc (rows * sizeof(float *));
    for (int i = 0; i < rows; i++)
        matrix[i] = & (data[i * cols]);
    return matrix;
}

int main (int argc, char *argv[])
{
    int   size, rank;
    int   i, j;
    const int root = 0;
    const int globalncols = 10, globalnrows = 10;
    int ncols, start;
    int *allncols, *allstarts;
    float **matrix;
    MPI_Datatype columnunsized, column;

    MPI_Init (&argc, &argv);
    MPI_Comm_size (MPI_COMM_WORLD, &size);
    MPI_Comm_rank (MPI_COMM_WORLD, &rank);

    /* everyone's number of columns and offsets */
    allncols = malloc(rank * sizeof(int));
    allstarts= malloc(rank * sizeof(int));

    /* everyone gets a global matrix */
    matrix = allocMat(globalnrows, globalncols);

    for (i = 0; i < globalnrows; i++)
        for (j = 0; j < globalncols; j++)
            matrix[i][j] = ( i == j? 1. : 0.);

    /* rank 0 print the results */
    if (rank == 0) {
        printf("Before:\n");
        printMatrix(matrix, globalnrows, globalncols);
    }


    /* how many columns are we responsble for? */
    ncols = (globalncols + rank)/size;

    MPI_Allgather(&ncols, 1, MPI_INT, allncols, 1, MPI_INT, MPI_COMM_WORLD);

    start = 0;
    for (int i=0; i<rank; i++)
        start += allncols[i];

    MPI_Allgather(&start, 1, MPI_INT, allstarts, 1, MPI_INT, MPI_COMM_WORLD);

    /* create the data type for a column of data */
    int sizes[2]    = {globalnrows, globalncols};
    int subsizes[2] = {globalnrows, 1};
    int starts[2]   = {0,0};
    MPI_Type_create_subarray (2, sizes, subsizes, starts, MPI_ORDER_C,
                              MPI_FLOAT, &columnunsized);
    MPI_Type_create_resized (columnunsized, 0, sizeof(float), &column);
    MPI_Type_commit(&column);


    /* everyone update their columns by adding their rank to all values */

    for (int row=0; row<globalnrows; row++)
        for (int col=start; col<start+ncols; col++)
            matrix[row][col] += rank;

    /* gather the updated columns */

    MPI_Allgatherv(&(matrix[0][start]), ncols, column,
                   &(matrix[0][0]), allncols, allstarts,
                   column, MPI_COMM_WORLD);

    /* rank 0 print the results */
    if (rank == 0) {
        printf("After:\n");
        printMatrix(matrix, globalnrows, globalncols);
    }

    MPI_Type_free (&column);
    free (matrix[0]);
    free (matrix);

    MPI_Finalize();
    return 0;
}

Running:

$ mpirun -np 4 ./columns2
Before:
  0:  1  0  0  0  0  0  0  0  0  0 
  1:  0  1  0  0  0  0  0  0  0  0 
  2:  0  0  1  0  0  0  0  0  0  0 
  3:  0  0  0  1  0  0  0  0  0  0 
  4:  0  0  0  0  1  0  0  0  0  0 
  5:  0  0  0  0  0  1  0  0  0  0 
  6:  0  0  0  0  0  0  1  0  0  0 
  7:  0  0  0  0  0  0  0  1  0  0 
  8:  0  0  0  0  0  0  0  0  1  0 
  9:  0  0  0  0  0  0  0  0  0  1 
After:
  0:  1  0  1  1  2  2  2  3  3  3 
  1:  0  1  1  1  2  2  2  3  3  3 
  2:  0  0  2  1  2  2  2  3  3  3 
  3:  0  0  1  2  2  2  2  3  3  3 
  4:  0  0  1  1  3  2  2  3  3  3 
  5:  0  0  1  1  2  3  2  3  3  3 
  6:  0  0  1  1  2  2  3  3  3  3 
  7:  0  0  1  1  2  2  2  4  3  3 
  8:  0  0  1  1  2  2  2  3  4  3 
  9:  0  0  1  1  2  2  2  3  3  4 

Answer 2

There are a number of possible solutions to this, I would recommend that you commit a new MPI data type using MPI_Type_struct(). If you have multiple matrices of different sizes, you'll need to commit a new data type for each (the data type size is static) or consider a more flexible approach.