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