Using MPI_Gatherv() and MPI_Datatype to 'gather' dynamically allocated 2D arrays in C++

Question

I think the easiest way to describe the problem is with a simple code. On each processor I have dynamically allocated '2D arrays' (achieved via the new*[rows],new[cols] formalism, see code below for clarification). Rightly or wrongly, I'm trying to use a committed MPI_Datatype to help me do MPI_Gatherv() to gather all the arrays into a single 2D array on the root processor.

Here's the code, and below it I highlight to salient points of it (it should be very easy to understand if compiled and ran - it asks for the dimensions of the array you desire):

#include <iostream>
#include <string>
#include <cmath>
#include <cstdlib>
#include <time.h>

#include "mpi.h" 


using namespace std;

// A function that prints out the 2D arrays to the terminal.
void print_2Darray(int **array_in,int dim_rows, int dim_cols) {
    cout << endl;
    for (int i=0;i<dim_rows;i++) {
        for (int j=0;j<dim_cols;j++) {
            cout << array_in[i][j] << " ";
            if (j==(dim_cols-1)) {
                cout << endl;
            }
        }
    }
    cout << endl;
}


int main(int argc, char *argv[]) {

    MPI::Init(argc, argv);


    // Typical MPI incantations...

    int size, rank;

    size = MPI::COMM_WORLD.Get_size(); 
    rank = MPI::COMM_WORLD.Get_rank();

    cout << "size = " << size << endl;
    cout << "rank = " << rank << endl;

    sleep(1);

    // Dynamically allocate a 2D square array of user-defined size 'dim'.

    int dim;
    if (rank == 0) {
        cout << "Please enter dimensions of 2D array ( dim x dim array ): ";
        cin >> dim;
        cout << "dim = " << dim << endl;
    }   

    MPI_Bcast(&dim,1,MPI_INT,0,MPI_COMM_WORLD);

    int **array2D;
    array2D = new int*[dim];
    for (int i=0; i<dim; i++) {
        array2D[i] = new int[dim](); // the extra '()' initializes to zero.
    }

    // Fill the arrays with i*j+rank where i and j are the indices.
    for (int i=0;i<dim;i++) {
        for (int j=0;j<dim;j++) {
            array2D[i][j] = i*j + rank;
        }
    }

    // Print out the arrays.
    print_2Darray(array2D,dim,dim);

    // Commit a MPI_Datatype for these arrays.
    MPI_Datatype MPI_ARRAYROW;
    MPI_Type_contiguous(dim, MPI_INT, &MPI_ARRAYROW);
    MPI_Type_commit(&MPI_ARRAYROW);

    // Declare 'all_array2D[][]' which will contain array2D[][] from all procs.
    int **all_array2D;
    all_array2D = new int*[size*dim];
    for (int i=0; i<size*dim; i++) {
        all_array2D[i] = new int[dim]();  // the extra '()' initializes to zero.
    }

    // Print out the arrays.
    print_2Darray(all_array2D,size*dim,dim);


    // Displacement vector for MPI_Gatherv() call.
    int *displace;
    displace = (int *)calloc(size,sizeof(int));
    int *dim_list;
    dim_list = (int *)calloc(size,sizeof(int));
    int j = 0;
    for (int i=0; i<size; i++) {
        displace[i] = j;
        cout << "displace[" << i << "] = " << displace[i] << endl;
        j += dim;
        dim_list[i] = dim;
    }

    // MPI_Gatherv call.
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Gatherv(array2D,dim,MPI_ARRAYROW,all_array2D,&dim_list[rank],&displace[rank],MPI_ARRAYROW,0,MPI_COMM_WORLD);

    // Print out the arrays.
    print_2Darray(all_array2D,size*dim,dim);

    MPI::Finalize();

    return 0;
}

The code compiles, but runs into segmentation faults (I compile with 'mpic++' and used 'mpirun -np 2' to use 2 processors):

