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