examples/smpi/NAS/ep.c

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <math.h>
   5 #include "smpi/mpi.h"
   6 #include "nas_common.h"
   7 #include "simgrid/instr.h" //TRACE_
   8
   9 char class;
  10 int nprocs;
  11
  12 #define TRUE 1
  13 #define FALSE 0
  14
  15 int main(int argc, char **argv) {
  16   double dum[3] = {1.,1.,1.};
  17   double x1, x2, sx, sy, tm, an, gc;
  18   double Mops;
  19   double epsilon=1.0E-8, a = 1220703125., s=271828183.;
  20   double t1, t2, t3, t4;
  21   double sx_verify_value, sy_verify_value, sx_err, sy_err;
  22
  23   int    m;
  24   int    mk=16;
  25   int    nk = (int)(pow(2,mk)),
  26          nq=10,
  27          np, node, no_nodes, i, ik, kk, l, k, nit, no_large_nodes, np_add, k_offset;
  28   int    verified;
  29   char   size[500]; // mind the size of the string to represent a big number
  30
  31   double *x = (double *) SMPI_SHARED_MALLOC (2*nk*sizeof(double));
  32   double *q = (double *) SMPI_SHARED_MALLOC (nq*sizeof(double));
  33
  34   MPI_Init( &argc, &argv );
  35   MPI_Comm_size( MPI_COMM_WORLD, &no_nodes);
  36   MPI_Comm_rank( MPI_COMM_WORLD, &node);
  37
  38   TRACE_smpi_set_category ("start");
  39
  40   get_info(argc, argv, &nprocs, &class);
  41   check_info(EP, nprocs, class);
  42
  43   if      (class == 'S')
  44     { m = 24; }
  45   else if (class == 'W')
  46     { m = 25; }
  47   else if (class == 'A')
  48     { m = 28; }
  49   else if (class == 'B')
  50     { m = 30; }
  51   else if (class == 'C')
  52     { m = 32; }
  53   else if (class == 'D')
  54     { m = 36; }
  55   else if (class == 'E')
  56     { m = 40; }
  57   else {
  58     printf("EP: Internal error: invalid class type %c\n", class);
  59     exit(1);
  60   }
  61
  62   int mm = m -mk;
  63   int nn = (int)(pow(2,mm));
  64
  65   int root = 0;
  66   if (node == root ) {
  67     /* Because the size of the problem is too large to store in a 32-bit integer for some classes, we put it into a
  68      * string (for printing). Have to strip off the decimal point put in there by the floating point print statement
  69      * (internal file) */
  70     fprintf(stdout," NAS Parallel Benchmarks 3.2 -- EP Benchmark");
  71     snprintf(size,500,"%lu",(unsigned long)pow(2,m+1));
  72     fprintf(stdout," Number of random numbers generated: %s\n",size);
  73     fprintf(stdout," Number of active processes: %d\n",no_nodes);
  74   }
  75
  76   /* Compute the number of "batches" of random number pairs generated per processor. Adjust if the number of processors
  77    * does not evenly divide the total number  */
  78   np = nn / no_nodes;
  79   no_large_nodes = nn % no_nodes;
  80   if (node < no_large_nodes)
  81     np_add = 1;
  82   else
  83     np_add = 0;
  84   np = np + np_add;
  85
  86   if (np == 0) {
  87     fprintf(stdout,"Too many nodes: %d  %d",no_nodes,nn);
  88     MPI_Abort(MPI_COMM_WORLD,1);
  89     exit(0);
  90   }
  91
  92   /* Call the random number generator functions and initialize the x-array to reduce the effects of paging the timings.
  93    Also, call all mathematical functions that are used. Make sure initializations cannot be eliminated as dead code. */
  94   vranlc(0, dum[0], dum[1], &(dum[2]));
  95
  96   dum[0] = randlc(&(dum[1]),&(dum[2]));
  97   for (i=0;i<2*nk;i++) {
  98     x[i] = -1e99;
  99   }
 100
 101   /* Synchronize before placing time stamp */
 102   MPI_Barrier( MPI_COMM_WORLD );
 103
 104   TRACE_smpi_set_category ("ep");
 105
 106   timer_clear(1);
 107   timer_clear(2);
 108   timer_clear(3);
 109   timer_start(1);
 110
 111   t1 = a;
 112   t1 = vranlc(0, t1, a, x);
 113
 114   /* Compute AN = A ^ (2 * NK) (mod 2^46). */
 115   t1 = a;
 116   for (i=1; i <= mk+1; i++) {
 117     t2 = randlc(&t1, &t1);
 118   }
 119   an = t1;
 120   gc = 0;
 121   sx = 0.;
 122   sy = 0.;
 123   for (i=0; i < nq ; i++) {
 124     q[i] = 0.;
 125   }
 126
 127 /* Each instance of this loop may be performed independently. We compute the k offsets separately to take into account
 128  * the fact that some nodes have more numbers to generate than others */
 129   if (np_add == 1)
 130     k_offset = node * np -1;
 131   else
 132     k_offset = no_large_nodes*(np+1) + (node-no_large_nodes)*np -1;
 133
 134   int stop = FALSE;
 135   SMPI_SAMPLE_GLOBAL(k = 1, k <= np, k++, 0.25 * np, 0.03){
 136     stop = FALSE;
 137     kk = k_offset + k ;
 138     t1 = s;
 139     t2 = an;
 140
 141     //       Find starting seed t1 for this kk.
 142     for (i=1;i<=100 && !stop;i++) {
 143       ik = kk / 2;
 144       if (2 * ik != kk)  {
 145         randlc(&t1, &t2);
 146       }
 147       if (ik==0)
 148         stop = TRUE;
 149       else {
 150         randlc(&t2, &t2);
 151         kk = ik;
 152       }
 153     }
 154     //       Compute uniform pseudorandom numbers.
 155
 156     timer_start(3);
 157
 158     t1 = vranlc(2 * nk, t1, a, x);
 159
 160     timer_stop(3);
 161
 162     /* Compute Gaussian deviates by acceptance-rejection method and tally counts in concentric square annuli.
 163      * This loop is not vectorizable. */
 164     timer_start(2);
 165
 166     for(i=1; i<=nk;i++) {
 167       x1 = 2. * x[2*i-2] -1.0;
 168       x2 = 2. * x[2*i-1] - 1.0;
 169       t1 = x1*x1 + x2*x2;
 170       if (t1 <= 1.) {
 171         t2   = sqrt(-2. * log(t1) / t1);
 172         t3   = (x1 * t2);
 173         t4   = (x2 * t2);
 174         l    = (int)(fabs(t3) > fabs(t4) ? fabs(t3) : fabs(t4));
 175         q[l] = q[l] + 1.;
 176         sx   = sx + t3;
 177         sy   = sy + t4;
 178       }
 179     }
 180     timer_stop(2);
 181   }
 182
 183   TRACE_smpi_set_category ("finalize");
 184
 185   MPI_Allreduce(&sx, x, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 186   sx = x[0];
 187   MPI_Allreduce(&sy, x, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 188   sy = x[0];
 189   MPI_Allreduce(q, x, nq, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 190
 191   for(i = 0; i < nq; i++) {
 192     q[i] = x[i];
 193   }
 194   for(i = 0; i < nq; i++) {
 195     gc += q[i];
 196   }
 197
 198   timer_stop(1);
 199   tm = timer_read(1);
 200   MPI_Allreduce(&tm, x, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
 201   tm = x[0];
 202
 203   SMPI_SHARED_FREE(x);
 204
 205   if(node == root) {
 206     nit = 0;
 207     verified = TRUE;
 208
 209     if(m == 24) {
 210       sx_verify_value = -3.247834652034740E3;
 211       sy_verify_value = -6.958407078382297E3;
 212     } else if(m == 25) {
 213       sx_verify_value = -2.863319731645753E3;
 214       sy_verify_value = -6.320053679109499E3;
 215     } else if(m == 28) {
 216       sx_verify_value = -4.295875165629892E3;
 217       sy_verify_value = -1.580732573678431E4;
 218     } else if(m == 30) {
 219       sx_verify_value =  4.033815542441498E4;
 220       sy_verify_value = -2.660669192809235E4;
 221     } else if(m == 32) {
 222       sx_verify_value =  4.764367927995374E4;
 223       sy_verify_value = -8.084072988043731E4;
 224     } else if(m == 36) {
 225       sx_verify_value =  1.982481200946593E5;
 226       sy_verify_value = -1.020596636361769E5;
 227     } else {
 228       verified = FALSE;
 229     }
 230
 231     if(verified) {
 232       sx_err = fabs((sx - sx_verify_value)/sx_verify_value);
 233       sy_err = fabs((sy - sy_verify_value)/sy_verify_value);
 234       verified = ((sx_err < epsilon) && (sy_err < epsilon));
 235     }
 236
 237     Mops = (pow(2.0, m+1))/tm/1000;
 238
 239     fprintf(stdout,"EP Benchmark Results:\n");
 240     fprintf(stdout,"CPU Time=%d\n",(int) tm);
 241     fprintf(stdout,"N = 2^%d\n",m);
 242     fprintf(stdout,"No. Gaussain Pairs =%d\n",(int) gc);
 243     fprintf(stdout,"Sum = %f %ld\n",sx,(long) sy);
 244     fprintf(stdout,"Count:");
 245     for(i = 0; i < nq; i++) {
 246       fprintf(stdout,"%d\t %ld\n",i,(long) q[i]);
 247     }
 248     c_print_results("EP", class, m+1, 0, 0, nit, nprocs, no_nodes, tm, Mops, "Random number generated",verified);
 249
 250     fprintf(stdout,"Total time:     %f\n",(timer_read(1)/1000));
 251     fprintf(stdout,"Gaussian pairs: %f\n",(timer_read(2)/1000));
 252     fprintf(stdout,"Random numbers: %f\n",(timer_read(3)/1000));
 253   }
 254
 255   SMPI_SHARED_FREE(q);
 256
 257   MPI_Finalize();
 258   return 0;
 259 }