[unknown-78-ca-39-b4-09-4f:02306] *** Process received signal ***
[unknown-78-ca-39-b4-09-4f:02306] Signal: Segmentation fault (11)
[unknown-78-ca-39-b4-09-4f:02306] Signal code: Address not mapped (1)
[unknown-78-ca-39-b4-09-4f:02306] Failing at address: 0x0
[unknown-78-ca-39-b4-09-4f:02306] [ 0] 2   libSystem.B.dylib                   0x00007fff844021ba _sigtramp + 26
[unknown-78-ca-39-b4-09-4f:02306] [ 1] 3   ???                                 0x0000000000000001 0x0 + 1
[unknown-78-ca-39-b4-09-4f:02306] [ 2] 4   gatherv2Darrays.x                   0x00000001000010c2 main + 1106
[unknown-78-ca-39-b4-09-4f:02306] [ 3] 5   gatherv2Darrays.x                   0x0000000100000a98 start + 52
[unknown-78-ca-39-b4-09-4f:02306] *** End of error message ***
mpirun noticed that job rank 0 with PID 2306 on node unknown-78-ca-39-b4-09-4f.home exited on signal 11 (Segmentation fault). 
1 additional process aborted (not shown)

The segmentation fault occurs upon execution of the 'print_2Darray(all_array2D,size*dim,dim)' function near the end of the code, where 'all_array2D' is 'supposed to' contain the gathered arrays. More specifically, the code seems to print the 'all_array2D' OK for the bit gathered from the master processor, but then gives the seg fault when the print_2Darray() function starts working on the bits from other processors.

Salient points of code:

1. I declare an MPI_Datatype that is a contiguous block of memory of sufficient size to store a single row of the 2D arrays. I then use MPI_Gatherv() to try and gathers these rows.
2. The code's sleep(1) call is just to help the user see the prompt for 'dims' more clearly, otherwise it get's buried between the 'size' and 'rank' couts.
3. The elements of the 2D array are initialized to values "i*j + rank" where i and j are the row and column indices respectively. My rationale is that the resulting numbers easily give away the rank of the processor that generated that array.

I guess it boils down to me not knowing how properly to MPI_Gatherv() dynamically allocated arrays... Should I be using MPI_Datatypes at all? It's quite important to me that the arrays are dynamically allocated.

I will be very grateful for any help/suggestions! I'm pretty much depleted of ideas!

Answer 1

I just wanted to summarise the solution which @Hristolliev and @JonathanDursi helped me get to.

1. MPI commands like MPI_Gatherv() work with contiguously allocated blocks of memory, hence use of 'new' to construct 2D arrays which then feed into MPI commands won't work since 'new' doesn't guarantee contiguous blocks. Use instead 'calloc' to make these arrays (see code below as an example).
2. An important point by @Hristolliev: The 1st and 4th arguments of MPI_Gatherv() must be pointers to the first elements of type MPI_ARRAYROW. Dereferencing the 2D arrays by one level e.g. array2D[0] will achieve this (again, see modified working code below).

The final, working code is given below:

#include <iostream>
#include <string>
#include <cmath>
#include <cstdlib>
#include <time.h>

#include "mpi.h" 


using namespace std;


void print_2Darray(int **array_in,int dim_rows, int dim_cols) {
    cout << endl;
    for (int i=0;i<dim_rows;i++) {
        for (int j=0;j<dim_cols;j++) {
            cout << array_in[i][j] << " ";
            if (j==(dim_cols-1)) {
                cout << endl;
            }
        }
    }
    cout << endl;
}


