Centre Blaise Pascal » Bench4GPU

//

// Estimation of Pi using Monte Carlo exploration process

// gcc -std=c99 -O3 -o Pi Pi.c -lm 

// Emmanuel Quemener <emmanuel.quemener@ens-lyon.fr>

// Needed for gethostname

#define _BSD_SOURCE

#include <sys/unistd.h>

#include <math.h>

#include <stdio.h>

#include <stdlib.h>

#include <limits.h>

#include <mpi.h>

#ifdef TIME

#include <sys/time.h>

#endif

// Marsaglia RNG very simple implementation

#define znew  ((z=36969*(z&65535)+(z>>16))<<16)

#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)

#define MWC   (znew+wnew)

#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))

#define CONG  (jcong=69069*jcong+1234567)

#define KISS  ((MWC^CONG)+SHR3)

#define MWCfp MWC * 2.328306435454494e-10f

#define KISSfp KISS * 2.328306435454494e-10f

#define ITERATIONS 1000000000

#ifdef LONG

#define LENGTH long long

#else

#define LENGTH int

#endif

unsigned int rotl(unsigned int value, int shift) {

    return (value << shift) | (value >> (sizeof(value) * CHAR_BIT - shift));

}

unsigned int rotr(unsigned int value, int shift) {

    return (value >> shift) | (value << (sizeof(value) * CHAR_BIT - shift));

}

LENGTH MainLoopGlobal(LENGTH iterations,unsigned int seed_w,unsigned int seed_z)

{

   unsigned int z=seed_z;

   unsigned int w=seed_w;

   LENGTH total=0;

   for (LENGTH i=0;i<iterations;i++) {

      float x=MWCfp ;

      float y=MWCfp ;

      // Matching test

      int inside=((x*x+y*y) < 1.0f) ? 1:0;

      total+=inside;

   }

   return(total);

}

int main(int argc, char *argv[]) {

  unsigned int seed_z=362436069,seed_w=52128862;

  LENGTH iterations=ITERATIONS,inside[1024],insides,part_inside,part_iterations;

  int numtasks,rank,rc,tag=1,i;

  float pi;

  char hostname[128];

  gethostname(hostname, sizeof hostname);

#ifdef TIME

  struct timeval start,end;

  long int useconds;

#endif

  MPI_Status Stat;

  rc = MPI_Init(&argc,&argv);

  if (rc != MPI_SUCCESS) {

    printf ("Error starting MPI program. Terminating.\n");

    MPI_Abort(MPI_COMM_WORLD, rc);

  }

  MPI_Comm_size(MPI_COMM_WORLD,&numtasks);

  MPI_Comm_rank(MPI_COMM_WORLD,&rank);

  if (rank==0) {

    if (argc > 1) {

      iterations=(LENGTH)atoll(argv[1]);

    }

    else {

      printf("\n\tPi : Estimate Pi with Monte Carlo exploration\n\n");

      printf("\t\t#1 : number of iterations (default 1 billion)\n\n");

    }

    printf ("Sizeof int = %lld bytes.\n", (long long)sizeof(int));

    printf ("Sizeof long = %lld bytes.\n", (long long)sizeof(long));

    printf ("Sizeof long long = %lld bytes.\n", (long long)sizeof(long long));

    printf ("Max int = %u\n", INT_MAX);

    printf ("Max long = %ld\n", LONG_MAX);

    printf ("Max long long = %lld\n", LLONG_MAX);

    part_iterations=iterations/numtasks+1;

    // Split part of code

    for (i=1;i<numtasks;i++) {

#ifdef LONG

      rc = MPI_Send(&part_iterations, 1, MPI_LONG_LONG, i, tag, 

                    MPI_COMM_WORLD);

#else

      rc = MPI_Send(&part_iterations, 1, MPI_INT, i, tag, 

                    MPI_COMM_WORLD);

#endif

    }

#ifdef TIME

    gettimeofday(&start,(struct timezone *)0);

#endif

    insides=MainLoopGlobal(part_iterations,seed_w,seed_z);

#ifdef TIME

    gettimeofday(&end,(struct timezone *)0);

    useconds=(end.tv_sec-start.tv_sec)*1000000+end.tv_usec-start.tv_usec;

      printf("\tOn %s with rank %i, find %lld inside in %lu useconds.\n",

             hostname,rank,(long long)insides,useconds);

#else

      printf("\tOn %s with rank %i, find %lld inside\n",hostname,rank,

             (long long)insides);

#endif

    // Join part of code

      for (i=1;i<numtasks;i++) {

#ifdef LONG

        rc = MPI_Recv(&inside[i], 1, MPI_LONG_LONG, i, tag, 

                      MPI_COMM_WORLD, &Stat);

#else

        rc = MPI_Recv(&inside[i], 1, MPI_INT, i, tag, 

                      MPI_COMM_WORLD, &Stat);

#endif

        printf("\tReceive %lu inside from rank %i\n",(unsigned long)inside[i],i);

        insides+=inside[i];

      }

      pi=4.*(float)insides/(float)((iterations/numtasks)*numtasks);

      printf("\n\tPi=%.40f\n\twith error %.40f\n\twith %lld iterations\n\n",pi,

             fabs(pi-4*atan(1.))/pi,(long long)iterations);

  }

  else

    {

      // Receive information from master

#ifdef LONG

      rc = MPI_Recv(&part_iterations, 1, MPI_LONG_LONG, 0, tag, 

                    MPI_COMM_WORLD, &Stat);

#else

      rc = MPI_Recv(&part_iterations, 1, MPI_INT, 0, tag, 

                    MPI_COMM_WORLD, &Stat);

#endif

      printf("\tOn %s with rank %i, receive from master %lld\n",

             hostname,rank,(long long)part_iterations);

#ifdef TIME

      gettimeofday(&start,(struct timezone *)0);

#endif

      part_inside=MainLoopGlobal(part_iterations,rotr(seed_w,rank),rotl(seed_z,rank));

#ifdef TIME

      gettimeofday(&end,(struct timezone *)0);

      useconds=(end.tv_sec-start.tv_sec)*1000000+end.tv_usec-start.tv_usec;

      printf("\tOn %s with rank %i, find %lld inside in %lu useconds.\n",

             hostname,rank,(long long)part_inside,useconds);

#else

      printf("\tOn %s with rank %i, find %lld inside\n",hostname,rank,

             (long long)part_inside);

#endif

#ifdef LONG

      rc = MPI_Send(&part_inside, 1, MPI_LONG_LONG, 0, tag, MPI_COMM_WORLD);

#else

      rc = MPI_Send(&part_inside, 1, MPI_INT, 0, tag, MPI_COMM_WORLD);

#endif

    }

  MPI_Finalize();

}