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