3

I have a 3D array Foo3D (50 x 100 x 100) in ranks 0 and 1. Foo3D is allocated like this:

int nx = 50;
int ny = 100;
int nz = 100;
typedef int nRarray[100][50];
nRarray *Foo3D;
if ((Foo3D = (nRarray *)malloc((nx*ny*nz)*sizeof(int))) == 0) {fprintf(stderr,"malloc1 Fail \n"); return 1;}

I assign some numbers to Foo3D in rank 0 and save it to a new 2D array (Foo2D) as:

if (myrank == 0) {
for (int j = 0; j < ny; j++) {
 for (int k = 0; k < nz; k++) {

   Foo3D[0][j][k] = j + k;
   Foo2D[j][k] = Foo3D[0][j][k];
}
 } 
}

Now, I'm interested to send Foo2D to rank 1, and place it in its position in Foo3D. In fact, I know I can send Foo2D to rank 1 as:

if (myrank == 0)
{
 MPI_Send(Foo2D,sizeof_Foo2D,MPI_INT,1,100,MPI_COMM_WORLD);
}
else if (myrank == 1)
{
 MPI_Recv(Foo2D,sizeof_Foo2D,MPI_INT,0,100,MPI_COMM_WORLD, &status);
}

And then assign the received Foo2D in rank 1 to its position in Foo3D as:

if (myrank == 1)
{
for (int j = 0; j < ny; j++) {
 for (int k = 0; k < nz; k++) {

   Foo3D[0][j][k] = Foo2D[j][k];
}
 } 
}

Instead of using this procedure and Foo2D as an intermediate variable, is it possible to send the slice of Foo3D in rank 0 directly to its equivalent position in rank 1 or not? In fact, I don't want to send the whole Foo3D to rank 1 because it's a really large array and I'm interested to send just a slice of it to rank 1.

1 Answers1

3

You are defining the array Foo3D as int[nx][ny][nz]. Because C/C++ are row-major languages, the elements of the grid face represented as Foo3D[0][j][k] are actually continuous in memory.

Memory Layout for i=0

So you can simply send the face where i=0 using:

if (myrank == 0) {
    MPI_Send(Foo3D, ny*nz, MPI_INT, 1, 100, MPI_COMM_WORLD);
}
else if (myrank == 1){
    MPI_Recv(Foo3D, ny*nz, MPI_INT, 0, 100, MPI_COMM_WORLD, &status);
}

On the other hand, if you want to send non-contiguous data without using a buffer (Zero-Copy), you can create a custom MPI datatype that represents the data you want to copy, and send directly from the source buffer, and MPI will read the data you specified even if it is non-contiguous in memory. This can be done using MPI_Type_vector.

For example, if you want to send the values where k=0, that is, the grid face Foo[i][j][0].

Memory Layout for k=0

First, create a data type representing the face of the grid that you want to send.

Diagram of MPI_Type_vector inputs

// Create a data type and save its size
MPI_Datatype cubeface;
int cubefacesize;
MPI_Type_vector(nx*ny, 1, nz, MPI_INT, &cubeface);
MPI_Type_commit(&cubeface);
MPI_Type_size(cubeface, &cubefacesize);

Then, you can send and receive using:

if (myrank == 0) {
    MPI_Send(Foo3D, 1, cubeface, 1, cubefacesize, MPI_COMM_WORLD);
} else if (myrank == 1) {
    MPI_Recv(Foo3D, 1, cubeface, 0, cubefacesize, MPI_COMM_WORLD, &status);
}
Aziz
  • 20,065
  • 8
  • 63
  • 69
  • Your derived datatype would be correct in Fortran, but this is C here. – Gilles Gouaillardet Aug 21 '19 at 12:20
  • @GillesGouaillardet You mean row and major column arrays? Cause in FORTRAN the arrays are column major but in C the arrays are row major, right? – Mithridates the Great Aug 21 '19 at 13:25
  • @GillesGouaillardet You are right! Thank you. I'll update my answer. – Aziz Aug 21 '19 at 15:50
  • @AloneProgrammer Yes, I believe this is what Gilles meant. I have updated the answer accordingly. – Aziz Aug 21 '19 at 16:52
  • @Aziz I think still something is not correct here... Do you flatten your 3D array by using this formula: i+j*n_{x}+k*n_{x}*n_{y}, if yes, so the XY plane is continuous not YZ plane... – Mithridates the Great Aug 21 '19 at 17:06
  • @AloneProgrammer If you're defining the array as `int[nx][ny][nz]`, then the order of elements in memory is `Foo3D[0][0][0]`, `Foo3D[0][0][1]`, `Foo3D[0][0][2]`, .. etc. Hence, the formula is `k + j*Nz + i*Nz*Ny`. This is because C/C++ are row-major. – Aziz Aug 21 '19 at 17:39
  • @Aziz I don't think so. The correct formula must be: `i+j*Nx+k*Nx*Ny`. See the discussion here: https://stackoverflow.com/a/16610069/10382271 – Mithridates the Great Aug 21 '19 at 18:04
  • @AloneProgrammer I understand the source of confusion. Your formula is correct **if** you define the array as `int[nz][ny][nx]` and access the elements as `Foo3D[k][j][i]`. But since you are using `int[nx][ny][nz]`, the mapping is done in the order `z,y,x` and not `x,y,z`. This example shows how elements are stored in memory https://godbolt.org/z/-oblH4. This is a C/C++ thing (because they're row-major), and because you are using the syntax `int[][][]`. – Aziz Aug 21 '19 at 18:18
  • 1
    Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/198260/discussion-between-aziz-and-alone-programmer). – Aziz Aug 21 '19 at 18:23