Particles Again

In this section we are going to have a yet another look at the problem of handling particles in MPI. But our particles are now going to be described as God meant them to be, i.e., they are going to have some character descriptors, some class descriptors (this is going to be an integer), and some floating point descriptors, e.g., a vector of coordinates and velocities. And the we are going to do various things with them in order to illustrate the use of the functions discussed in this section.

The following is not a real, complete code. Instead these are various excerpts from a code that you may wish to write one day, so the ideas presented here will come handy.

So let us assume that we have done our #include and main and now we get down to the interesting bits. First some declarations. Let a structure that represents particles be defined as follows:

struct Partstruct 
   { 
   int    class;  /* particle class */ 
   double d[6];   /* particle coordinates */ 
   char   b[7];   /* some additional information */ 
   };

Now we define an array that is going to house 1000 of these particles:

struct Partstruct    particle[1000];

In order to represent the particles to MPI we have to prepare the ground for calling MPI_Struct. So we begin from the following definitions:

MPI_Datatype Particletype; 
MPI_Datatype type[3] = {MPI_INT, MPI_DOUBLE, MPI_CHAR}; 
int          blocklen[3] = {1, 6, 7}; 
MPI_Aint     disp[3]; 
int          base;

Given a particle structure, as implemented by the compiler, and we really don't know a priori how that is going to be done, and we want to write a portable program, we obtain the information about data layout within the structure Partstruct by interrogating this structure with function MPI_Address:

 
MPI_Address( particle, disp); 
MPI_Address( particle[0].d, disp+1); 
MPI_Address( particle[0].b, disp+2); 
base = disp[0]; 
for (i=0; i <3; i++) disp[i] -= base;

We have subtracted disp[0] from all other disps so as to obtain displacements from the beginning of a buffer, as opposed to absolute addresses. This way we don't have to assume while writing the program that data is going to be laid out in some specific way.

Now we are ready to construct an MPI structure that matches the C-defined structure:

MPI_Type_struct( 3, blocklen, disp, type, &Particletype);

Looking at the definitions of the variables used in this call we are constructing a polymorphic data type which is made of 3 blocks of data. The first block contains 1 variable of type MPI_INT, and its relative displacement is 0. The second block contains 6 variables of type MPI_DOUBLE and its displacement is given by disp[1], which we have calculated by looking at the absolute address of this block for a real object of type Partstruct. The third block of data comprises 7 characters, and its displacement is given by disp[2].

Now we have to commit this newly created MPI data type and from this point onwards we can use it in order to send and receive particles:

MPI_Type_commit( &Particletype); 
MPI_Send( particle, 1000, Particletype, dest, tag, comm);

But we can do better. Say that some of the particles happen to be of class zero and that we want to send only particles in this class. How to do this? Let us begin by adding some more definitions for a new data type:

MPI_Datatype Zparticles;
MPI_Datatype Ztype; 
MPI_Aint     zdisp[1000]; 
int zblock[1000], j, k; 
int zzblock[2] = {1,1}; 
MPI_Aint     zzdisp[2]; 
MPI_Datatype zztype[2];

When we get through to the part of the code that is going to construct the new data type that corresponds to particles of class zero we begin by computing displacements for these particles:

j = 0; 
for(i=0; i < 1000; i++) 
  if (particle[i].class==0) 
     { 
     zdisp[j] = i; 
     zblock[j] = 1; 
     j++; 
     }

Since we have already defined an MPI type for a particle we can use a simpler function, namely MPI_Type_indexed in order to define the type for particles of class zero:

MPI_Type_indexed( j, zblock, zdisp, Particletype, &Zparticles);

This data type comprises j blocks, each of which contains one particle. The blocks are picked up from locations given by the array zdisp. The newly constructured data type is written on Zparticles.

We could now send the lot across to another process in just one operation. But it would be nice to also add information about the number of particles we are going to send across. This number sits inside j, and we are going to prepend it to Zparticles and create another compound MPI data type, which comprises one integer, namely the number of particles that we are going to send, and then all those particles.

First we find what are absolute addresses of j and of the first particle, then we write those addresses on zzdisp, which are the displacements for the new data type. The array that specifies the corresponding types

MPI_Address(&j, zzdisp); 
MPI_Address(particle, zzdisp+1);

Before we can call MPI_Type_struct we also have to construct an array of types:

zztype[0] = MPI_INT; 
zztype[1] = Zparticles;

Finally we construct this new type that comprises an integer followed by all particles of class zero, we commit this type (observe that we didn't have to commit the intermediate type Zparticles) and send the lot across. Because we have used absolute addresses, the pointer to the buffer is MPI_BOTTOM.

MPI_Type_struct(2, zzblock, zzdisp, zztype, &Ztype); 
MPI_Type_commit( &Ztype); 
MPI_Send( MPI_BOTTOM, 1, Ztype, dest, tag, comm);

Now let us say that we want to send the first two coordinates for all particles, i.e., particle[j].d[0] and particle[j].d[1], for all j. Here are the definitions that we have to introduce to accomplish this task:

MPI_Datatype Allpairs; 
MPI_Aint sizeofentry;

Since we have already constructed the MPI type for a single particle, that type is Particletype, we can use this in order to find stride in bytes:

MPI_Type_extent( Particletype, &sizeofentry);

The type Allpairs can be now defined by calling MPI_Type_hvector:

MPI_Type_hvector( 1000, 2, sizeofentry, MPI_DOUBLE, &Allpairs);

Here we are taking 1000 blocks of 2 double precision floating point numbers each from an otherwise unspecified array of data (it may be polymorphic) and the blocks are separated by a stride of sizeofentry. We commit this newly constructed MPI data type, then locate the first coordinate of the first particle and send the whole lot in one go:

MPI_Type_commit( &Allpairs); 
MPI_Send( particle[0].d, 1, Allpairs, dest, tag, comm);