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