examples/smpi/NAS/EP/ep.c

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <string.h>
   4 #include <math.h>
   5
   6 #include "mpi.h"
   7 #include "npbparams.h"
   8
   9 #include "simgrid/instr.h" //TRACE_
  10
  11 #include "randlc.h"
  12
  13 #ifndef CLASS
  14 #define CLASS 'S'
  15 #define NUM_PROCS            1
  16 #endif
  17 #define true 1
  18 #define false 0
  19
  20 //---NOTE : all the timers function have been modified to
  21 //          avoid global timers (privatize these).
  22       // ----------------------- timers ---------------------
  23       void timer_clear(double *onetimer) {
  24             //elapsed[n] = 0.0;
  25             *onetimer = 0.0;
  26       }
  27
  28       void timer_start(double *onetimer) {
  29             *onetimer = MPI_Wtime();
  30       }
  31
  32       void timer_stop(int n,double *elapsed,double *start) {
  33             double t, now;
  34
  35             now = MPI_Wtime();
  36             t = now - start[n];
  37             elapsed[n] += t;
  38       }
  39
  40       double timer_read(int n, double *elapsed) {  /* ok, useless, but jsut to keep function call */
  41             return(elapsed[n]);
  42       }
  43       /********************************************************************
  44        *****************            V R A N L C          ******************
  45        *****************                                 *****************/
  46       double vranlc(int n, double x, double a, double *y)
  47       {
  48         int i;
  49         long  i246m1=0x00003FFFFFFFFFFF;
  50     long  LLx, Lx, La;
  51         double d2m46;
  52
  53 // This doesn't work, because the compiler does the calculation in 32
  54 // bits and overflows. No standard way (without f90 stuff) to specify
  55 // that the rhs should be done in 64 bit arithmetic.
  56 //     parameter(i246m1=2**46-1)
  57
  58       d2m46=pow(0.5,46);
  59
  60 // c Note that the v6 compiler on an R8000 does something stupid with
  61 // c the above. Using the following instead (or various other things)
  62 // c makes the calculation run almost 10 times as fast.
  63 //
  64 // c     save d2m46
  65 // c      data d2m46/0.0d0/
  66 // c      if (d2m46 .eq. 0.0d0) then
  67 // c         d2m46 = 0.5d0**46
  68 // c      endif
  69
  70       Lx = (long)x;
  71       La = (long)a;
  72       //fprintf(stdout,("================== Vranlc ================");
  73       //fprintf(stdout,("Before Loop: Lx = " + Lx + ", La = " + La);
  74   LLx = Lx;
  75   for (i=0; i< n; i++) {
  76       Lx   = Lx*La & i246m1 ;
  77       LLx = Lx;
  78       y[i] = d2m46 * (double)LLx;
  79       /*
  80          if(i == 0) {
  81          fprintf(stdout,("After loop 0:");
  82          fprintf(stdout,("Lx = " + Lx + ", La = " + La);
  83          fprintf(stdout,("d2m46 = " + d2m46);
  84          fprintf(stdout,("LLX(Lx) = " + LLX.doubleValue());
  85          fprintf(stdout,("Y[0]" + y[0]);
  86          }
  87        */
  88   }
  89
  90       x = (double)LLx;
  91       /*
  92       fprintf(stdout,("Change: Lx = " + Lx);
  93       fprintf(stdout,("=============End   Vranlc ================");
  94       */
  95       return x;
  96     }
  97
  98
  99
 100 //-------------- the core (unique function) -----------
 101       void doTest(int argc, char **argv) {
 102       double dum[3] = {1.,1.,1.};
 103       double x1, x2, sx, sy, tm, an, tt, gc;
 104       double Mops;
 105       double epsilon=1.0E-8, a = 1220703125., s=271828183.;
 106       double t1, t2, t3, t4;
 107       double sx_verify_value, sy_verify_value, sx_err, sy_err;
 108
 109 #include "npbparams.h"
 110       int    mk=16,
 111          // --> set by make : in npbparams.h
 112          //m=28, // for CLASS=A
 113          //m=30, // for CLASS=B
 114          //npm=2, // NPROCS
 115          mm = m-mk,
 116          nn = (int)(pow(2,mm)),
 117          nk = (int)(pow(2,mk)),
 118          nq=10,
 119          np,
 120          node,
 121          no_nodes,
 122          i,
 123          ik,
 124          kk,
 125          l,
 126          k, nit, no_large_nodes,
 127          np_add, k_offset, j;
 128       int    me, nprocs, root=0, dp_type;
 129       int verified,
 130           timers_enabled=true;
 131       char  size[500]; // mind the size of the string to represent a big number
 132
 133       //Use in randlc..
 134       int KS = 0;
 135       double R23, R46, T23, T46;
 136
 137       double *qq = (double *) malloc (10000*sizeof(double));
 138       double *start = (double *) malloc (64*sizeof(double));
 139       double *elapsed = (double *) malloc (64*sizeof(double));
 140
 141       double *x = (double *) malloc (2*nk*sizeof(double));
 142       double *q = (double *) malloc (nq*sizeof(double));
 143
 144       TRACE_smpi_set_category ("start");
 145
 146       MPI_Init( &argc, &argv );
 147       MPI_Comm_size( MPI_COMM_WORLD, &no_nodes);
 148       MPI_Comm_rank( MPI_COMM_WORLD, &node);
 149
 150 #ifdef USE_MPE
 151     MPE_Init_log();
 152 #endif
 153       root = 0;
 154       if (node == root ) {
 155
 156           /*   Because the size of the problem is too large to store in a 32-bit
 157            *   integer for some classes, we put it into a string (for printing).
 158            *   Have to strip off the decimal point put in there by the floating
 159            *   point print statement (internal file)
 160            */
 161           fprintf(stdout," NAS Parallel Benchmarks 3.2 -- EP Benchmark");
 162           sprintf(size,"%d",(int)pow(2,m+1));
 163           //size = size.replace('.', ' ');
 164           fprintf(stdout," Number of random numbers generated: %s\n",size);
 165           fprintf(stdout," Number of active processes: %d\n",no_nodes);
 166
 167       }
 168       verified = false;
 169
 170       /* c   Compute the number of "batches" of random number pairs generated
 171          c   per processor. Adjust if the number of processors does not evenly
 172          c   divide the total number
 173 */
 174
 175        np = nn / no_nodes;
 176        no_large_nodes = nn % no_nodes;
 177        if (node < no_large_nodes) np_add = 1;
 178        else np_add = 0;
 179        np = np + np_add;
 180
 181        if (np == 0) {
 182              fprintf(stdout,"Too many nodes: %d  %d",no_nodes,nn);
 183              MPI_Abort(MPI_COMM_WORLD,1);
 184              exit(0);
 185        }
 186
 187 /* c   Call the random number generator functions and initialize
 188    c   the x-array to reduce the effects of paging on the timings.
 189    c   Also, call all mathematical functions that are used. Make
 190    c   sure these initializations cannot be eliminated as dead code.
 191 */
 192
 193    //call vranlc(0, dum[1], dum[2], dum[3]);
 194    // Array indexes start at 1 in Fortran, 0 in Java
 195    vranlc(0, dum[0], dum[1], &(dum[2]));
 196
 197    dum[0] = randlc(&(dum[1]),&(dum[2]));
 198    /////////////////////////////////
 199    for (i=0;i<2*nk;i++) {
 200        x[i] = -1e99;
 201    }
 202    Mops = log(sqrt(abs(1)));
 203
 204    /*
 205       c---------------------------------------------------------------------
 206       c    Synchronize before placing time stamp
 207       c---------------------------------------------------------------------
 208     */
 209         MPI_Barrier( MPI_COMM_WORLD );
 210
 211         TRACE_smpi_set_category ("ep");
 212
 213         timer_clear(&(elapsed[1]));
 214         timer_clear(&(elapsed[2]));
 215         timer_clear(&(elapsed[3]));
 216         timer_start(&(start[1]));
 217
 218         t1 = a;
 219   //fprintf(stdout,("(ep.f:160) t1 = " + t1);
 220         t1 = vranlc(0, t1, a, x);
 221   //fprintf(stdout,("(ep.f:161) t1 = " + t1);
 222
 223
 224 /* c   Compute AN = A ^ (2 * NK) (mod 2^46). */
 225
 226         t1 = a;
 227   //fprintf(stdout,("(ep.f:165) t1 = " + t1);
 228         for (i=1; i <= mk+1; i++) {
 229                t2 = randlc(&t1, &t1);
 230          //fprintf(stdout,("(ep.f:168)[loop i=" + i +"] t1 = " + t1);
 231         }
 232         an = t1;
 233   //fprintf(stdout,("(ep.f:172) s = " + s);
 234         tt = s;
 235         gc = 0.;
 236         sx = 0.;
 237         sy = 0.;
 238         for (i=0; i < nq ; i++) {
 239                q[i] = 0.;
 240         }
 241
 242 /*
 243     Each instance of this loop may be performed independently. We compute
 244     the k offsets separately to take into account the fact that some nodes
 245     have more numbers to generate than others
 246 */
 247
 248       if (np_add == 1)
 249          k_offset = node * np -1;
 250       else
 251          k_offset = no_large_nodes*(np+1) + (node-no_large_nodes)*np -1;
 252
 253       int stop = false;
 254       for(k = 1; k <= np; k++) {
 255          stop = false;
 256          kk = k_offset + k ;
 257          t1 = s;
 258          //fprintf(stdout,("(ep.f:193) t1 = " + t1);
 259          t2 = an;
 260
 261 //       Find starting seed t1 for this kk.
 262
 263          for (i=1;i<=100 && !stop;i++) {
 264             ik = kk / 2;
 265       //fprintf(stdout,("(ep.f:199) ik = " +ik+", kk = " + kk);
 266             if (2 * ik != kk)  {
 267                 t3 = randlc(&t1, &t2);
 268                 //fprintf(stdout,("(ep.f:200) t1= " +t1 );
 269             }
 270             if (ik==0)
 271                 stop = true;
 272             else {
 273                t3 = randlc(&t2, &t2);
 274                kk = ik;
 275            }
 276          }
 277 //       Compute uniform pseudorandom numbers.
 278
 279          //if (timers_enabled)  timer_start(3);
 280    timer_start(&(start[3]));
 281          //call vranlc(2 * nk, t1, a, x)  --> t1 and y are modified
 282
 283   //fprintf(stdout,">>>>>>>>>>>Before vranlc(l.210)<<<<<<<<<<<<<");
 284   //fprintf(stdout,"2*nk = " + (2*nk));
 285   //fprintf(stdout,"t1 = " + t1);
 286   //fprintf(stdout,"a  = " + a);
 287   //fprintf(stdout,"x[0] = " + x[0]);
 288   //fprintf(stdout,">>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<");
 289
 290   t1 = vranlc(2 * nk, t1, a, x);
 291
 292   //fprintf(stdout,(">>>>>>>>>>>After  Enter vranlc (l.210)<<<<<<");
 293   //fprintf(stdout,("2*nk = " + (2*nk));
 294   //fprintf(stdout,("t1 = " + t1);
 295   //fprintf(stdout,("a  = " + a);
 296   //fprintf(stdout,("x[0] = " + x[0]);
 297   //fprintf(stdout,(">>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<");
 298
 299          //if (timers_enabled)  timer_stop(3);
 300    timer_stop(3,elapsed,start);
 301
 302 /*       Compute Gaussian deviates by acceptance-rejection method and
 303  *       tally counts in concentric square annuli.  This loop is not
 304  *       vectorizable.
 305  */
 306          //if (timers_enabled) timer_start(2);
 307     timer_start(&(start[2]));
 308          for(i=1; i<=nk;i++) {
 309             x1 = 2. * x[2*i-2] -1.0;
 310             x2 = 2. * x[2*i-1] - 1.0;
 311             t1 = x1*x1 + x2*x2;
 312             if (t1 <= 1.) {
 313                t2   = sqrt(-2. * log(t1) / t1);
 314                t3   = (x1 * t2);
 315                t4   = (x2 * t2);
 316                l    = (int)(abs(t3) > abs(t4) ? abs(t3) : abs(t4));
 317                q[l] = q[l] + 1.;
 318                sx   = sx + t3;
 319                sy   = sy + t4;
 320              }
 321     /*
 322        if(i == 1) {
 323                 fprintf(stdout,"x1 = " + x1);
 324                 fprintf(stdout,"x2 = " + x2);
 325                 fprintf(stdout,"t1 = " + t1);
 326                 fprintf(stdout,"t2 = " + t2);
 327                 fprintf(stdout,"t3 = " + t3);
 328                 fprintf(stdout,"t4 = " + t4);
 329                 fprintf(stdout,"l = " + l);
 330                 fprintf(stdout,"q[l] = " + q[l]);
 331                 fprintf(stdout,"sx = " + sx);
 332                 fprintf(stdout,"sy = " + sy);
 333        }
 334     */
 335            }
 336          //if (timers_enabled)  timer_stop(2);
 337     timer_stop(2,elapsed,start);
 338       }
 339
 340     TRACE_smpi_set_category ("finalize");
 341
 342       //int MPI_Allreduce(void *sbuf, void *rbuf, int count, MPI_Datatype dtype, MPI_Op op, MPI_Comm comm)
 343   MPI_Allreduce(&sx, x, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 344   sx = x[0]; //FIXME :  x[0] or x[1] => x[0] because fortran starts with 1
 345       MPI_Allreduce(&sy, x, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 346       sy = x[0];
 347       MPI_Allreduce(q, x, nq, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
 348
 349       for(i = 0; i < nq; i++) {
 350     q[i] = x[i];
 351   }
 352   for(i = 0; i < nq; i++) {
 353     gc += q[i];
 354   }
 355
 356   timer_stop(1,elapsed,start);
 357       tm = timer_read(1,elapsed);
 358   MPI_Allreduce(&tm, x, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
 359   tm = x[0];
 360
 361   if(node == root) {
 362     nit = 0;
 363     verified = true;
 364
 365     if(m == 24) {
 366              sx_verify_value = -3.247834652034740E3;
 367                   sy_verify_value = -6.958407078382297E3;
 368     } else if(m == 25) {
 369                 sx_verify_value = -2.863319731645753E3;
 370       sy_verify_value = -6.320053679109499E3;
 371     } else if(m == 28) {
 372             sx_verify_value = -4.295875165629892E3;
 373       sy_verify_value = -1.580732573678431E4;
 374     } else if(m == 30) {
 375             sx_verify_value =  4.033815542441498E4;
 376                   sy_verify_value = -2.660669192809235E4;
 377     } else if(m == 32) {
 378                   sx_verify_value =  4.764367927995374E4;
 379                      sy_verify_value = -8.084072988043731E4;
 380     } else if(m == 36) {
 381             sx_verify_value =  1.982481200946593E5;
 382             sy_verify_value = -1.020596636361769E5;
 383     } else {
 384       verified = false;
 385     }
 386
 387     /*
 388     fprintf(stdout,("sx        = " + sx);
 389     fprintf(stdout,("sx_verify = " + sx_verify_value);
 390     fprintf(stdout,("sy        = " + sy);
 391     fprintf(stdout,("sy_verify = " + sy_verify_value);
 392     */
 393     if(verified) {
 394       sx_err = abs((sx - sx_verify_value)/sx_verify_value);
 395       sy_err = abs((sy - sy_verify_value)/sy_verify_value);
 396       /*
 397       fprintf(stdout,("sx_err = " + sx_err);
 398       fprintf(stdout,("sy_err = " + sx_err);
 399       fprintf(stdout,("epsilon= " + epsilon);
 400       */
 401       verified = ((sx_err < epsilon) && (sy_err < epsilon));
 402     }
 403
 404     Mops = (pow(2.0, m+1))/tm/1000;
 405
 406     fprintf(stdout,"EP Benchmark Results:\n");
 407     fprintf(stdout,"CPU Time=%d\n",(int) tm);
 408     fprintf(stdout,"N = 2^%d\n",m);
 409     fprintf(stdout,"No. Gaussain Pairs =%d\n",(int) gc);
 410     fprintf(stdout,"Sum = %f %ld\n",sx,(long) sy);
 411     fprintf(stdout,"Count:");
 412     for(i = 0; i < nq; i++) {
 413       fprintf(stdout,"%d\t %ld\n",i,(long) q[i]);
 414     }
 415
 416     /*
 417     print_results("EP", _class, m+1, 0, 0, nit, npm, no_nodes, tm, Mops,
 418         "Random numbers generated", verified, npbversion,
 419         compiletime, cs1, cs2, cs3, cs4, cs5, cs6, cs7) */
 420     fprintf(stdout,"\nEP Benchmark Completed\n");
 421             fprintf(stdout,"Class           = %s\n", _class);
 422     fprintf(stdout,"Size            = %s\n", size);
 423     fprintf(stdout,"Iteration       = %d\n", nit);
 424     fprintf(stdout,"Time in seconds = %f\n",(tm/1000));
 425     fprintf(stdout,"Total processes = %d\n",no_nodes);
 426     fprintf(stdout,"Mops/s total    = %f\n",Mops);
 427     fprintf(stdout,"Mops/s/process  = %f\n", Mops/no_nodes);
 428     fprintf(stdout,"Operation type  = Random number generated\n");
 429     if(verified) {
 430       fprintf(stdout,"Verification    = SUCCESSFUL\n");
 431     } else {
 432       fprintf(stdout,"Verification    = UNSUCCESSFUL\n");
 433     }
 434              fprintf(stdout,"Total time:     %f\n",(timer_read(1,elapsed)/1000));
 435              fprintf(stdout,"Gaussian pairs: %f\n",(timer_read(2,elapsed)/1000));
 436              fprintf(stdout,"Random numbers: %f\n",(timer_read(3,elapsed)/1000));
 437          }
 438 #ifdef USE_MPE
 439     MPE_Finish_log(argv[0]);
 440 #endif
 441
 442        MPI_Finalize();
 443       }
 444
 445     int main(int argc, char **argv) {
 446        doTest(argc,argv);
 447     }