Robust/src/Benchmarks/Prefetch/Moldyn/java/JGFMolDynBench.java

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