The listing of ninth.c

#include <stdio.h>
#include <string.h>
#include <sys/types.h>
#include <stdlib.h>
#include <math.h>
#include "mpi.h"

int main ( int argc, char **argv )
{
    int pool_size, my_rank, host_rank, destination, source, node_name_length,
        *host_rank_ptr, found_flag, int_buffer[BUFSIZ];
    char char_buffer[BUFSIZ], my_node_name[BUFSIZ], console_name[BUFSIZ],
         host_name[BUFSIZ]; 
    MPI_Status status;
    boolean_t i_am_the_master = B_FALSE; 

#include "preamble.c"

    {

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

       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;

       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, host_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 ();
}