Robust/src/Benchmarks/oooJava/moldyn/JGFMolDynBench.java

   1 /**************************************************************************
   2  * * Java Grande Forum Benchmark Suite - Thread Version 1.0 * * produced by * *
   3  * Java Grande Benchmarking Project * * at * * Edinburgh Parallel Computing
   4  * Centre * * email: epcc-javagrande@epcc.ed.ac.uk * * * This version copyright
   5  * (c) The University of Edinburgh, 2001. * All rights reserved. * *
   6  **************************************************************************/
   7 public class JGFMolDynBench {
   8   public int ITERS;
   9   public double LENGTH;
  10   public double m;
  11   public double mu;
  12   public double kb;
  13   public double TSIM;
  14   public double deltat;
  15
  16   public int PARTSIZE;
  17
  18   public double[] epot;
  19   public double[] vir;
  20   public double[] ek;
  21
  22   int size, mm;
  23   int[] datasizes;
  24
  25   public int interactions;
  26   public int[] interacts;
  27
  28   public int nthreads;
  29   public int workload;
  30
  31   public JGFMolDynBench(int nthreads,int workload) {
  32     this.nthreads = nthreads;
  33     this.workload=workload;
  34   }
  35
  36   public void JGFsetsize(int size) {
  37     this.size = size;
  38   }
  39
  40   public void JGFinitialise() {
  41     interactions = 0;
  42     datasizes = new int[3];
  43     datasizes[0] = 8;
  44     datasizes[1] = 13;
  45     datasizes[2] = 11;
  46
  47     mm = datasizes[size];
  48     PARTSIZE = mm * mm * mm * 4;
  49     ITERS = 100;
  50     LENGTH = 50e-10;
  51     m = 4.0026;
  52     mu = 1.66056e-27;
  53     kb = 1.38066e-23;
  54     TSIM = 50;
  55     deltat = 5e-16;
  56   }
  57
  58   public static void JGFapplication(JGFMolDynBench mold) {
  59     double sh_force[][];
  60     double sh_force2[][][];
  61     int partsize, numthreads;
  62     partsize = mold.PARTSIZE;
  63     numthreads = mold.nthreads;
  64
  65     sh_force = new double[3][partsize];
  66     sh_force2 = new double[3][numthreads][partsize];
  67     mold.epot = new double[numthreads];
  68     mold.vir = new double[numthreads];
  69     mold.ek = new double[numthreads];
  70     mold.interacts = new int[numthreads];
  71     // for(int i=0;i<numthreads;i++) {
  72     // mold.epot[i]=new double();
  73     // mold.vir[i]=new double();
  74     // mold.ek[i]=new double();
  75     // mold.interacts[i]=new IntWrapper();
  76     // }
  77
  78     // spawn threads
  79     MDWrap[] thobjects = new MDWrap[numthreads];
  80
  81     for (int i = 0; i < numthreads; i++) {
  82       thobjects[i] = new MDWrap(new mdRunner(i, mold.mm, sh_force, sh_force2, mold.nthreads, mold,mold.workload));
  83     }
  84
  85     /*
  86      * boolean waitfordone=true; while(waitfordone) { if (mybarr.done)
  87      * waitfordone=false; }
  88      */
  89
  90     for (int i = 0; i < numthreads; i++) {
  91       // thobjects[i].md.start(mid[i]);
  92       thobjects[i].md.run();
  93     }
  94   }
  95
  96   public void JGFvalidate() {
  97     double[] refval = new double[2];
  98     refval[0] = 1731.4306625334357;
  99     refval[1] = 7397.392307839352;
 100     double dev = Math.fabs(ek[0] - refval[size]);
 101     if (dev > 1.0e-10) {
 102       // System.printString("Validation failed\n");
 103       // System.printString("Kinetic Energy = " + (long)ek[0] + "  " + (long)dev
 104       // + "  " + size + "\n");
 105     }
 106   }
 107 }
 108
 109 class mdRunner {
 110
 111   double count;
 112   int id, i, j, k, lg, mm;
 113   double l, rcoff, rcoffs, side, sideh, hsq, hsq2, vel, velt;
 114   double a, r, sum, tscale, sc, ekin, ts, sp;
 115   double den;
 116   double tref;
 117   double h;
 118   double vaver, vaverh, rand;
 119   double etot, temp, pres, rp;
 120   double u1, u2, s, xx, yy, zz;
 121   double xvelocity, yvelocity, zvelocity;
 122
 123   double[][] sh_force;
 124   double[][][] sh_force2;
 125
 126   int ijk, npartm, iseed, tint;
 127   int irep;
 128   int istop;
 129   int iprint;
 130
 131   JGFMolDynBench mymd;
 132   int nthreads;
 133   int workload;
 134
 135   public mdRunner(int id, int mm, double[][] sh_force, double[][][] sh_force2, int nthreads,
 136       JGFMolDynBench mymd, int workload) {
 137     this.id = id;
 138     this.mm = mm;
 139     this.sh_force = sh_force;
 140     this.sh_force2 = sh_force2;
 141     this.nthreads = nthreads;
 142     this.mymd = mymd;
 143     count = 0.0;
 144     den = 0.83134;
 145     tref = 0.722;
 146     h = 0.064;
 147     irep = 10;
 148     istop = 19;
 149     iprint = 10;
 150     this.workload=workload;
 151   }
 152
 153   public void init(particle[] one, int mdsize) {
 154     int id = this.id;
 155     for (int lg = 0; lg <= 1; lg++) {
 156       for (int i = 0; i < mm; i++) {
 157         for (int j = 0; j < mm; j++) {
 158           for (int k = 0; k < mm; k++) {
 159             one[ijk] =
 160                 new particle((i * a + lg * a * 0.5), (j * a + lg * a * 0.5), (k * a), xvelocity,
 161                     yvelocity, zvelocity, sh_force, sh_force2, id, one);
 162             ijk = ijk + 1;
 163           }
 164         }
 165       }
 166     }
 167
 168     for (int lg = 1; lg <= 2; lg++) {
 169       for (int i = 0; i < mm; i++) {
 170         for (int j = 0; j < mm; j++) {
 171           for (int k = 0; k < mm; k++) {
 172             one[ijk] =
 173                 new particle((i * a + (2 - lg) * a * 0.5), (j * a + (lg - 1) * a * 0.5),
 174                     (k * a + a * 0.5), xvelocity, yvelocity, zvelocity, sh_force, sh_force2, id,
 175                     one);
 176             ijk = ijk + 1;
 177           }
 178         }
 179       }
 180     }
 181
 182     /* Initialise velocities */
 183
 184     iseed = 0;
 185     double v1 = 0.0;
 186     double v2 = 0.0;
 187     random randnum = new random(iseed, v1, v2);
 188
 189     for (int i = 0; i < mdsize; i += 2) {
 190       r = randnum.seed();
 191       one[i].xvelocity = r * randnum.v1;
 192       one[i + 1].xvelocity = r * randnum.v2;
 193     }
 194
 195     for (int i = 0; i < mdsize; i += 2) {
 196       r = randnum.seed();
 197       one[i].yvelocity = r * randnum.v1;
 198       one[i + 1].yvelocity = r * randnum.v2;
 199     }
 200
 201     for (int i = 0; i < mdsize; i += 2) {
 202       r = randnum.seed();
 203       one[i].zvelocity = r * randnum.v1;
 204       one[i + 1].zvelocity = r * randnum.v2;
 205     }
 206
 207     /* velocity scaling */
 208
 209     ekin = 0.0;
 210     sp = 0.0;
 211
 212     for (int i = 0; i < mdsize; i++) {
 213       sp = sp + one[i].xvelocity;
 214     }
 215     sp = sp / mdsize;
 216
 217     for (int i = 0; i < mdsize; i++) {
 218       one[i].xvelocity = one[i].xvelocity - sp;
 219       ekin = ekin + one[i].xvelocity * one[i].xvelocity;
 220     }
 221
 222     sp = 0.0;
 223     for (int i = 0; i < mdsize; i++) {
 224       sp = sp + one[i].yvelocity;
 225     }
 226     sp = sp / mdsize;
 227
 228     for (int i = 0; i < mdsize; i++) {
 229       one[i].yvelocity = one[i].yvelocity - sp;
 230       ekin = ekin + one[i].yvelocity * one[i].yvelocity;
 231     }
 232
 233     sp = 0.0;
 234     for (int i = 0; i < mdsize; i++) {
 235       sp = sp + one[i].zvelocity;
 236     }
 237     sp = sp / mdsize;
 238
 239     for (int i = 0; i < mdsize; i++) {
 240       one[i].zvelocity = one[i].zvelocity - sp;
 241       ekin = ekin + one[i].zvelocity * one[i].zvelocity;
 242     }
 243
 244     ts = tscale * ekin;
 245     sc = h * Math.sqrt(tref / ts);
 246
 247     for (int i = 0; i < mdsize; i++) {
 248
 249       one[i].xvelocity = one[i].xvelocity * sc;
 250       one[i].yvelocity = one[i].yvelocity * sc;
 251       one[i].zvelocity = one[i].zvelocity * sc;
 252
 253     }
 254
 255   }
 256
 257   public void doinit(int mdsize) {
 258     for (int j = 0; j < 3; j++) {
 259       double[] sh = sh_force[j];
 260       for (int i = 0; i < mdsize; i++) {
 261         sh[i] = 0.0;
 262       }
 263     }
 264   }
 265
 266   public void doinit2(int mdsize) {
 267     for (int k = 0; k < 3; k++) {
 268       double[] sh = sh_force[k];
 269       double[][] sha = sh_force2[k];
 270       for (int j = 0; j < nthreads; j++) {
 271         double[] sha2 = sha[j];
 272         for (int i = 0; i < mdsize; i++) {
 273           sh[i] += sha2[i];
 274         }
 275       }
 276     }
 277
 278     for (int k = 0; k < 3; k++) {
 279       double[][] sh1 = sh_force2[k];
 280       for (int j = 0; j < nthreads; j++) {
 281         double[] sh2 = sh1[j];
 282         for (int i = 0; i < mdsize; i++) {
 283
 284           sh2[i] = 0.0;
 285         }
 286       }
 287     }
 288
 289     for (int j = 1; j < nthreads; j++) {
 290       mymd.epot[0] += mymd.epot[j];
 291       mymd.vir[0] += mymd.vir[j];
 292     }
 293     for (int j = 1; j < nthreads; j++) {
 294       mymd.epot[j] = mymd.epot[0];
 295       mymd.vir[j] = mymd.vir[0];
 296     }
 297     for (int j = 0; j < nthreads; j++) {
 298       mymd.interactions += mymd.interacts[j];
 299     }
 300
 301     for (int j = 0; j < 3; j++) {
 302       double sh[] = sh_force[j];
 303       for (int i = 0; i < mdsize; i++) {
 304         sh[i] = sh[i] * hsq2;
 305       }
 306     }
 307   }
 308
 309   public void run() {
 310     /* Parameter determination */
 311
 312     int mdsize;
 313     double tmpden;
 314     int movemx = 50;
 315     particle[] one;
 316     int id;
 317     id = this.id;
 318     mdsize = mymd.PARTSIZE;
 319     one = new particle[mdsize];
 320     l = mymd.LENGTH;
 321     tmpden = den;
 322     side = Math.pow((mdsize / tmpden), 0.3333333);
 323     rcoff = mm / 4.0;
 324
 325     a = side / mm;
 326     sideh = side * 0.5;
 327     hsq = h * h;
 328     hsq2 = hsq * 0.5;
 329     npartm = mdsize - 1;
 330     rcoffs = rcoff * rcoff;
 331     tscale = 16.0 / (1.0 * mdsize - 1.0);
 332     vaver = 1.13 * Math.sqrt(tref / 24.0);
 333     vaverh = vaver * h;
 334
 335     /* Particle Generation */
 336
 337     xvelocity = 0.0;
 338     yvelocity = 0.0;
 339     zvelocity = 0.0;
 340     ijk = 0;
 341     init(one, mdsize);
 342
 343     /* Synchronise threads and start timer before MD simulation */
 344
 345     /* MD simulation */
 346
 347     int numP= (mdsize / workload)+1;
 348
 349     double scratchpad[][][];
 350     scratchpad=new double[numP][3][mdsize];
 351
 352     for (int move = 0; move < movemx; move++) {
 353       /* move the particles and update velocities */
 354       for (int i = 0; i < mdsize; i++) {
 355         one[i].domove(side, i);
 356       }
 357
 358       if (id == 0) {
 359         doinit(mdsize);
 360       }
 361
 362       mymd.epot[id] = 0.0;
 363       mymd.vir[id] = 0.0;
 364       mymd.interacts[id] = 0;
 365
 366       /* compute forces */
 367       int numThread = nthreads;
 368       int lworkload = workload;
 369       // for (int i=0+id;i<mdsize;i+=numThread) {
 370       int scratch_idx=0;
 371       for (int i = 0 ; i < mdsize; i += lworkload) {
 372
 373         int ilow = i;
 374         int iupper = i + lworkload;
 375         if (iupper > mdsize) {
 376           iupper = mdsize;
 377         }
 378         int l_size = iupper - ilow;
 379
 380         double workingpad[][]=scratchpad[scratch_idx++];
 381         sese parallel{
 382           for(int j=0;j<3;j++){
 383             for(int l=0;l<mdsize;l++){
 384               workingpad[j][l]=0;
 385             }
 386           }
 387           MDStore store = new MDStore();
 388           for(int idx=ilow;idx<iupper;idx++){
 389 //            one[i].force(side, rcoff, mdsize, i, xx, yy, zz, mymd, worker);
 390             one[idx].force(side, rcoff, mdsize, idx, xx, yy, zz, mymd, store,workingpad);
 391           }
 392         }
 393
 394         sese serial{
 395           for (int k = 0; k < 3; k++) {
 396             for (int j = 0; j < mdsize; j++) {
 397 //              sh_force[k][j] += worker.sh_force2[k][j];
 398               sh_force[k][j] += workingpad[k][j];
 399             }
 400           }
 401           mymd.epot[0] += store.epot;
 402           mymd.vir[0] += store.vir;
 403           mymd.interactions += store.interacts;
 404         }
 405
 406       }
 407
 408       for (int k = 0; k < 3; k++) {
 409         for (int j = 0; j < mdsize; j++) {
 410           sh_force[k][j] = sh_force[k][j] * hsq2;
 411         }
 412       }
 413
 414       /* update force arrays */
 415       // if(id == 0) {
 416       // doinit2(mdsize);
 417       // }
 418
 419       /* scale forces, update velocities */
 420       sum = 0.0;
 421       for (int i = 0; i < mdsize; i++) {
 422         sum = sum + one[i].mkekin(hsq2, i);
 423       }
 424
 425       ekin = sum / hsq;
 426
 427       vel = 0.0;
 428       count = 0.0;
 429
 430       /* average velocity */
 431
 432       for (int i = 0; i < mdsize; i++) {
 433         velt = one[i].velavg(vaverh, h);
 434         if (velt > vaverh) {
 435           count = count + 1.0;
 436         }
 437         vel = vel + velt;
 438       }
 439
 440       vel = vel / h;
 441
 442       /* temperature scale if required */
 443
 444       if ((move < istop) && (((move + 1) % irep) == 0)) {
 445         sc = Math.sqrt(tref / (tscale * ekin));
 446         for (int i = 0; i < mdsize; i++) {
 447           one[i].dscal(sc, 1);
 448         }
 449         ekin = tref / tscale;
 450       }
 451
 452       /* sum to get full potential energy and virial */
 453
 454       if (((move + 1) % iprint) == 0) {
 455         mymd.ek[id] = 24.0 * ekin;
 456         mymd.epot[id] = 4.0 * mymd.epot[id];
 457         etot = mymd.ek[id] + mymd.epot[id];
 458         temp = tscale * ekin;
 459         pres = tmpden * 16.0 * (ekin - mymd.vir[id]) / mdsize;
 460         vel = vel / mdsize;
 461         rp = (count / mdsize) * 100.0;
 462       }
 463
 464       // if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run",
 465       // instr.timers);
 466     }
 467   }
 468 }
 469
 470 class particle {
 471
 472   public double xcoord, ycoord, zcoord;
 473   public double xvelocity, yvelocity, zvelocity;
 474   int part_id;
 475   int id;
 476   double[][] sh_force;
 477   double[][][] sh_force2;
 478   particle[] one;
 479
 480   public particle(double xcoord, double ycoord, double zcoord, double xvelocity, double yvelocity,
 481       double zvelocity, double[][] sh_force, double[][][] sh_force2, int id, particle[] one) {
 482
 483     this.xcoord = xcoord;
 484     this.ycoord = ycoord;
 485     this.zcoord = zcoord;
 486     this.xvelocity = xvelocity;
 487     this.yvelocity = yvelocity;
 488     this.zvelocity = zvelocity;
 489     this.sh_force = sh_force;
 490     this.sh_force2 = sh_force2;
 491     this.id = id;
 492     this.one = one;
 493   }
 494
 495   public void domove(double side, int part_id) {
 496
 497     xcoord = xcoord + xvelocity + sh_force[0][part_id];
 498     ycoord = ycoord + yvelocity + sh_force[1][part_id];
 499     zcoord = zcoord + zvelocity + sh_force[2][part_id];
 500
 501     if (xcoord < 0) {
 502       xcoord = xcoord + side;
 503     }
 504     if (xcoord > side) {
 505       xcoord = xcoord - side;
 506     }
 507     if (ycoord < 0) {
 508       ycoord = ycoord + side;
 509     }
 510     if (ycoord > side) {
 511       ycoord = ycoord - side;
 512     }
 513     if (zcoord < 0) {
 514       zcoord = zcoord + side;
 515     }
 516     if (zcoord > side) {
 517       zcoord = zcoord - side;
 518     }
 519
 520     xvelocity = xvelocity + sh_force[0][part_id];
 521     yvelocity = yvelocity + sh_force[1][part_id];
 522     zvelocity = zvelocity + sh_force[2][part_id];
 523
 524   }
 525
 526   // public void force(double side, double rcoff,int mdsize,int x, double xx,
 527   // double yy, double zz, JGFMolDynBench mymd) {
 528   public void force(double side, double rcoff, int mdsize, int x, double xx, double yy, double zz,
 529       JGFMolDynBench mymd, MDStore store,double workingpad[][]) {
 530
 531     double sideh;
 532     double rcoffs;
 533
 534     double fxi, fyi, fzi;
 535     double rd, rrd, rrd2, rrd3, rrd4, rrd6, rrd7, r148;
 536     double forcex, forcey, forcez;
 537     int id = this.id;
 538     sideh = 0.5 * side;
 539     rcoffs = rcoff * rcoff;
 540
 541     fxi = 0.0;
 542     fyi = 0.0;
 543     fzi = 0.0;
 544
 545     for (int i = x + 1; i < mdsize; i++) {
 546       xx = this.xcoord - one[i].xcoord;
 547       yy = this.ycoord - one[i].ycoord;
 548       zz = this.zcoord - one[i].zcoord;
 549
 550       if (xx < (-sideh)) {
 551         xx = xx + side;
 552       }
 553       if (xx > (sideh)) {
 554         xx = xx - side;
 555       }
 556       if (yy < (-sideh)) {
 557         yy = yy + side;
 558       }
 559       if (yy > (sideh)) {
 560         yy = yy - side;
 561       }
 562       if (zz < (-sideh)) {
 563         zz = zz + side;
 564       }
 565       if (zz > (sideh)) {
 566         zz = zz - side;
 567       }
 568
 569       rd = xx * xx + yy * yy + zz * zz;
 570
 571       if (rd <= rcoffs) {
 572         rrd = 1.0 / rd;
 573         rrd2 = rrd * rrd;
 574         rrd3 = rrd2 * rrd;
 575         rrd4 = rrd2 * rrd2;
 576         rrd6 = rrd2 * rrd4;
 577         rrd7 = rrd6 * rrd;
 578         // mymd.epot[id] += (rrd6 - rrd3);
 579         store.epot += (rrd6 - rrd3);
 580         r148 = rrd7 - 0.5 * rrd4;
 581         // mymd.vir[id] += - rd*r148;
 582         store.vir += -rd * r148;
 583         forcex = xx * r148;
 584         fxi = fxi + forcex;
 585
 586         // sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
 587 //        worker.sh_force2[0][i] = worker.sh_force2[0][i] - forcex;
 588         workingpad[0][i] = workingpad[0][i] - forcex;
 589
 590         forcey = yy * r148;
 591         fyi = fyi + forcey;
 592
 593         // sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
 594 //        worker.sh_force2[1][i] = worker.sh_force2[1][i] - forcey;
 595         workingpad[1][i] = workingpad[1][i] - forcey;
 596
 597         forcez = zz * r148;
 598         fzi = fzi + forcez;
 599
 600         // sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
 601 //        worker.sh_force2[2][i] = worker.sh_force2[2][i] - forcez;
 602         workingpad[2][i] = workingpad[2][i] - forcez;
 603
 604         // mymd.interacts[id]++;
 605         store.interacts++;
 606       }
 607
 608     }
 609
 610     // sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
 611     // sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
 612     // sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
 613
 614 //    worker.sh_force2[0][x] = worker.sh_force2[0][x] + fxi;
 615 //    worker.sh_force2[1][x] = worker.sh_force2[1][x] + fyi;
 616 //    worker.sh_force2[2][x] = worker.sh_force2[2][x] + fzi;
 617
 618     workingpad[0][x] = workingpad[0][x] + fxi;
 619     workingpad[1][x] = workingpad[1][x] + fyi;
 620     workingpad[2][x] = workingpad[2][x] + fzi;
 621
 622   }
 623
 624   public double mkekin(double hsq2, int part_id) {
 625
 626     double sumt = 0.0;
 627
 628     xvelocity = xvelocity + sh_force[0][part_id];
 629     yvelocity = yvelocity + sh_force[1][part_id];
 630     zvelocity = zvelocity + sh_force[2][part_id];
 631
 632     sumt = (xvelocity * xvelocity) + (yvelocity * yvelocity) + (zvelocity * zvelocity);
 633     return sumt;
 634   }
 635
 636   public double velavg(double vaverh, double h) {
 637
 638     double velt;
 639     double sq;
 640
 641     sq = Math.sqrt(xvelocity * xvelocity + yvelocity * yvelocity + zvelocity * zvelocity);
 642
 643     velt = sq;
 644     return velt;
 645   }
 646
 647   public void dscal(double sc, int incx) {
 648     xvelocity = xvelocity * sc;
 649     yvelocity = yvelocity * sc;
 650     zvelocity = zvelocity * sc;
 651   }
 652 }
 653
 654 class random {
 655
 656   public int iseed;
 657   public double v1, v2;
 658
 659   public random(int iseed, double v1, double v2) {
 660     this.iseed = iseed;
 661     this.v1 = v1;
 662     this.v2 = v2;
 663   }
 664
 665   public double update() {
 666
 667     double rand;
 668     double scale = 4.656612875e-10;
 669
 670     int is1, is2, iss2;
 671     int imult = 16807;
 672     int imod = 2147483647;
 673
 674     if (iseed <= 0) {
 675       iseed = 1;
 676     }
 677
 678     is2 = iseed % 32768;
 679     is1 = (iseed - is2) / 32768;
 680     iss2 = is2 * imult;
 681     is2 = iss2 % 32768;
 682     is1 = (is1 * imult + (iss2 - is2) / 32768) % (65536);
 683
 684     iseed = (is1 * 32768 + is2) % imod;
 685
 686     rand = scale * iseed;
 687
 688     return rand;
 689
 690   }
 691
 692   public double seed() {
 693
 694     double s, u1, u2, r;
 695     s = 1.0;
 696     do {
 697       u1 = update();
 698       u2 = update();
 699
 700       v1 = 2.0 * u1 - 1.0;
 701       v2 = 2.0 * u2 - 1.0;
 702       s = v1 * v1 + v2 * v2;
 703
 704     } while (s >= 1.0);
 705
 706     r = Math.sqrt(-2.0 * Math.log(s) / s);
 707
 708     return r;
 709
 710   }
 711 }