The code is straight forward and self-explainitory.
#include <stdio.h> #include <mpi.h> #include <math.h> #include <stddef.h> #include <stdlib.h> #define d_sphere 5 /** the number of dimensions of the sphere ***/ #define num_trials 1000000 /** the number of trials for each processor ***/ void srand48(); double drand48(); void main(int argc, char *argv[]) { int numprocs; /* Number of processors */ int procnum; /* Processor number */ int j,i,*int_arr; double volume,pi,r2,avg,sigma,volume_arr[4],sigma_arr[4]; MPI_Status status; MPI_Init(&argc, &argv); MPI_Comm_rank(MPI_COMM_WORLD, &procnum); MPI_Comm_size(MPI_COMM_WORLD, &numprocs); srand48 ( (long) procnum); int_arr = (int *) malloc((num_trials + 1) * sizeof(int)); for(i=1;i<=num_trials;i++) { /*** calculate one point ***/ r2 = 0.0; for(j=1;j<=d_sphere;j++) r2 += pow(drand48(),2.0); if(r2 <= 1.0) { int_arr[i] = 1; } else { int_arr[i] = 0; } } /*** now get average and standard deviation ***/ avg = 0.0; for(i=1;i<=num_trials;i++) avg += ((double) int_arr[i]); avg /= ((double) num_trials); volume = avg * pow(2.0,d_sphere); sigma = 0.0; for(i=1;i<=num_trials;i++) sigma += pow(((double) int_arr[i]) - avg,2.0); sigma = sigma / ((double) num_trials); sigma = sqrt(sigma); /*** send and receive global data ***/ if(procnum == 0) { volume_arr[0] = volume; sigma_arr[0] = sigma; for(i=1;i<numprocs;i++) { MPI_Recv(&(volume_arr[i]),1,MPI_DOUBLE,i,0,MPI_COMM_WORLD,&status); MPI_Recv(&(sigma_arr[i]),1,MPI_DOUBLE,i,0,MPI_COMM_WORLD,&status); printf("%d vol is %25.16e \n",i,volume_arr[i]); } } else { MPI_Send(&volume,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD); MPI_Send(&sigma,1,MPI_DOUBLE,0,0,MPI_COMM_WORLD); } /*** compute and print out global information ***/ if(procnum == 0) { /*** compute global averages ***/ sigma = 0.0; volume = 0.0; for(i=0;i<numprocs;i++) { volume += volume_arr[i]; sigma += sigma_arr[i]; } volume /= ((double) numprocs); sigma /= ((double) numprocs); printf("--------------------------------------------- \n"); printf("VALUES AVERAGED OVER ALL PROCESSORS: \n"); printf(" volume of %d sphere is %25.16e \n",d_sphere,volume); if(d_sphere == 2) printf("in 2 dim, volume = pi = %25.16e \n",volume); printf(" st. dev. is %25.16e \n",sigma); printf(" total number of points: %d \n",4*num_trials); printf("--------------------------------------------- \n"); } MPI_Finalize(); }