micro/bh.d

   1 /**\r
   2 A D implementation of the _bh_ Olden benchmark.\r
   3 The Olden benchmark implements the Barnes-Hut benchmark\r
   4 that is decribed in:\r
   5 \r
   6 J. Barnes and P. Hut, "A hierarchical o(N log N) force-calculation algorithm",\r
   7 Nature, 324:446-449, Dec. 1986\r
   8 \r
   9 The original code in the Olden benchmark suite is derived from the\r
  10 ftp://hubble.ifa.hawaii.edu/pub/barnes/treecode\r
  11 source distributed by Barnes.\r
  12 \r
  13 Java code converted to D by leonardo maffi, V.1.0, Dec 25 2009.\r
  14 \r
  15 Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.\r
  16 Downloaded from http://www.fantascienza.net/leonardo/js/\r
  17                 http://www.fantascienza.net/leonardo/js/dolden_bh.zip\r
  18 */\r
  19 \r
  20 \r
  21 version (Tango) {\r
  22     import tango.stdc.stdio: printf, sprintf, fprintf, stderr;\r
  23     import tango.stdc.stdlib: exit, malloc, free;\r
  24     import tango.stdc.time: CLOCKS_PER_SEC, clock;\r
  25     import tango.math.Math: sqrt, floor, PI, pow;\r
  26 \r
  27     import Integer = tango.text.convert.Integer;\r
  28     alias Integer.parse toInt;\r
  29 } else {\r
  30     import std.c.stdio: printf, sprintf, fprintf, stderr;\r
  31     import std.c.stdlib: exit, malloc, free;\r
  32     import std.c.time: CLOCKS_PER_SEC, clock;\r
  33     import std.math: sqrt, floor, PI, pow;\r
  34 \r
  35     version (D_Version2) {\r
  36         import std.conv: to;\r
  37         alias to!(int, char[]) toInt;\r
  38     } else {\r
  39         import std.conv: toInt;\r
  40     }\r
  41 }\r
  42 \r
  43 \r
  44 double myclock() {\r
  45     return clock() / cast(double)CLOCKS_PER_SEC;\r
  46 }\r
  47 \r
  48 \r
  49 /*\r
  50 Basic uniform random generator: Minimal Standard in Park and\r
  51 Miller (1988): "Random Number Generators: Good Ones Are Hard to\r
  52 Find", Comm. of the ACM, 31, 1192-1201.\r
  53 Parameters: m = 2^31-1, a=48271.\r
  54 \r
  55 Adapted from Pascal code by Jesper Lund:\r
  56 http://www.gnu-pascal.de/crystal/gpc/en/mail1390.html\r
  57 */\r
  58 struct Random {\r
  59     const int m = int.max;\r
  60     const int a = 48_271;\r
  61     const int q = m / a;\r
  62     const int r = m % a;\r
  63     int seed;\r
  64 \r
  65     public double uniform(double min, double max) {\r
  66         int k = seed / q;\r
  67         seed = a * (seed - k * q) - r * k;\r
  68         if (seed < 1)\r
  69             seed += m;\r
  70         double r = cast(double)seed / m;\r
  71         return r * (max - min) + min;\r
  72     }\r
  73 }\r
  74 \r
  75 \r
  76 interface Enumeration {\r
  77     bool hasMoreElements();\r
  78     Object nextElement();\r
  79 }\r
  80 \r
  81 \r
  82 /**\r
  83 A class representing a three dimensional vector that implements\r
  84 several math operations.  To improve speed we implement the\r
  85 vector as an array of doubles rather than use the exising\r
  86 code in the java.util.Vec3 class.\r
  87 */\r
  88 final class Vec3 {\r
  89     /// The number of dimensions in the vector\r
  90     // Can't be changed because of crossProduct()\r
  91     const int NDIM = 3;\r
  92 \r
  93     /// An array containing the values in the vector.\r
  94     private double data[];\r
  95 \r
  96     /// Construct an empty 3 dimensional vector for use in Barnes-Hut algorithm.\r
  97     this() {\r
  98         data.length = NDIM;\r
  99         data[] = 0.0;\r
 100     }\r
 101 \r
 102     /// Create a copy of the vector. Return a clone of the math vector\r
 103     public Vec3 clone() {\r
 104         Vec3 v = new Vec3;\r
 105         v.data.length = NDIM;\r
 106         v.data[] = this.data[];\r
 107         return v;\r
 108     }\r
 109 \r
 110     /**\r
 111     Return the value at the i'th index of the vector.\r
 112     @param i the vector index\r
 113     @return the value at the i'th index of the vector.\r
 114     */\r
 115     final double value(int i) {\r
 116         return data[i];\r
 117     }\r
 118 \r
 119     /**\r
 120     Set the value of the i'th index of the vector.\r
 121     @param i the vector index\r
 122     @param v the value to store\r
 123     */\r
 124     final void value(int i, double v) {\r
 125         data[i] = v;\r
 126     }\r
 127 \r
 128     /**\r
 129     Set one of the dimensions of the vector to 1.0\r
 130     param j the dimension to set.\r
 131     */\r
 132     final void unit(int j) {\r
 133         for (int i = 0; i < NDIM; i++)\r
 134             data[i] = (i == j) ? 1.0 : 0.0;\r
 135     }\r
 136 \r
 137     /**\r
 138     Add two vectors and the result is placed in this vector.\r
 139     @param u the other operand of the addition\r
 140     */\r
 141     final void addition(Vec3 u) {\r
 142         for (int i = 0; i < NDIM; i++)\r
 143             data[i] += u.data[i];\r
 144     }\r
 145 \r
 146     /**\r
 147    * Subtract two vectors and the result is placed in this vector.\r
 148    * This vector contain the first operand.\r
 149    * @param u the other operand of the subtraction.\r
 150    **/\r
 151     final void subtraction(Vec3 u) {\r
 152         for (int i = 0; i < NDIM; i++)\r
 153             data[i] -= u.data[i];\r
 154     }\r
 155 \r
 156     /**\r
 157     Subtract two vectors and the result is placed in this vector.\r
 158     @param u the first operand of the subtraction.\r
 159     @param v the second opernd of the subtraction\r
 160     */\r
 161     final void subtraction(Vec3 u, Vec3 v) {\r
 162         for (int i = 0; i < NDIM; i++)\r
 163             data[i] = u.data[i] - v.data[i];\r
 164     }\r
 165 \r
 166     /**\r
 167     Multiply the vector times a scalar.\r
 168     @param s the scalar value\r
 169     **/\r
 170     final void multScalar(double s) {\r
 171         for (int i = 0; i < NDIM; i++)\r
 172             data[i] *= s;\r
 173     }\r
 174 \r
 175     /**\r
 176     * Multiply the vector times a scalar and place the result in this vector.\r
 177     * @param u the vector\r
 178     * @param s the scalar value\r
 179     **/\r
 180     final void multScalar(Vec3 u, double s) {\r
 181         for (int i = 0; i < NDIM; i++)\r
 182             data[i] = u.data[i] * s;\r
 183     }\r
 184 \r
 185     /**\r
 186     * Divide each element of the vector by a scalar value.\r
 187     * @param s the scalar value.\r
 188     **/\r
 189     final void divScalar(double s) {\r
 190         for (int i = 0; i < NDIM; i++)\r
 191             data[i] /= s;\r
 192     }\r
 193 \r
 194     /**\r
 195     * Return the dot product of a vector.\r
 196     * @return the dot product of a vector.\r
 197     **/\r
 198     final double dotProduct() {\r
 199         double s = 0.0;\r
 200         for (int i = 0; i < NDIM; i++)\r
 201             s += data[i] * data[i];\r
 202         return s;\r
 203     }\r
 204 \r
 205     final double absolute() {\r
 206         double tmp = 0.0;\r
 207         for (int i = 0; i < NDIM; i++)\r
 208             tmp += data[i] * data[i];\r
 209         return sqrt(tmp);\r
 210     }\r
 211 \r
 212     final double distance(Vec3 v) {\r
 213         double tmp = 0.0;\r
 214         for (int i = 0; i < NDIM; i++)\r
 215             tmp += ((data[i] - v.data[i]) * (data[i] - v.data[i]));\r
 216         return sqrt(tmp);\r
 217     }\r
 218 \r
 219     final void crossProduct(Vec3 u, Vec3 w) {\r
 220         data[0] = u.data[1] * w.data[2] - u.data[2] * w.data[1];\r
 221         data[1] = u.data[2] * w.data[0] - u.data[0] * w.data[2];\r
 222         data[2] = u.data[0] * w.data[1] - u.data[1] * w.data[0];\r
 223     }\r
 224 \r
 225     final void incrementalAdd(Vec3 u) {\r
 226         for (int i = 0; i < NDIM; i++)\r
 227             data[i] += u.data[i];\r
 228     }\r
 229 \r
 230     final void incrementalSub(Vec3 u) {\r
 231         for (int i = 0; i < NDIM; i++)\r
 232             data[i] -= u.data[i];\r
 233     }\r
 234 \r
 235     final void incrementalMultScalar(double s) {\r
 236         for (int i = 0; i < NDIM; i++)\r
 237             data[i] *= s;\r
 238     }\r
 239 \r
 240     final void incrementalDivScalar(double s) {\r
 241         for (int i = 0; i < NDIM; i++)\r
 242             data[i] /= s;\r
 243     }\r
 244 \r
 245     final void addScalar(Vec3 u, double s) {\r
 246         for (int i = 0; i < NDIM; i++)\r
 247             data[i] = u.data[i] + s;\r
 248     }\r
 249 \r
 250     public char* toCString() {\r
 251         // this is not generic code at all\r
 252         char* result = cast(char*)malloc(100);\r
 253         if (result == null) exit(1);\r
 254         sprintf(result, "%.17f %.17f %.17f ", data[0], data[1], data[2]);\r
 255         return result;\r
 256     }\r
 257 }\r
 258 \r
 259 \r
 260 /// A class that represents the common fields of a cell or body data structure.\r
 261 abstract class Node {\r
 262     // mass of the node\r
 263     double mass;\r
 264 \r
 265     // Position of the node\r
 266     Vec3 pos;\r
 267 \r
 268     // highest bit of int coord\r
 269     const int IMAX = 1073741824;\r
 270 \r
 271     // potential softening parameter\r
 272     const double EPS = 0.05;\r
 273 \r
 274     /// Construct an empty node\r
 275     protected this() {\r
 276         mass = 0.0;\r
 277         pos = new Vec3();\r
 278     }\r
 279 \r
 280     abstract Node loadTree(Body p, Vec3 xpic, int l, Tree root);\r
 281 \r
 282     abstract double hackcofm();\r
 283 \r
 284     abstract HG walkSubTree(double dsq, HG hg);\r
 285 \r
 286     static final int oldSubindex(Vec3 ic, int l) {\r
 287         int i = 0;\r
 288         for (int k = 0; k < Vec3.NDIM; k++)\r
 289             if ((cast(int)ic.value(k) & l) != 0)\r
 290                 i += Cell.NSUB >> (k + 1);\r
 291         return i;\r
 292     }\r
 293 \r
 294     public char* toCString() {\r
 295         char* result = cast(char*)malloc(130);\r
 296         if (result == null) exit(1);\r
 297         char* pos_str = pos.toCString();\r
 298         sprintf(result, "%f : %s", mass, pos_str);\r
 299         free(pos_str);\r
 300         return result;\r
 301     }\r
 302 \r
 303     /// Compute a single body-body or body-cell interaction\r
 304     final HG gravSub(HG hg) {\r
 305         Vec3 dr = new Vec3();\r
 306         dr.subtraction(pos, hg.pos0);\r
 307 \r
 308         double drsq = dr.dotProduct() + (EPS * EPS);\r
 309         double drabs = sqrt(drsq);\r
 310 \r
 311         double phii = mass / drabs;\r
 312         hg.phi0 -= phii;\r
 313         double mor3 = phii / drsq;\r
 314         dr.multScalar(mor3);\r
 315         hg.acc0.addition(dr);\r
 316         return hg;\r
 317     }\r
 318 \r
 319     /**\r
 320     * A sub class which is used to compute and save information during the\r
 321     * gravity computation phase.\r
 322     **/\r
 323     static protected final class HG {\r
 324         /// Body to skip in force evaluation\r
 325         Body pskip;\r
 326 \r
 327         /// Point at which to evaluate field\r
 328         Vec3 pos0;\r
 329 \r
 330         /// Computed potential at pos0\r
 331 \r
 332         double phi0;\r
 333 \r
 334         /// computed acceleration at pos0\r
 335         Vec3 acc0;\r
 336 \r
 337         /**\r
 338         * Create a HG  object.\r
 339         * @param b the body object\r
 340         * @param p a vector that represents the body\r
 341         **/\r
 342         this(Body b, Vec3 p) {\r
 343             pskip = b;\r
 344             pos0 = p.clone();\r
 345             phi0 = 0.0;\r
 346             acc0 = new Vec3();\r
 347         }\r
 348     }\r
 349 }\r
 350 \r
 351 \r
 352 /// A class used to representing particles in the N-body simulation.\r
 353 final class Body : Node {\r
 354     Vec3 vel, acc, newAcc;\r
 355     double phi = 0.0;\r
 356     private Body next, procNext;\r
 357 \r
 358     /// Create an empty body.\r
 359     this() {\r
 360         vel = new Vec3();\r
 361         acc = new Vec3();\r
 362         newAcc = new Vec3();\r
 363     }\r
 364 \r
 365     /**\r
 366     * Set the next body in the list.\r
 367     * @param n the body\r
 368     **/\r
 369     final void setNext(Body n) {\r
 370         next = n;\r
 371     }\r
 372 \r
 373     /**\r
 374     * Get the next body in the list.\r
 375     * @return the next body\r
 376     **/\r
 377     final Body getNext() {\r
 378         return next;\r
 379     }\r
 380 \r
 381     /**\r
 382     * Set the next body in the list.\r
 383     * @param n the body\r
 384     **/\r
 385     final void setProcNext(Body n) {\r
 386         procNext = n;\r
 387     }\r
 388 \r
 389     /**\r
 390     * Get the next body in the list.\r
 391     * @return the next body\r
 392     **/\r
 393     final Body getProcNext() {\r
 394         return procNext;\r
 395     }\r
 396 \r
 397     /**\r
 398     * Enlarge cubical "box", salvaging existing tree structure.\r
 399     * @param tree the root of the tree.\r
 400     * @param nsteps the current time step\r
 401     **/\r
 402     final void expandBox(Tree tree, int nsteps) {\r
 403         Vec3 rmid = new Vec3();\r
 404 \r
 405         bool inbox = icTest(tree);\r
 406         while (!inbox) {\r
 407             double rsize = tree.rsize;\r
 408             rmid.addScalar(tree.rmin, 0.5 * rsize);\r
 409 \r
 410             for (int k = 0; k < Vec3.NDIM; k++)\r
 411                 if (pos.value(k) < rmid.value(k)) {\r
 412                     double rmin = tree.rmin.value(k);\r
 413                     tree.rmin.value(k, rmin - rsize);\r
 414                 }\r
 415 \r
 416             tree.rsize = 2.0 * rsize;\r
 417             if (tree.root !is null) {\r
 418                 Vec3 ic = tree.intcoord(rmid);\r
 419                 if (ic is null)\r
 420                     throw new Exception("Value is out of bounds");\r
 421                 int k = oldSubindex(ic, IMAX >> 1);\r
 422                 Cell newt = new Cell();\r
 423                 newt.subp[k] = tree.root;\r
 424                 tree.root = newt;\r
 425                 inbox = icTest(tree);\r
 426             }\r
 427         }\r
 428     }\r
 429 \r
 430     /**\r
 431     * Check the bounds of the body and return true if it isn't in the\r
 432     * correct bounds.\r
 433     **/\r
 434     final bool icTest(Tree tree) {\r
 435         double pos0 = pos.value(0);\r
 436         double pos1 = pos.value(1);\r
 437         double pos2 = pos.value(2);\r
 438 \r
 439         // by default, it is in bounds\r
 440         bool result = true;\r
 441 \r
 442         double xsc = (pos0 - tree.rmin.value(0)) / tree.rsize;\r
 443         if (!(0.0 < xsc && xsc < 1.0))\r
 444             result = false;\r
 445 \r
 446         xsc = (pos1 - tree.rmin.value(1)) / tree.rsize;\r
 447         if (!(0.0 < xsc && xsc < 1.0))\r
 448             result = false;\r
 449 \r
 450         xsc = (pos2 - tree.rmin.value(2)) / tree.rsize;\r
 451         if (!(0.0 < xsc && xsc < 1.0))\r
 452             result = false;\r
 453 \r
 454         return result;\r
 455     }\r
 456 \r
 457     /**\r
 458     * Descend Tree and insert particle.  We're at a body so we need to\r
 459     * create a cell and attach this body to the cell.\r
 460     * @param p the body to insert\r
 461     * @param xpic\r
 462     * @param l\r
 463     * @param tree the root of the data structure\r
 464     * @return the subtree with the new body inserted\r
 465     **/\r
 466     final Node loadTree(Body p, Vec3 xpic, int l, Tree tree) {\r
 467         // create a Cell\r
 468         Cell retval = new Cell();\r
 469         int si = subindex(tree, l);\r
 470         // attach this Body node to the cell\r
 471         retval.subp[si] = this;\r
 472 \r
 473         // move down one level\r
 474         si = oldSubindex(xpic, l);\r
 475         Node rt = retval.subp[si];\r
 476         if (rt !is null)\r
 477             retval.subp[si] = rt.loadTree(p, xpic, l >> 1, tree);\r
 478         else\r
 479             retval.subp[si] = p;\r
 480         return retval;\r
 481     }\r
 482 \r
 483     /**\r
 484     * Descend tree finding center of mass coordinates\r
 485     * @return the mass of this node\r
 486     **/\r
 487     final double hackcofm() {\r
 488         return mass;\r
 489     }\r
 490 \r
 491     /// iteration of the bodies\r
 492     int opApply(int delegate(ref Body) dg) {\r
 493         int result;\r
 494         Body current = this;\r
 495         while (current !is null) {\r
 496             result = dg(current);\r
 497             current = current.next;\r
 498             if (result)\r
 499                 break;\r
 500         }\r
 501         return result;\r
 502     }\r
 503 \r
 504     /// inverse iteration of the bodies\r
 505     int opApplyReverse(int delegate(ref Body) dg) {\r
 506         int result;\r
 507         Body current = this;\r
 508         while (current !is null) {\r
 509             result = dg(current);\r
 510             current = current.procNext;\r
 511             if (result)\r
 512                 break;\r
 513         }\r
 514         return result;\r
 515     }\r
 516 \r
 517 \r
 518     /**\r
 519     * Return an enumeration of the bodies\r
 520     * @return an enumeration of the bodies\r
 521     **/\r
 522     final Enumeration elements() {\r
 523         // a local class that implements the enumerator\r
 524         static final class Enumerate : Enumeration {\r
 525             private Body current;\r
 526             public this(Body outer) {\r
 527                 //this.current = Body.this;\r
 528                 this.current = outer;\r
 529             }\r
 530             public bool hasMoreElements() {\r
 531                 return current !is null;\r
 532             }\r
 533             public Object nextElement() {\r
 534                 Object retval = current;\r
 535                 current = current.next;\r
 536                 return retval;\r
 537             }\r
 538         }\r
 539 \r
 540         return new Enumerate(this);\r
 541     }\r
 542 \r
 543     final Enumeration elementsRev() {\r
 544         // a local class that implements the enumerator\r
 545         static class Enumerate : Enumeration {\r
 546             private Body current;\r
 547             public this(Body outer) {\r
 548                 //this.current = Body.this;\r
 549                 this.current = outer;\r
 550             }\r
 551             public bool hasMoreElements() {\r
 552                 return current !is null;\r
 553             }\r
 554             public Object nextElement() {\r
 555                 Object retval = current;\r
 556                 current = current.procNext;\r
 557                 return retval;\r
 558             }\r
 559         }\r
 560 \r
 561         return new Enumerate(this);\r
 562     }\r
 563 \r
 564     /**\r
 565     * Determine which subcell to select.\r
 566     * Combination of intcoord and oldSubindex.\r
 567     * @param t the root of the tree\r
 568     **/\r
 569     final int subindex(Tree tree, int l) {\r
 570         Vec3 xp = new Vec3();\r
 571 \r
 572         double xsc = (pos.value(0) - tree.rmin.value(0)) / tree.rsize;\r
 573         xp.value(0, floor(IMAX * xsc));\r
 574 \r
 575         xsc = (pos.value(1) - tree.rmin.value(1)) / tree.rsize;\r
 576         xp.value(1, floor(IMAX * xsc));\r
 577 \r
 578         xsc = (pos.value(2) - tree.rmin.value(2)) / tree.rsize;\r
 579         xp.value(2, floor(IMAX * xsc));\r
 580 \r
 581         int i = 0;\r
 582         for (int k = 0; k < Vec3.NDIM; k++)\r
 583             if ((cast(int)xp.value(k) & l) != 0)\r
 584                 i += Cell.NSUB >> (k + 1);\r
 585         return i;\r
 586     }\r
 587 \r
 588     /**\r
 589     * Evaluate gravitational field on the body.\r
 590     * The original olden version calls a routine named "walkscan",\r
 591     * but we use the same name that is in the Barnes code.\r
 592     **/\r
 593     final void hackGravity(double rsize, Node root) {\r
 594         Vec3 pos0 = pos.clone();\r
 595         HG hg = new HG(this, pos);\r
 596         hg = root.walkSubTree(rsize * rsize, hg);\r
 597         phi = hg.phi0;\r
 598         newAcc = hg.acc0;\r
 599     }\r
 600 \r
 601     /// Recursively walk the tree to do hackwalk calculation\r
 602     final HG walkSubTree(double dsq, HG hg) {\r
 603         if (this != hg.pskip)\r
 604             hg = gravSub(hg);\r
 605         return hg;\r
 606     }\r
 607 \r
 608     /**\r
 609     * Return a string represenation of a body.\r
 610     * @return a string represenation of a body.\r
 611     **/\r
 612     public char* toCString() {\r
 613         char* result = cast(char*)malloc(140);\r
 614         if (result == null) exit(1);\r
 615         char* super_str = super.toCString();\r
 616         sprintf(result, "Body %s", super_str);\r
 617         free(super_str);\r
 618         return result;\r
 619     }\r
 620 }\r
 621 \r
 622 \r
 623 \r
 624 /// A class used to represent internal nodes in the tree\r
 625 final class Cell : Node {\r
 626     // subcells per cell\r
 627     const int NSUB = 1 << Vec3.NDIM;\r
 628 \r
 629     /**\r
 630     * The children of this cell node.  Each entry may contain either\r
 631     * another cell or a body.\r
 632     **/\r
 633     Node[] subp;\r
 634     Cell next;\r
 635 \r
 636     this() {\r
 637         subp.length = NSUB;\r
 638     }\r
 639 \r
 640     /**\r
 641     * Descend Tree and insert particle.  We're at a cell so\r
 642     * we need to move down the tree.\r
 643     * @param p the body to insert into the tree\r
 644     * @param xpic\r
 645     * @param l\r
 646     * @param tree the root of the tree\r
 647     * @return the subtree with the new body inserted\r
 648     **/\r
 649     Node loadTree(Body p, Vec3 xpic, int l, Tree tree) {\r
 650         // move down one level\r
 651         int si = oldSubindex(xpic, l);\r
 652         Node rt = subp[si];\r
 653         if (rt !is null)\r
 654             subp[si] = rt.loadTree(p, xpic, l >> 1, tree);\r
 655         else\r
 656             subp[si] = p;\r
 657         return this;\r
 658     }\r
 659 \r
 660     /**\r
 661     * Descend tree finding center of mass coordinates\r
 662     * @return the mass of this node\r
 663     **/\r
 664     double hackcofm() {\r
 665         double mq = 0.0;\r
 666         Vec3 tmpPos = new Vec3();\r
 667         Vec3 tmpv = new Vec3();\r
 668         for (int i = 0; i < NSUB; i++) {\r
 669             Node r = subp[i];\r
 670             if (r !is null) {\r
 671                 double mr = r.hackcofm();\r
 672                 mq = mr + mq;\r
 673                 tmpv.multScalar(r.pos, mr);\r
 674                 tmpPos.addition(tmpv);\r
 675             }\r
 676         }\r
 677         mass = mq;\r
 678         pos = tmpPos;\r
 679         pos.divScalar(mass);\r
 680 \r
 681         return mq;\r
 682     }\r
 683 \r
 684     /// Recursively walk the tree to do hackwalk calculation\r
 685     HG walkSubTree(double dsq, HG hg) {\r
 686         if (subdivp(dsq, hg)) {\r
 687             for (int k = 0; k < Cell.NSUB; k++) {\r
 688                 Node r = subp[k];\r
 689                 if (r !is null)\r
 690                     hg = r.walkSubTree(dsq / 4.0, hg);\r
 691             }\r
 692         } else {\r
 693             hg = gravSub(hg);\r
 694         }\r
 695         return hg;\r
 696     }\r
 697 \r
 698     /**\r
 699     * Decide if the cell is too close to accept as a single term.\r
 700     * @return true if the cell is too close.\r
 701     **/\r
 702     bool subdivp(double dsq, HG hg) {\r
 703         Vec3 dr = new Vec3();\r
 704         dr.subtraction(pos, hg.pos0);\r
 705         double drsq = dr.dotProduct();\r
 706 \r
 707         // in the original olden version drsp is multiplied by 1.0\r
 708         return drsq < dsq;\r
 709     }\r
 710 \r
 711     /**\r
 712     * Return a string represenation of a cell.\r
 713     * @return a string represenation of a cell.\r
 714     **/\r
 715     public char* toCString() {\r
 716         char* result = cast(char*)malloc(140);\r
 717         if (result == null) exit(1);\r
 718         char* super_str = super.toCString();\r
 719         sprintf(result, "Cell %s", super_str);\r
 720         free(super_str);\r
 721         return result;\r
 722     }\r
 723 }\r
 724 \r
 725 \r
 726 /**\r
 727 * A class that represents the root of the data structure used\r
 728 * to represent the N-bodies in the Barnes-Hut algorithm.\r
 729 **/\r
 730 final class Tree {\r
 731     Vec3 rmin;\r
 732     double rsize;\r
 733 \r
 734     /// A reference to the root node.\r
 735     Node root;\r
 736 \r
 737     /// The complete list of bodies that have been created.\r
 738     private Body bodyTab;\r
 739 \r
 740     /// The complete list of bodies that have been created - in reverse.\r
 741     private Body bodyTabRev;\r
 742 \r
 743     /// Construct the root of the data structure that represents the N-bodies.\r
 744     this() {\r
 745         rmin = new Vec3();\r
 746         rsize = -2.0 * -2.0;\r
 747         root = null;\r
 748         bodyTab = null;\r
 749         bodyTabRev = null;\r
 750         rmin.value(0, -2.0);\r
 751         rmin.value(1, -2.0);\r
 752         rmin.value(2, -2.0);\r
 753     }\r
 754 \r
 755     /**\r
 756     * Return an enumeration of the bodies.\r
 757     * @return an enumeration of the bodies.\r
 758     **/\r
 759     final Enumeration bodies() {\r
 760         return bodyTab.elements();\r
 761     }\r
 762 \r
 763     /**\r
 764     * Return an enumeration of the bodies - in reverse.\r
 765     * @return an enumeration of the bodies - in reverse.\r
 766     **/\r
 767     final Enumeration bodiesRev() {\r
 768         return bodyTabRev.elementsRev();\r
 769     }\r
 770 \r
 771     /**\r
 772     * Create the testdata used in the benchmark.\r
 773     * @param nbody the number of bodies to create\r
 774     **/\r
 775     final void createTestData(int nbody) {\r
 776         Vec3 cmr = new Vec3();\r
 777         Vec3 cmv = new Vec3();\r
 778         Body head = new Body();\r
 779         Body prev = head;\r
 780 \r
 781         double rsc = 3.0 * PI / 16.0;\r
 782         double vsc = sqrt(1.0 / rsc);\r
 783         int seed = 123;\r
 784         Random rnd = Random(seed);\r
 785 \r
 786         for (int i = 0; i < nbody; i++) {\r
 787             Body p = new Body();\r
 788 \r
 789             prev.setNext(p);\r
 790             prev = p;\r
 791             p.mass = 1.0 / cast(double)nbody;\r
 792 \r
 793             double t1 = rnd.uniform(0.0, 0.999);\r
 794             t1 = pow(t1, (-2.0 / 3.0)) - 1.0;\r
 795             double r = 1.0 / sqrt(t1);\r
 796 \r
 797             double coeff = 4.0;\r
 798             for (int k = 0; k < Vec3.NDIM; k++) {\r
 799                 r = rnd.uniform(0.0, 0.999);\r
 800                 p.pos.value(k, coeff * r);\r
 801             }\r
 802 \r
 803             cmr.addition(p.pos);\r
 804 \r
 805             double x, y;\r
 806             do {\r
 807                 x = rnd.uniform(0.0, 1.0);\r
 808                 y = rnd.uniform(0.0, 0.1);\r
 809             } while (y > (x * x * pow(1.0 - x * x, 3.5)));\r
 810 \r
 811             double v = sqrt(2.0) * x / pow(1 + r * r, 0.25);\r
 812 \r
 813             double rad = vsc * v;\r
 814             double rsq;\r
 815             do {\r
 816                 for (int k = 0; k < Vec3.NDIM; k++)\r
 817                     p.vel.value(k, rnd.uniform(-1.0, 1.0));\r
 818                 rsq = p.vel.dotProduct();\r
 819             } while (rsq > 1.0);\r
 820             double rsc1 = rad / sqrt(rsq);\r
 821             p.vel.multScalar(rsc1);\r
 822             cmv.addition(p.vel);\r
 823         }\r
 824 \r
 825         // mark end of list\r
 826         prev.setNext(null);\r
 827 \r
 828         // toss the dummy node at the beginning and set a reference to the first element\r
 829         bodyTab = head.getNext();\r
 830 \r
 831         cmr.divScalar(cast(double)nbody);\r
 832         cmv.divScalar(cast(double)nbody);\r
 833 \r
 834         prev = null;\r
 835 \r
 836         for (Enumeration e = bodyTab.elements(); e.hasMoreElements();) {\r
 837             Body b = cast(Body)e.nextElement();\r
 838             b.pos.subtraction(cmr);\r
 839             b.vel.subtraction(cmv);\r
 840             b.setProcNext(prev);\r
 841             prev = b;\r
 842         }\r
 843 \r
 844         // set the reference to the last element\r
 845         bodyTabRev = prev;\r
 846     }\r
 847 \r
 848 \r
 849     /**\r
 850     * Advance the N-body system one time-step.\r
 851     * @param nstep the current time step\r
 852     **/\r
 853     void stepSystem(int nstep) {\r
 854         // free the tree\r
 855         root = null;\r
 856 \r
 857         makeTree(nstep);\r
 858 \r
 859         // compute the gravity for all the particles\r
 860         for (Enumeration e = bodyTabRev.elementsRev(); e.hasMoreElements();) {\r
 861             Body b = cast(Body)e.nextElement();\r
 862             b.hackGravity(rsize, root);\r
 863         }\r
 864         vp(bodyTabRev, nstep);\r
 865     }\r
 866 \r
 867     /**\r
 868     * Initialize the tree structure for hack force calculation.\r
 869     * @param nsteps the current time step\r
 870     **/\r
 871     private void makeTree(int nstep) {\r
 872         for (Enumeration e = bodiesRev(); e.hasMoreElements();) {\r
 873             Body q = cast(Body)e.nextElement();\r
 874             if (q.mass != 0.0) {\r
 875                 q.expandBox(this, nstep);\r
 876                 Vec3 xqic = intcoord(q.pos);\r
 877                 if (root is null) {\r
 878                     root = q;\r
 879                 } else {\r
 880                     root = root.loadTree(q, xqic, Node.IMAX >> 1, this);\r
 881                 }\r
 882             }\r
 883         }\r
 884         root.hackcofm();\r
 885     }\r
 886 \r
 887     /**\r
 888     * Compute integerized coordinates.\r
 889     * @return the coordinates or null if rp is out of bounds\r
 890     **/\r
 891     final Vec3 intcoord(Vec3 vp) {\r
 892         Vec3 xp = new Vec3();\r
 893 \r
 894         double xsc = (vp.value(0) - rmin.value(0)) / rsize;\r
 895         if (0.0 <= xsc && xsc < 1.0)\r
 896             xp.value(0, floor(Node.IMAX * xsc));\r
 897         else\r
 898             return null;\r
 899 \r
 900         xsc = (vp.value(1) - rmin.value(1)) / rsize;\r
 901         if (0.0 <= xsc && xsc < 1.0)\r
 902             xp.value(1, floor(Node.IMAX * xsc));\r
 903         else\r
 904             return null;\r
 905 \r
 906         xsc = (vp.value(2) - rmin.value(2)) / rsize;\r
 907         if (0.0 <= xsc && xsc < 1.0)\r
 908             xp.value(2, floor(Node.IMAX * xsc));\r
 909         else\r
 910             return null;\r
 911 \r
 912         return xp;\r
 913     }\r
 914 \r
 915     static final private void vp(Body p, int nstep) {\r
 916         Vec3 dacc = new Vec3();\r
 917         Vec3 dvel = new Vec3();\r
 918         double dthf = 0.5 * BH.DTIME;\r
 919 \r
 920         for (Enumeration e = p.elementsRev(); e.hasMoreElements();) {\r
 921             Body b = cast(Body)e.nextElement();\r
 922             Vec3 acc1 = b.newAcc.clone();\r
 923             if (nstep > 0) {\r
 924                 dacc.subtraction(acc1, b.acc);\r
 925                 dvel.multScalar(dacc, dthf);\r
 926                 dvel.addition(b.vel);\r
 927                 b.vel = dvel.clone();\r
 928             }\r
 929 \r
 930             b.acc = acc1.clone();\r
 931             dvel.multScalar(b.acc, dthf);\r
 932 \r
 933             Vec3 vel1 = b.vel.clone();\r
 934             vel1.addition(dvel);\r
 935             Vec3 dpos = vel1.clone();\r
 936             dpos.multScalar(BH.DTIME);\r
 937             dpos.addition(b.pos);\r
 938             b.pos = dpos.clone();\r
 939             vel1.addition(dvel);\r
 940             b.vel = vel1.clone();\r
 941         }\r
 942     }\r
 943 }\r
 944 \r
 945 \r
 946 final class BH {\r
 947     /// The user specified number of bodies to create.\r
 948     private static int nbody = 0;\r
 949 \r
 950     /// The maximum number of time steps to take in the simulation\r
 951     private static int nsteps = 10;\r
 952 \r
 953     /// Should we print information messsages\r
 954     private static bool printMsgs = false;\r
 955 \r
 956     /// Should we print detailed results\r
 957     private static bool printResults = false;\r
 958 \r
 959     const double DTIME = 0.0125;\r
 960     private const double TSTOP = 2.0;\r
 961 \r
 962     public static final void main(char[][] args) {\r
 963         parseCmdLine(args);\r
 964 \r
 965         if (printMsgs)\r
 966             printf("nbody = %d\n", nbody);\r
 967 \r
 968         auto start0 = myclock();\r
 969         Tree root = new Tree();\r
 970         root.createTestData(nbody);\r
 971         auto end0 = myclock();\r
 972         if (printMsgs)\r
 973               printf("Bodies created\n");\r
 974 \r
 975         auto start1 = myclock();\r
 976         double tnow = 0.0;\r
 977         int i = 0;\r
 978         while ((tnow < TSTOP + 0.1 * DTIME) && (i < nsteps)) {\r
 979             root.stepSystem(i++);\r
 980             tnow += DTIME;\r
 981         }\r
 982         auto end1 = myclock();\r
 983 \r
 984         if (printResults) {\r
 985             int j = 0;\r
 986             for (Enumeration e = root.bodies(); e.hasMoreElements();) {\r
 987                 Body b = cast(Body)e.nextElement();\r
 988                 char* str_ptr = b.pos.toCString();\r
 989                 printf("body %d: %s\n", j++, str_ptr);\r
 990                 free(str_ptr);\r
 991             }\r
 992         }\r
 993 \r
 994         if (printMsgs) {\r
 995             printf("Build time %.2f\n", end0 - start0);\r
 996             printf("Compute Time %.2f\n", end1 - start1);\r
 997             printf("Total Time %.2f\n", end1 - start0);\r
 998             printf("Done!\n");\r
 999         }\r
1000     }\r
1001 \r
1002     private static final void parseCmdLine(char[][] args) {\r
1003         int i = 1;\r
1004         while (i < args.length && args[i][0] == '-') {\r
1005             char[] arg = args[i++];\r
1006 \r
1007             // check for options that require arguments\r
1008             if (arg == "-b") {\r
1009                 if (i < args.length)\r
1010                     nbody = toInt(args[i++]);\r
1011                 else\r
1012                     throw new Exception("-l requires the number of levels");\r
1013             } else if (arg == "-s") {\r
1014                 if (i < args.length)\r
1015                     nsteps = toInt(args[i++]);\r
1016                 else\r
1017                     throw new Exception("-l requires the number of levels");\r
1018             } else if (arg == "-m") {\r
1019                 printMsgs = true;\r
1020             } else if (arg == "-p") {\r
1021                 printResults = true;\r
1022             } else if (arg == "-h") {\r
1023                 usage();\r
1024             }\r
1025         }\r
1026 \r
1027         if (nbody == 0)\r
1028             usage();\r
1029     }\r
1030 \r
1031     /// The usage routine which describes the program options.\r
1032     private static final void usage() {\r
1033         fprintf(stderr, "usage: BH -b <size> [-s <steps>] [-p] [-m] [-h]\n");\r
1034         fprintf(stderr, "  -b the number of bodies\n");\r
1035         fprintf(stderr, "  -s the max. number of time steps (default=10)\n");\r
1036         fprintf(stderr, "  -p (print detailed results)\n");\r
1037         fprintf(stderr, "  -m (print information messages\n");\r
1038         exit(0);\r
1039     }\r
1040 }\r
1041 \r
1042 \r
1043 void main(char[][] args) {\r
1044     BH.main(args);\r
1045 }\r