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