int main(int argc, char *argv[]) {

    MPI::Init(argc, argv);


    // Typical MPI incantations...

    int size, rank;

    size = MPI::COMM_WORLD.Get_size(); 
    rank = MPI::COMM_WORLD.Get_rank();

    cout << "size = " << size << endl;
    cout << "rank = " << rank << endl;

    sleep(1);

    // Dynamically allocate a 2D square array of user-defined size 'dim'.

    int dim;
    if (rank == 0) {
        cout << "Please enter dimensions of 2D array ( dim x dim array ): ";
        cin >> dim;
        cout << "dim = " << dim << endl;
    }   

    MPI_Bcast(&dim,1,MPI_INT,0,MPI_COMM_WORLD);


    // Use another way of declaring the 2D array which ensures it is contiguous in memory.
    int **array2D;
    array2D = (int **) calloc(dim,sizeof(int *));
    array2D[0] = (int *) calloc(dim*dim,sizeof(int));
    for (int i=1;i<dim;i++) {
        array2D[i] = array2D[0] + i*dim;
    }

    // Fill the arrays with i*j+rank where i and j are the indices.
    for (int i=0;i<dim;i++) {
        for (int j=0;j<dim;j++) {
            array2D[i][j] = i*j + rank;
        }
    }

    // Print out the arrays.
    print_2Darray(array2D,dim,dim);

    // Commit a MPI_Datatype for these arrays.
    MPI_Datatype MPI_ARRAYROW;
    MPI_Type_contiguous(dim, MPI_INT, &MPI_ARRAYROW);
    MPI_Type_commit(&MPI_ARRAYROW);


    // Use another way of declaring the 2D array which ensures it is contiguous in memory.
    int **all_array2D;
    all_array2D = (int **) calloc(size*dim,sizeof(int *));
    all_array2D[0] = (int *) calloc(dim*dim,sizeof(int));
    for (int i=1;i<size*dim;i++) {
        all_array2D[i] = all_array2D[0] + i*dim;
    }

    // Print out the arrays.
    print_2Darray(all_array2D,size*dim,dim);


    // Displacement vector for MPI_Gatherv() call.
    int *displace;
    displace = (int *)calloc(size,sizeof(int));
    int *dim_list;
    dim_list = (int *)calloc(size,sizeof(int));
    int j = 0;
    for (int i=0; i<size; i++) {
        displace[i] = j;
        cout << "displace[" << i << "] = " << displace[i] << endl;
        j += dim;
        dim_list[i] = dim;
        cout << "dim_list[" << i << "] = " << dim_list[i] << endl;
    }

    // MPI_Gatherv call.
    MPI_Barrier(MPI_COMM_WORLD);
    cout << "array2D[0] = " << array2D[0] << endl;
    MPI_Gatherv(array2D[0],dim,MPI_ARRAYROW,all_array2D[0],&dim_list[rank],&displace[rank],MPI_ARRAYROW,0,MPI_COMM_WORLD);


    // Print out the arrays.
    print_2Darray(all_array2D,size*dim,dim);


    MPI::Finalize();

    return 0;
}

Compile with mpic++.

Answer 2

MPI_Gatherv, MPI_Scatterv, and in fact all other MPI communication calls that take array arguments, expect that array elements are laid out consecutively in memory. This means that in the call MPI_Gatherv(array2D, dim, MPI_ARRAYROW, ...), MPI expects that the first element of type MPI_ARRAYROW starts at the memory location that array2D points to, the second element starts at (BYTE*)array2D + extent_of(MPI_ARRAYROW), the third element starts at (BYTE*)array2D + 2*extent_of(MPI_ARRAYROW), and so on. Here extent_of() is the extent of the MPI_ARRAYROW type, which can be obtained by calling MPI_Type_get_extent.

Clearly the rows of your 2D array are not consecutive in memory since each of them is allocated by a separate invocation of the new operator. Also array2D is not a pointer to the data, but rather a pointer to the vector of pointers to each row. This doesn't work in MPI and there are countless of other questions here on StackOverflow, where this fact is discussed - just search for MPI 2D and see for yourself.

The solution is to use a big chunk of singly allocated memory block with an accompanying dope vector - see this question and the arralloc() function mentioned in the answer.

Answer 3

This problem, involving array allocations, comes up all the time in dealing with C/C++ and MPI. This:

int **array2D;
array2D = new int*[dim];
for (int i=0; i<dim; i++) {
    array2D[i] = new int[dim](); // the extra '()' initializes to zero.
}

allocates dim 1d arrays, each dim ints in length. However, there's no reason at all why these should be laid out next to each other - the dim arrays are likely scattered across memory. So even sending dim*dim ints from array2D[0] won't work. The all_array2D is the same; you are creating size*dim arrays, each of size dim, but where they are in relation to each other who knows, making your displacements likely wrong.

To make the arrays contiguous in memory, you need to do something like

int **array2D;
array2D = new int*[dim];
array2D[0] = new int[dim*dim];
for (int i=1; i<dim; i++) {
    array2D[i] = &(array2D[dim*i]);
}

and similarly for all_array2D. Only then can you start reasoning about memory layouts.