# Is there an MPI All Gather operation for matrices?

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.

• How is the matrix stored? Sparse or dense? How did you end up with a system distributed like this and what are you going to do with it once you have a redundant copy everywhere? – Jed Brown Mar 9 '12 at 13:55
• it is related to a previous question I asked: scicomp.stackexchange.com/questions/1283/…. – Paul Mar 13 '12 at 3:26
• @JedBrown: Let's say I have 100 large linear systems on each of 10 processors (so 1000 systems in total). The solution of each linear system produces 1 value of interest, so each processor ends up with 100 values of interest. These 'values of interest' form a distributed vector, which I am using an all-gather operation so that each processor gets the entire vector. Each processor is then going to assemble and solve the same linear system (redundant work) to avoid communication costs. The algorithm iteratively repeats this process. – Paul Mar 13 '12 at 3:32
• From this description, it sounds like you have no use for "gathering" distributed sparse matrices, you just need to gather the small (dense) vector. – Jed Brown Mar 13 '12 at 15:03
• For a 1D problem, I'm only gathering a 1D vector. I'm trying to apply the code to a 2D problem, which would require a 2D matrix (at least, that how I want to structure the decomposition of the problem). – Paul Mar 13 '12 at 23:24

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    = {globalnrows, globalncols};
int subsizes = {globalnrows, 1};
int starts   = {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[start]), ncols, column,
&(matrix), 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);
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

• How different would things be if the matrix were sparse, symmetric, and structured? – Paul Mar 13 '12 at 3:37
• It all comes down to data layout. How is the data structured? If it's a (say) tridiagonal matrix, you could do the above with nrows=2 and be just about there. More complex structure would result in more complex code.. – user389 Mar 13 '12 at 11:51

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.

• How do you mean by a more 'flexible' approach? – Paul Mar 13 '12 at 3:37