next up previous index
Next: The Discussion Up: Interacting Particles Previous: Interacting Particles

The Code

#include <stdio.h>
#include <math.h>
#include <mpi.h>

#define FALSE            0
#define TRUE             1
#define MASTER_RANK      0

#define MAX_PARTICLES 1000
#define MAX_PROCS      128
#define EPSILON          1.0E-10


int main ( int argc, char **argv )
{
   int pool_size, my_rank;
   int i_am_the_master = FALSE; 
   extern double drand48();
   extern void srand48();

   typedef struct {
      double x, y, z;
      double mass;
   } Particle;

   typedef struct {
      double vx, vy, vz;
   } ParticleV;

   Particle  particles[MAX_PARTICLES];  /* Particles on all nodes */
   ParticleV vector[MAX_PARTICLES]; /* Particle velocity */
   int       counts[MAX_PROCS];         /* Number of ptcls on each proc */
   int       displacements[MAX_PROCS];  /* Offsets into particles */
   int       offsets[MAX_PROCS];        /* Offsets used by the master */
   int       particle_number, i, j, 
             my_offset;                 /* Location of local particles */
   int       total_particles;           /* Total number of particles */
   int       count;                     /* Number of times in loop */

   MPI_Datatype particle_type;

   double    time;                      /* Computation time */
   double    dt, dt_old;                /* Integration time step */
   double    a0, a1, a2;
   double    start_time, end_time;

   MPI_Init(&argc, &argv);
   MPI_Comm_size(MPI_COMM_WORLD, &pool_size);
   MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

   if (my_rank == MASTER_RANK) i_am_the_master = TRUE;

   particle_number = MAX_PARTICLES / pool_size;

   if (i_am_the_master)
      printf ("%d particles per processor\n", particle_number);

   MPI_Type_contiguous ( 4, MPI_DOUBLE, &particle_type );
   MPI_Type_commit ( &particle_type );

   MPI_Allgather ( &particle_number, 1, MPI_INT, counts, 1, MPI_INT,
                   MPI_COMM_WORLD );

   displacements[0] = 0;
   for (i = 1; i < pool_size; i++)
      displacements[i] = displacements[i-1] + counts[i-1];
   total_particles = displacements[pool_size - 1]
                     + counts[pool_size - 1];
       
   if (i_am_the_master)
      printf ("total number of particles = %d\n", total_particles);

   my_offset = displacements[my_rank];

   MPI_Gather ( &my_offset, 1, MPI_INT, offsets, 1, MPI_INT, MASTER_RANK,
                MPI_COMM_WORLD );

   if (i_am_the_master) {
      printf ("offsets: ");
      for (i = 0; i < pool_size; i++)
         printf ("%d ", offsets[i]);
      printf("\n");
   }

   srand48((long) my_rank);

   for (i = 0; i < particle_number; i++) {
      particles[my_offset + i].x = drand48();
      particles[my_offset + i].y = drand48();
      particles[my_offset + i].z = drand48();
      particles[my_offset + i].mass = 1.0;
   }

   start_time = MPI_Wtime();

   MPI_Allgatherv ( particles + my_offset, particle_number,
                    particle_type,
                    particles, counts, displacements, particle_type,
                    MPI_COMM_WORLD );

   end_time = MPI_Wtime();

   if (i_am_the_master) {
      printf ("Communicating = %8.5f seconds\n", end_time - start_time);
      printf ("particles[offsets[i]].x: ");
      for (i = 0; i < pool_size; i++)
         printf ("%8.5f ", particles[offsets[i]].x);
      printf("\n"); fflush(stdout);
   }

   start_time = MPI_Wtime();

   count = 0;
   for (i = 0; i < particle_number; i++) {

      vector[my_offset + i].vx = 0.0;
      vector[my_offset + i].vy = 0.0;
      vector[my_offset + i].vz = 0.0;
          
      for (j = 0; j < total_particles; j++) {
         if (j != i) {

            double dx, dy, dz, r2, r, mimj_by_r3;

            dx = particles[my_offset + i].x - particles[j].x;
            dy = particles[my_offset + i].y - particles[j].y;
            dz = particles[my_offset + i].z - particles[j].z;

            r2 = dx * dx + dy * dy + dz * dz; r = sqrt(r2);

            if (r2 < EPSILON) mimj_by_r3 = 0.0;
            else
               mimj_by_r3 = particles[my_offset + i].mass 
                            * particles[j].mass / (r2 * r);

            vector[my_offset + i].vx = vector[my_offset + i].vx +
                                       mimj_by_r3 * dx;
            vector[my_offset + i].vy = vector[my_offset + i].vy +
                                       mimj_by_r3 * dy;
            vector[my_offset + i].vz = vector[my_offset + i].vz +
                                       mimj_by_r3 * dz;
            count = count + 1;
         }
      }
   }

   end_time = MPI_Wtime();
   if (i_am_the_master)
      printf ("done my job in %8.5f seconds, waiting for slow \
processes...\n", end_time - start_time);

   MPI_Barrier (MPI_COMM_WORLD);

   end_time = MPI_Wtime();

   if (i_am_the_master)
          printf ("evaluated %d 3D interactions in %8.5f seconds\n",
                  count * pool_size, end_time - start_time);

   MPI_Finalize ();
}

And here is how this code compiles and runs:

gustav@sp20:../MPI 18:43:30 !588 $ mpcc -o particles particles.c -lm
gustav@sp20:../MPI 18:43:47 !589 $ cat particles.ll
# @ job_type = parallel
# @ environment = COPY_ALL; MP_EUILIB=ip; MP_INFOLEVEL=2
# @ requirements = (Adapter == "hps_ip") && (Machine != "sp43") && (Machine != "sp42")
# @ min_processors = 6
# @ max_processors = 10
# @ class = test
# @ output = particles.out
# @ error = particles.err
# @ executable = /usr/bin/poe
# @ arguments = particles
# @ queue
gustav@sp20:../MPI 18:43:52 !590 $ llsubmit particles.ll
submit: The job "sp20.123" has been submitted.
gustav@sp20:../MPI 18:43:59 !591 $ llq | grep gustav
gustav@sp20:../MPI 18:44:39 !592 $ cat particles.out
100 particles per processor
total number of particles = 1000
offsets: 0 100 200 300 400 500 600 700 800 900 
Communicating =  0.00802 seconds
particles[offsets[i]].x:  0.17083  0.04163  0.91243  0.78323  0.65404  0.52484  0.39564  0.26644  0.13725  0.00805 
done my job in  0.09365 seconds, waiting for slow processes...
evaluated 999000 3D interactions in  0.09555 seconds
gustav@sp20:../MPI 18:44:43 !593 $

The messages that the program has written on standard output may not make that much sense to you right now, but they'll become clearer as we discuss the code in detail.



Zdzislaw Meglicki
2001-02-26