GC := cdgc
-O := build/$(GC)
-DC := dmd
-LD := dmd
+B := build
+O := $B/$(GC)
+DC := ../dmd/src/dmd
+LD := ../dmd/src/dmd
LN := ln
-GNUPLOT := gnuplot
+TIME := /usr/bin/time
override DFLAGS ?= -release -inline -O -gc
LDFLAGS := -defaultlib=tango-$(GC) -debuglib=tango-$(GC)
LIBTANGO := ../lib/libtango-$(GC).a
P = @
P_DC = @printf ' DC %- 40s <- %s\n' '$@' '$(if \
$(subst 1,,$(words $^)),$< ... $(lastword $^),$<)';
-P_PLOT = @printf ' PLOT %- 40s <- %s\n' '$@' '$(filter %.csv,$^)';
+P_PLOT = @printf ' PLOT %- 40s <- %s\n' '$@' '$(notdir $(filter %.csv,$^))';
P_AWK = @printf ' AWK %- 40s <- %s\n' '$@' '$<';
P_RUN = @printf ' RUN $< $(args)\n';
P_MAKE = @printf ' MAKE $@\n';
# create build directories if they don't already exist
dummy_mkdir := $(shell mkdir -p $O $O/bin $O/time $O/stats)
+
+########################################################
+# general rules that doesn't depend on the GC variable #
+########################################################
+
.PHONY: all
-all: cdgc basic
+all: basic cdgc
+
+.PHONY: micro-time
+micro: micro-time
+
+.PHONY: dil
+dil:
+ # TODO
+
+.PHONY: micro-time
+micro-time: $B/time.eps $B/time.svg
+$B/time.%: $(patsubst micro/%.d,$B/time-%.csv,$(wildcard micro/*.d)) \
+ time-plot.tpl-gpi time-plot.sh templite.py
+ $(P_PLOT) ./time-plot.sh $* $@ $(filter %.csv,$^)
+
+.PRECIOUS: $B/time-%.csv
+$B/time-%.csv: $B/basic/time/%.csv | basic cdgc
+ $P echo -n > $@
+ $P for t in basic cdgc; do \
+ (echo -n $$t,; ./stats.py < $B/$$t/time/$*.csv) >> $@; \
+ echo " STATS `tail -n1 $@` >> $@"; \
+ done
.PHONY: cdgc basic
cdgc basic:
- $(P_MAKE) $(MAKE) --no-print-directory micro-time dil-build GC=$@
+ $(P_MAKE) $(MAKE) --no-print-directory micro-gc-build dil-gc-build GC=$@
+.PHONY: clean-all clean-cdgc clean-basic
+clean-all: clean-cdgc clean-basic
+clean-cdgc:
+ $(P_MAKE) $(MAKE) --no-print-directory clean GC=cdgc
+clean-basic:
+ $(P_MAKE) $(MAKE) --no-print-directory clean GC=basic
+
+
+#########################################
+# rules that depends on the GC variable #
+#########################################
+
# micro
########
micro-src := $(wildcard micro/*.d)
-.PHONY: micro-build
-micro-build: $(patsubst micro/%.d,$O/bin/%,$(wildcard micro/*.d))
+.PHONY: micro-gc-build
+micro-gc-build: $(patsubst micro/%.d,$O/bin/%,$(wildcard micro/*.d))
.PRECIOUS: $O/bin/%
-$O/bin/%: $O/micro/%.o
+$O/bin/%: $O/micro/%.o $(LIBTANGO)
$(P_DC) $(DC) $(LDFLAGS) -of$@ $^
-.PHONY: micro-time
-micro-time: $O/time/stats.csv
+.PHONY: micro-gc-time
+micro-gc-time: $(patsubst micro/%.d,$O/time/%.csv,$(wildcard micro/*.d))
-.PHONY: micro-stats
-micro-stats: $(patsubst micro/%.d,$O/stats/%.eps,$(wildcard micro/*.d))
+.PHONY: micro-gc-stats
+micro-gc-stats: $(patsubst micro/%.d,$O/stats/%.eps,$(wildcard micro/*.d))
# special command line arguments for 'shootout_binarytrees' micro benchmark
-$O/time/shootout_binarytrees.t.csv $O/time/shootout_binarytrees.s.csv \
- $O/stats/shootout_binarytrees.c.csv \
+$O/time/shootout_binarytrees.csv $O/stats/shootout_binarytrees.c.csv \
$O/stats/shootout_binarytrees.a.csv: \
override args := 16
# special command line arguments for 'split' micro benchmark
-$O/time/split.t.csv $O/time/split.s.csv \
- $O/stats/split.c.csv $O/stats/split.a.csv: \
- override args := micro/bible.txt
+$O/time/split.csv $O/stats/split.c.csv $O/stats/split.a.csv: \
+ override args := micro/bible.txt 2
# special command line arguments for 'voronoi' micro benchmark
-$O/time/voronoi.t.csv $O/time/voronoi.s.csv \
- $O/stats/voronoi.c.csv $O/voronoi/split.a.csv: \
+$O/time/voronoi.csv $O/stats/voronoi.c.csv $O/voronoi/split.a.csv: \
override args := -n 30000
DIL_SRC = $(wildcard dil/src/*.d dil/src/cmd/*.d dil/src/util/*.d \
dil/src/dil/*.d dil/src/dil/*/*.d)
-.PHONY: dil-nop-stats
-dil-nop-stats: $O/stats/dil-nop.eps
-
-.PHONY: dil-build
-dil-build: $O/bin/dil
+.PHONY: dil-gc-build
+dil-gc-build: $O/bin/dil
$O/bin/dil: override DFLAGS += -Idil/src
$O/bin/dil: $(patsubst %.d,$O/%.o,$(DIL_SRC)) $(LIBTANGO)
$(P_DC) $(DC) $(LDFLAGS) -L-lmpfr -L-lgmp -of$@ $^
-$O/bin/dil-nop: $O/bin/dil
- @$(P_LN) $(LN) -sf $(<F) $@
+.PHONY: dil-gc-time
+dil-gc-time: $O/time/dil.csv
# common rules
$O/%.o: %.d
$(P_DC) $(DC) -c $(DFLAGS) -of$@ $^
-I := 10
-.PRECIOUS: $O/time/%.t.csv
+I := 3
+.PRECIOUS: $O/time/%.csv
ifeq ($F,1)
.PHONY: $(patsubst micro/%.d,$O/bin/%,$(wildcard micro/*.d))
endif
-$O/time/%.t.csv: $O/bin/%
+$O/time/%.csv: $O/bin/%
$P echo -n ' RUN $* $(args) > $@ ($I)'
$P echo -n > $@
$P for i in `seq $I`; do \
echo -n " $$i"; \
- time -f%e -a -o $@ ./$< $(args); \
- done; echo
-
-.PRECIOUS: $O/time/stats.csv
-$O/time/stats.csv: $(patsubst micro/%.d,$O/time/%.t.csv,$(micro-src))
- $P echo -n > $@
- $P for t in $^; do \
- (echo -n `basename $$t`,; ./stats.py < $$t) >> $@; \
- echo " STATS `tail -n1 $@` >> $@"; \
- done
+ $(TIME) -f%e -a -o $@ ./$< $(args); \
+ done; echo
.PRECIOUS: $O/stats/%.c.csv $O/stats/%.a.csv
$O/stats/%.c.csv $O/stats/%.a.csv: $O/bin/%
clean: $O/
$(P_RM) $(RM) -r $^
-.PHONY: clean-all
-clean-all: clean-cdgc clean-basic
-clean-cdgc:
- $(P_MAKE) $(MAKE) --no-print-directory clean GC=cdgc
-clean-basic:
- $(P_MAKE) $(MAKE) --no-print-directory clean GC=basic
-
--- /dev/null
+/**\r
+A D implementation of the _bh_ Olden benchmark.\r
+The Olden benchmark implements the Barnes-Hut benchmark\r
+that is decribed in:\r
+\r
+J. Barnes and P. Hut, "A hierarchical o(N log N) force-calculation algorithm",\r
+Nature, 324:446-449, Dec. 1986\r
+\r
+The original code in the Olden benchmark suite is derived from the\r
+ftp://hubble.ifa.hawaii.edu/pub/barnes/treecode\r
+source distributed by Barnes.\r
+\r
+Java code converted to D by leonardo maffi, V.1.0, Dec 25 2009.\r
+\r
+Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.\r
+Downloaded from http://www.fantascienza.net/leonardo/js/\r
+ http://www.fantascienza.net/leonardo/js/dolden_bh.zip\r
+*/\r
+\r
+\r
+version (Tango) {\r
+ import tango.stdc.stdio: printf, sprintf, fprintf, stderr;\r
+ import tango.stdc.stdlib: exit, malloc, free;\r
+ import tango.stdc.time: CLOCKS_PER_SEC, clock;\r
+ import tango.math.Math: sqrt, floor, PI, pow;\r
+\r
+ import Integer = tango.text.convert.Integer;\r
+ alias Integer.parse toInt;\r
+} else {\r
+ import std.c.stdio: printf, sprintf, fprintf, stderr;\r
+ import std.c.stdlib: exit, malloc, free;\r
+ import std.c.time: CLOCKS_PER_SEC, clock;\r
+ import std.math: sqrt, floor, PI, pow;\r
+\r
+ version (D_Version2) {\r
+ import std.conv: to;\r
+ alias to!(int, char[]) toInt;\r
+ } else {\r
+ import std.conv: toInt;\r
+ }\r
+}\r
+\r
+\r
+double myclock() {\r
+ return clock() / cast(double)CLOCKS_PER_SEC;\r
+}\r
+\r
+\r
+/*\r
+Basic uniform random generator: Minimal Standard in Park and\r
+Miller (1988): "Random Number Generators: Good Ones Are Hard to\r
+Find", Comm. of the ACM, 31, 1192-1201.\r
+Parameters: m = 2^31-1, a=48271.\r
+\r
+Adapted from Pascal code by Jesper Lund:\r
+http://www.gnu-pascal.de/crystal/gpc/en/mail1390.html\r
+*/\r
+struct Random {\r
+ const int m = int.max;\r
+ const int a = 48_271;\r
+ const int q = m / a;\r
+ const int r = m % a;\r
+ int seed;\r
+\r
+ public double uniform(double min, double max) {\r
+ int k = seed / q;\r
+ seed = a * (seed - k * q) - r * k;\r
+ if (seed < 1)\r
+ seed += m;\r
+ double r = cast(double)seed / m;\r
+ return r * (max - min) + min;\r
+ }\r
+}\r
+\r
+\r
+interface Enumeration {\r
+ bool hasMoreElements();\r
+ Object nextElement();\r
+}\r
+\r
+\r
+/**\r
+A class representing a three dimensional vector that implements\r
+several math operations. To improve speed we implement the\r
+vector as an array of doubles rather than use the exising\r
+code in the java.util.Vec3 class.\r
+*/\r
+final class Vec3 {\r
+ /// The number of dimensions in the vector\r
+ // Can't be changed because of crossProduct()\r
+ const int NDIM = 3;\r
+\r
+ /// An array containing the values in the vector.\r
+ private double data[];\r
+\r
+ /// Construct an empty 3 dimensional vector for use in Barnes-Hut algorithm.\r
+ this() {\r
+ data.length = NDIM;\r
+ data[] = 0.0;\r
+ }\r
+\r
+ /// Create a copy of the vector. Return a clone of the math vector\r
+ public Vec3 clone() {\r
+ Vec3 v = new Vec3;\r
+ v.data.length = NDIM;\r
+ v.data[] = this.data[];\r
+ return v;\r
+ }\r
+\r
+ /**\r
+ Return the value at the i'th index of the vector.\r
+ @param i the vector index\r
+ @return the value at the i'th index of the vector.\r
+ */\r
+ final double value(int i) {\r
+ return data[i];\r
+ }\r
+\r
+ /**\r
+ Set the value of the i'th index of the vector.\r
+ @param i the vector index\r
+ @param v the value to store\r
+ */\r
+ final void value(int i, double v) {\r
+ data[i] = v;\r
+ }\r
+\r
+ /**\r
+ Set one of the dimensions of the vector to 1.0\r
+ param j the dimension to set.\r
+ */\r
+ final void unit(int j) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = (i == j) ? 1.0 : 0.0;\r
+ }\r
+\r
+ /**\r
+ Add two vectors and the result is placed in this vector.\r
+ @param u the other operand of the addition\r
+ */\r
+ final void addition(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] += u.data[i];\r
+ }\r
+\r
+ /**\r
+ * Subtract two vectors and the result is placed in this vector.\r
+ * This vector contain the first operand.\r
+ * @param u the other operand of the subtraction.\r
+ **/\r
+ final void subtraction(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] -= u.data[i];\r
+ }\r
+\r
+ /**\r
+ Subtract two vectors and the result is placed in this vector.\r
+ @param u the first operand of the subtraction.\r
+ @param v the second opernd of the subtraction\r
+ */\r
+ final void subtraction(Vec3 u, Vec3 v) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = u.data[i] - v.data[i];\r
+ }\r
+\r
+ /**\r
+ Multiply the vector times a scalar.\r
+ @param s the scalar value\r
+ **/\r
+ final void multScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] *= s;\r
+ }\r
+\r
+ /**\r
+ * Multiply the vector times a scalar and place the result in this vector.\r
+ * @param u the vector\r
+ * @param s the scalar value\r
+ **/\r
+ final void multScalar(Vec3 u, double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = u.data[i] * s;\r
+ }\r
+\r
+ /**\r
+ * Divide each element of the vector by a scalar value.\r
+ * @param s the scalar value.\r
+ **/\r
+ final void divScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] /= s;\r
+ }\r
+\r
+ /**\r
+ * Return the dot product of a vector.\r
+ * @return the dot product of a vector.\r
+ **/\r
+ final double dotProduct() {\r
+ double s = 0.0;\r
+ for (int i = 0; i < NDIM; i++)\r
+ s += data[i] * data[i];\r
+ return s;\r
+ }\r
+\r
+ final double absolute() {\r
+ double tmp = 0.0;\r
+ for (int i = 0; i < NDIM; i++)\r
+ tmp += data[i] * data[i];\r
+ return sqrt(tmp);\r
+ }\r
+\r
+ final double distance(Vec3 v) {\r
+ double tmp = 0.0;\r
+ for (int i = 0; i < NDIM; i++)\r
+ tmp += ((data[i] - v.data[i]) * (data[i] - v.data[i]));\r
+ return sqrt(tmp);\r
+ }\r
+\r
+ final void crossProduct(Vec3 u, Vec3 w) {\r
+ data[0] = u.data[1] * w.data[2] - u.data[2] * w.data[1];\r
+ data[1] = u.data[2] * w.data[0] - u.data[0] * w.data[2];\r
+ data[2] = u.data[0] * w.data[1] - u.data[1] * w.data[0];\r
+ }\r
+\r
+ final void incrementalAdd(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] += u.data[i];\r
+ }\r
+\r
+ final void incrementalSub(Vec3 u) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] -= u.data[i];\r
+ }\r
+\r
+ final void incrementalMultScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] *= s;\r
+ }\r
+\r
+ final void incrementalDivScalar(double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] /= s;\r
+ }\r
+\r
+ final void addScalar(Vec3 u, double s) {\r
+ for (int i = 0; i < NDIM; i++)\r
+ data[i] = u.data[i] + s;\r
+ }\r
+\r
+ public char* toCString() {\r
+ // this is not generic code at all\r
+ char* result = cast(char*)malloc(100);\r
+ if (result == null) exit(1);\r
+ sprintf(result, "%.17f %.17f %.17f ", data[0], data[1], data[2]);\r
+ return result;\r
+ }\r
+}\r
+\r
+\r
+/// A class that represents the common fields of a cell or body data structure.\r
+abstract class Node {\r
+ // mass of the node\r
+ double mass;\r
+\r
+ // Position of the node\r
+ Vec3 pos;\r
+\r
+ // highest bit of int coord\r
+ const int IMAX = 1073741824;\r
+\r
+ // potential softening parameter\r
+ const double EPS = 0.05;\r
+\r
+ /// Construct an empty node\r
+ protected this() {\r
+ mass = 0.0;\r
+ pos = new Vec3();\r
+ }\r
+\r
+ abstract Node loadTree(Body p, Vec3 xpic, int l, Tree root);\r
+\r
+ abstract double hackcofm();\r
+\r
+ abstract HG walkSubTree(double dsq, HG hg);\r
+\r
+ static final int oldSubindex(Vec3 ic, int l) {\r
+ int i = 0;\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ if ((cast(int)ic.value(k) & l) != 0)\r
+ i += Cell.NSUB >> (k + 1);\r
+ return i;\r
+ }\r
+\r
+ public char* toCString() {\r
+ char* result = cast(char*)malloc(130);\r
+ if (result == null) exit(1);\r
+ char* pos_str = pos.toCString();\r
+ sprintf(result, "%f : %s", mass, pos_str);\r
+ free(pos_str);\r
+ return result;\r
+ }\r
+\r
+ /// Compute a single body-body or body-cell interaction\r
+ final HG gravSub(HG hg) {\r
+ Vec3 dr = new Vec3();\r
+ dr.subtraction(pos, hg.pos0);\r
+\r
+ double drsq = dr.dotProduct() + (EPS * EPS);\r
+ double drabs = sqrt(drsq);\r
+\r
+ double phii = mass / drabs;\r
+ hg.phi0 -= phii;\r
+ double mor3 = phii / drsq;\r
+ dr.multScalar(mor3);\r
+ hg.acc0.addition(dr);\r
+ return hg;\r
+ }\r
+\r
+ /**\r
+ * A sub class which is used to compute and save information during the\r
+ * gravity computation phase.\r
+ **/\r
+ static protected final class HG {\r
+ /// Body to skip in force evaluation\r
+ Body pskip;\r
+\r
+ /// Point at which to evaluate field\r
+ Vec3 pos0;\r
+\r
+ /// Computed potential at pos0\r
+\r
+ double phi0;\r
+\r
+ /// computed acceleration at pos0\r
+ Vec3 acc0;\r
+\r
+ /**\r
+ * Create a HG object.\r
+ * @param b the body object\r
+ * @param p a vector that represents the body\r
+ **/\r
+ this(Body b, Vec3 p) {\r
+ pskip = b;\r
+ pos0 = p.clone();\r
+ phi0 = 0.0;\r
+ acc0 = new Vec3();\r
+ }\r
+ }\r
+}\r
+\r
+\r
+/// A class used to representing particles in the N-body simulation.\r
+final class Body : Node {\r
+ Vec3 vel, acc, newAcc;\r
+ double phi = 0.0;\r
+ private Body next, procNext;\r
+\r
+ /// Create an empty body.\r
+ this() {\r
+ vel = new Vec3();\r
+ acc = new Vec3();\r
+ newAcc = new Vec3();\r
+ }\r
+\r
+ /**\r
+ * Set the next body in the list.\r
+ * @param n the body\r
+ **/\r
+ final void setNext(Body n) {\r
+ next = n;\r
+ }\r
+\r
+ /**\r
+ * Get the next body in the list.\r
+ * @return the next body\r
+ **/\r
+ final Body getNext() {\r
+ return next;\r
+ }\r
+\r
+ /**\r
+ * Set the next body in the list.\r
+ * @param n the body\r
+ **/\r
+ final void setProcNext(Body n) {\r
+ procNext = n;\r
+ }\r
+\r
+ /**\r
+ * Get the next body in the list.\r
+ * @return the next body\r
+ **/\r
+ final Body getProcNext() {\r
+ return procNext;\r
+ }\r
+\r
+ /**\r
+ * Enlarge cubical "box", salvaging existing tree structure.\r
+ * @param tree the root of the tree.\r
+ * @param nsteps the current time step\r
+ **/\r
+ final void expandBox(Tree tree, int nsteps) {\r
+ Vec3 rmid = new Vec3();\r
+\r
+ bool inbox = icTest(tree);\r
+ while (!inbox) {\r
+ double rsize = tree.rsize;\r
+ rmid.addScalar(tree.rmin, 0.5 * rsize);\r
+\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ if (pos.value(k) < rmid.value(k)) {\r
+ double rmin = tree.rmin.value(k);\r
+ tree.rmin.value(k, rmin - rsize);\r
+ }\r
+\r
+ tree.rsize = 2.0 * rsize;\r
+ if (tree.root !is null) {\r
+ Vec3 ic = tree.intcoord(rmid);\r
+ if (ic is null)\r
+ throw new Exception("Value is out of bounds");\r
+ int k = oldSubindex(ic, IMAX >> 1);\r
+ Cell newt = new Cell();\r
+ newt.subp[k] = tree.root;\r
+ tree.root = newt;\r
+ inbox = icTest(tree);\r
+ }\r
+ }\r
+ }\r
+\r
+ /**\r
+ * Check the bounds of the body and return true if it isn't in the\r
+ * correct bounds.\r
+ **/\r
+ final bool icTest(Tree tree) {\r
+ double pos0 = pos.value(0);\r
+ double pos1 = pos.value(1);\r
+ double pos2 = pos.value(2);\r
+\r
+ // by default, it is in bounds\r
+ bool result = true;\r
+\r
+ double xsc = (pos0 - tree.rmin.value(0)) / tree.rsize;\r
+ if (!(0.0 < xsc && xsc < 1.0))\r
+ result = false;\r
+\r
+ xsc = (pos1 - tree.rmin.value(1)) / tree.rsize;\r
+ if (!(0.0 < xsc && xsc < 1.0))\r
+ result = false;\r
+\r
+ xsc = (pos2 - tree.rmin.value(2)) / tree.rsize;\r
+ if (!(0.0 < xsc && xsc < 1.0))\r
+ result = false;\r
+\r
+ return result;\r
+ }\r
+\r
+ /**\r
+ * Descend Tree and insert particle. We're at a body so we need to\r
+ * create a cell and attach this body to the cell.\r
+ * @param p the body to insert\r
+ * @param xpic\r
+ * @param l\r
+ * @param tree the root of the data structure\r
+ * @return the subtree with the new body inserted\r
+ **/\r
+ final Node loadTree(Body p, Vec3 xpic, int l, Tree tree) {\r
+ // create a Cell\r
+ Cell retval = new Cell();\r
+ int si = subindex(tree, l);\r
+ // attach this Body node to the cell\r
+ retval.subp[si] = this;\r
+\r
+ // move down one level\r
+ si = oldSubindex(xpic, l);\r
+ Node rt = retval.subp[si];\r
+ if (rt !is null)\r
+ retval.subp[si] = rt.loadTree(p, xpic, l >> 1, tree);\r
+ else\r
+ retval.subp[si] = p;\r
+ return retval;\r
+ }\r
+\r
+ /**\r
+ * Descend tree finding center of mass coordinates\r
+ * @return the mass of this node\r
+ **/\r
+ final double hackcofm() {\r
+ return mass;\r
+ }\r
+\r
+ /// iteration of the bodies\r
+ int opApply(int delegate(ref Body) dg) {\r
+ int result;\r
+ Body current = this;\r
+ while (current !is null) {\r
+ result = dg(current);\r
+ current = current.next;\r
+ if (result)\r
+ break;\r
+ }\r
+ return result;\r
+ }\r
+\r
+ /// inverse iteration of the bodies\r
+ int opApplyReverse(int delegate(ref Body) dg) {\r
+ int result;\r
+ Body current = this;\r
+ while (current !is null) {\r
+ result = dg(current);\r
+ current = current.procNext;\r
+ if (result)\r
+ break;\r
+ }\r
+ return result;\r
+ }\r
+\r
+\r
+ /**\r
+ * Return an enumeration of the bodies\r
+ * @return an enumeration of the bodies\r
+ **/\r
+ final Enumeration elements() {\r
+ // a local class that implements the enumerator\r
+ static final class Enumerate : Enumeration {\r
+ private Body current;\r
+ public this(Body outer) {\r
+ //this.current = Body.this;\r
+ this.current = outer;\r
+ }\r
+ public bool hasMoreElements() {\r
+ return current !is null;\r
+ }\r
+ public Object nextElement() {\r
+ Object retval = current;\r
+ current = current.next;\r
+ return retval;\r
+ }\r
+ }\r
+\r
+ return new Enumerate(this);\r
+ }\r
+\r
+ final Enumeration elementsRev() {\r
+ // a local class that implements the enumerator\r
+ static class Enumerate : Enumeration {\r
+ private Body current;\r
+ public this(Body outer) {\r
+ //this.current = Body.this;\r
+ this.current = outer;\r
+ }\r
+ public bool hasMoreElements() {\r
+ return current !is null;\r
+ }\r
+ public Object nextElement() {\r
+ Object retval = current;\r
+ current = current.procNext;\r
+ return retval;\r
+ }\r
+ }\r
+\r
+ return new Enumerate(this);\r
+ }\r
+\r
+ /**\r
+ * Determine which subcell to select.\r
+ * Combination of intcoord and oldSubindex.\r
+ * @param t the root of the tree\r
+ **/\r
+ final int subindex(Tree tree, int l) {\r
+ Vec3 xp = new Vec3();\r
+\r
+ double xsc = (pos.value(0) - tree.rmin.value(0)) / tree.rsize;\r
+ xp.value(0, floor(IMAX * xsc));\r
+\r
+ xsc = (pos.value(1) - tree.rmin.value(1)) / tree.rsize;\r
+ xp.value(1, floor(IMAX * xsc));\r
+\r
+ xsc = (pos.value(2) - tree.rmin.value(2)) / tree.rsize;\r
+ xp.value(2, floor(IMAX * xsc));\r
+\r
+ int i = 0;\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ if ((cast(int)xp.value(k) & l) != 0)\r
+ i += Cell.NSUB >> (k + 1);\r
+ return i;\r
+ }\r
+\r
+ /**\r
+ * Evaluate gravitational field on the body.\r
+ * The original olden version calls a routine named "walkscan",\r
+ * but we use the same name that is in the Barnes code.\r
+ **/\r
+ final void hackGravity(double rsize, Node root) {\r
+ Vec3 pos0 = pos.clone();\r
+ HG hg = new HG(this, pos);\r
+ hg = root.walkSubTree(rsize * rsize, hg);\r
+ phi = hg.phi0;\r
+ newAcc = hg.acc0;\r
+ }\r
+\r
+ /// Recursively walk the tree to do hackwalk calculation\r
+ final HG walkSubTree(double dsq, HG hg) {\r
+ if (this != hg.pskip)\r
+ hg = gravSub(hg);\r
+ return hg;\r
+ }\r
+\r
+ /**\r
+ * Return a string represenation of a body.\r
+ * @return a string represenation of a body.\r
+ **/\r
+ public char* toCString() {\r
+ char* result = cast(char*)malloc(140);\r
+ if (result == null) exit(1);\r
+ char* super_str = super.toCString();\r
+ sprintf(result, "Body %s", super_str);\r
+ free(super_str);\r
+ return result;\r
+ }\r
+}\r
+\r
+\r
+\r
+/// A class used to represent internal nodes in the tree\r
+final class Cell : Node {\r
+ // subcells per cell\r
+ const int NSUB = 1 << Vec3.NDIM;\r
+\r
+ /**\r
+ * The children of this cell node. Each entry may contain either\r
+ * another cell or a body.\r
+ **/\r
+ Node[] subp;\r
+ Cell next;\r
+\r
+ this() {\r
+ subp.length = NSUB;\r
+ }\r
+\r
+ /**\r
+ * Descend Tree and insert particle. We're at a cell so\r
+ * we need to move down the tree.\r
+ * @param p the body to insert into the tree\r
+ * @param xpic\r
+ * @param l\r
+ * @param tree the root of the tree\r
+ * @return the subtree with the new body inserted\r
+ **/\r
+ Node loadTree(Body p, Vec3 xpic, int l, Tree tree) {\r
+ // move down one level\r
+ int si = oldSubindex(xpic, l);\r
+ Node rt = subp[si];\r
+ if (rt !is null)\r
+ subp[si] = rt.loadTree(p, xpic, l >> 1, tree);\r
+ else\r
+ subp[si] = p;\r
+ return this;\r
+ }\r
+\r
+ /**\r
+ * Descend tree finding center of mass coordinates\r
+ * @return the mass of this node\r
+ **/\r
+ double hackcofm() {\r
+ double mq = 0.0;\r
+ Vec3 tmpPos = new Vec3();\r
+ Vec3 tmpv = new Vec3();\r
+ for (int i = 0; i < NSUB; i++) {\r
+ Node r = subp[i];\r
+ if (r !is null) {\r
+ double mr = r.hackcofm();\r
+ mq = mr + mq;\r
+ tmpv.multScalar(r.pos, mr);\r
+ tmpPos.addition(tmpv);\r
+ }\r
+ }\r
+ mass = mq;\r
+ pos = tmpPos;\r
+ pos.divScalar(mass);\r
+\r
+ return mq;\r
+ }\r
+\r
+ /// Recursively walk the tree to do hackwalk calculation\r
+ HG walkSubTree(double dsq, HG hg) {\r
+ if (subdivp(dsq, hg)) {\r
+ for (int k = 0; k < Cell.NSUB; k++) {\r
+ Node r = subp[k];\r
+ if (r !is null)\r
+ hg = r.walkSubTree(dsq / 4.0, hg);\r
+ }\r
+ } else {\r
+ hg = gravSub(hg);\r
+ }\r
+ return hg;\r
+ }\r
+\r
+ /**\r
+ * Decide if the cell is too close to accept as a single term.\r
+ * @return true if the cell is too close.\r
+ **/\r
+ bool subdivp(double dsq, HG hg) {\r
+ Vec3 dr = new Vec3();\r
+ dr.subtraction(pos, hg.pos0);\r
+ double drsq = dr.dotProduct();\r
+\r
+ // in the original olden version drsp is multiplied by 1.0\r
+ return drsq < dsq;\r
+ }\r
+\r
+ /**\r
+ * Return a string represenation of a cell.\r
+ * @return a string represenation of a cell.\r
+ **/\r
+ public char* toCString() {\r
+ char* result = cast(char*)malloc(140);\r
+ if (result == null) exit(1);\r
+ char* super_str = super.toCString();\r
+ sprintf(result, "Cell %s", super_str);\r
+ free(super_str);\r
+ return result;\r
+ }\r
+}\r
+\r
+\r
+/**\r
+* A class that represents the root of the data structure used\r
+* to represent the N-bodies in the Barnes-Hut algorithm.\r
+**/\r
+final class Tree {\r
+ Vec3 rmin;\r
+ double rsize;\r
+\r
+ /// A reference to the root node.\r
+ Node root;\r
+\r
+ /// The complete list of bodies that have been created.\r
+ private Body bodyTab;\r
+\r
+ /// The complete list of bodies that have been created - in reverse.\r
+ private Body bodyTabRev;\r
+\r
+ /// Construct the root of the data structure that represents the N-bodies.\r
+ this() {\r
+ rmin = new Vec3();\r
+ rsize = -2.0 * -2.0;\r
+ root = null;\r
+ bodyTab = null;\r
+ bodyTabRev = null;\r
+ rmin.value(0, -2.0);\r
+ rmin.value(1, -2.0);\r
+ rmin.value(2, -2.0);\r
+ }\r
+\r
+ /**\r
+ * Return an enumeration of the bodies.\r
+ * @return an enumeration of the bodies.\r
+ **/\r
+ final Enumeration bodies() {\r
+ return bodyTab.elements();\r
+ }\r
+\r
+ /**\r
+ * Return an enumeration of the bodies - in reverse.\r
+ * @return an enumeration of the bodies - in reverse.\r
+ **/\r
+ final Enumeration bodiesRev() {\r
+ return bodyTabRev.elementsRev();\r
+ }\r
+\r
+ /**\r
+ * Create the testdata used in the benchmark.\r
+ * @param nbody the number of bodies to create\r
+ **/\r
+ final void createTestData(int nbody) {\r
+ Vec3 cmr = new Vec3();\r
+ Vec3 cmv = new Vec3();\r
+ Body head = new Body();\r
+ Body prev = head;\r
+\r
+ double rsc = 3.0 * PI / 16.0;\r
+ double vsc = sqrt(1.0 / rsc);\r
+ int seed = 123;\r
+ Random rnd = Random(seed);\r
+\r
+ for (int i = 0; i < nbody; i++) {\r
+ Body p = new Body();\r
+\r
+ prev.setNext(p);\r
+ prev = p;\r
+ p.mass = 1.0 / cast(double)nbody;\r
+\r
+ double t1 = rnd.uniform(0.0, 0.999);\r
+ t1 = pow(t1, (-2.0 / 3.0)) - 1.0;\r
+ double r = 1.0 / sqrt(t1);\r
+\r
+ double coeff = 4.0;\r
+ for (int k = 0; k < Vec3.NDIM; k++) {\r
+ r = rnd.uniform(0.0, 0.999);\r
+ p.pos.value(k, coeff * r);\r
+ }\r
+\r
+ cmr.addition(p.pos);\r
+\r
+ double x, y;\r
+ do {\r
+ x = rnd.uniform(0.0, 1.0);\r
+ y = rnd.uniform(0.0, 0.1);\r
+ } while (y > (x * x * pow(1.0 - x * x, 3.5)));\r
+\r
+ double v = sqrt(2.0) * x / pow(1 + r * r, 0.25);\r
+\r
+ double rad = vsc * v;\r
+ double rsq;\r
+ do {\r
+ for (int k = 0; k < Vec3.NDIM; k++)\r
+ p.vel.value(k, rnd.uniform(-1.0, 1.0));\r
+ rsq = p.vel.dotProduct();\r
+ } while (rsq > 1.0);\r
+ double rsc1 = rad / sqrt(rsq);\r
+ p.vel.multScalar(rsc1);\r
+ cmv.addition(p.vel);\r
+ }\r
+\r
+ // mark end of list\r
+ prev.setNext(null);\r
+\r
+ // toss the dummy node at the beginning and set a reference to the first element\r
+ bodyTab = head.getNext();\r
+\r
+ cmr.divScalar(cast(double)nbody);\r
+ cmv.divScalar(cast(double)nbody);\r
+\r
+ prev = null;\r
+\r
+ for (Enumeration e = bodyTab.elements(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ b.pos.subtraction(cmr);\r
+ b.vel.subtraction(cmv);\r
+ b.setProcNext(prev);\r
+ prev = b;\r
+ }\r
+\r
+ // set the reference to the last element\r
+ bodyTabRev = prev;\r
+ }\r
+\r
+\r
+ /**\r
+ * Advance the N-body system one time-step.\r
+ * @param nstep the current time step\r
+ **/\r
+ void stepSystem(int nstep) {\r
+ // free the tree\r
+ root = null;\r
+\r
+ makeTree(nstep);\r
+\r
+ // compute the gravity for all the particles\r
+ for (Enumeration e = bodyTabRev.elementsRev(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ b.hackGravity(rsize, root);\r
+ }\r
+ vp(bodyTabRev, nstep);\r
+ }\r
+\r
+ /**\r
+ * Initialize the tree structure for hack force calculation.\r
+ * @param nsteps the current time step\r
+ **/\r
+ private void makeTree(int nstep) {\r
+ for (Enumeration e = bodiesRev(); e.hasMoreElements();) {\r
+ Body q = cast(Body)e.nextElement();\r
+ if (q.mass != 0.0) {\r
+ q.expandBox(this, nstep);\r
+ Vec3 xqic = intcoord(q.pos);\r
+ if (root is null) {\r
+ root = q;\r
+ } else {\r
+ root = root.loadTree(q, xqic, Node.IMAX >> 1, this);\r
+ }\r
+ }\r
+ }\r
+ root.hackcofm();\r
+ }\r
+\r
+ /**\r
+ * Compute integerized coordinates.\r
+ * @return the coordinates or null if rp is out of bounds\r
+ **/\r
+ final Vec3 intcoord(Vec3 vp) {\r
+ Vec3 xp = new Vec3();\r
+\r
+ double xsc = (vp.value(0) - rmin.value(0)) / rsize;\r
+ if (0.0 <= xsc && xsc < 1.0)\r
+ xp.value(0, floor(Node.IMAX * xsc));\r
+ else\r
+ return null;\r
+\r
+ xsc = (vp.value(1) - rmin.value(1)) / rsize;\r
+ if (0.0 <= xsc && xsc < 1.0)\r
+ xp.value(1, floor(Node.IMAX * xsc));\r
+ else\r
+ return null;\r
+\r
+ xsc = (vp.value(2) - rmin.value(2)) / rsize;\r
+ if (0.0 <= xsc && xsc < 1.0)\r
+ xp.value(2, floor(Node.IMAX * xsc));\r
+ else\r
+ return null;\r
+\r
+ return xp;\r
+ }\r
+\r
+ static final private void vp(Body p, int nstep) {\r
+ Vec3 dacc = new Vec3();\r
+ Vec3 dvel = new Vec3();\r
+ double dthf = 0.5 * BH.DTIME;\r
+\r
+ for (Enumeration e = p.elementsRev(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ Vec3 acc1 = b.newAcc.clone();\r
+ if (nstep > 0) {\r
+ dacc.subtraction(acc1, b.acc);\r
+ dvel.multScalar(dacc, dthf);\r
+ dvel.addition(b.vel);\r
+ b.vel = dvel.clone();\r
+ }\r
+\r
+ b.acc = acc1.clone();\r
+ dvel.multScalar(b.acc, dthf);\r
+\r
+ Vec3 vel1 = b.vel.clone();\r
+ vel1.addition(dvel);\r
+ Vec3 dpos = vel1.clone();\r
+ dpos.multScalar(BH.DTIME);\r
+ dpos.addition(b.pos);\r
+ b.pos = dpos.clone();\r
+ vel1.addition(dvel);\r
+ b.vel = vel1.clone();\r
+ }\r
+ }\r
+}\r
+\r
+\r
+final class BH {\r
+ /// The user specified number of bodies to create.\r
+ private static int nbody = 0;\r
+\r
+ /// The maximum number of time steps to take in the simulation\r
+ private static int nsteps = 10;\r
+\r
+ /// Should we print information messsages\r
+ private static bool printMsgs = false;\r
+\r
+ /// Should we print detailed results\r
+ private static bool printResults = false;\r
+\r
+ const double DTIME = 0.0125;\r
+ private const double TSTOP = 2.0;\r
+\r
+ public static final void main(char[][] args) {\r
+ parseCmdLine(args);\r
+\r
+ if (printMsgs)\r
+ printf("nbody = %d\n", nbody);\r
+\r
+ auto start0 = myclock();\r
+ Tree root = new Tree();\r
+ root.createTestData(nbody);\r
+ auto end0 = myclock();\r
+ if (printMsgs)\r
+ printf("Bodies created\n");\r
+\r
+ auto start1 = myclock();\r
+ double tnow = 0.0;\r
+ int i = 0;\r
+ while ((tnow < TSTOP + 0.1 * DTIME) && (i < nsteps)) {\r
+ root.stepSystem(i++);\r
+ tnow += DTIME;\r
+ }\r
+ auto end1 = myclock();\r
+\r
+ if (printResults) {\r
+ int j = 0;\r
+ for (Enumeration e = root.bodies(); e.hasMoreElements();) {\r
+ Body b = cast(Body)e.nextElement();\r
+ char* str_ptr = b.pos.toCString();\r
+ printf("body %d: %s\n", j++, str_ptr);\r
+ free(str_ptr);\r
+ }\r
+ }\r
+\r
+ if (printMsgs) {\r
+ printf("Build time %.2f\n", end0 - start0);\r
+ printf("Compute Time %.2f\n", end1 - start1);\r
+ printf("Total Time %.2f\n", end1 - start0);\r
+ printf("Done!\n");\r
+ }\r
+ }\r
+\r
+ private static final void parseCmdLine(char[][] args) {\r
+ int i = 1;\r
+ while (i < args.length && args[i][0] == '-') {\r
+ char[] arg = args[i++];\r
+\r
+ // check for options that require arguments\r
+ if (arg == "-b") {\r
+ if (i < args.length)\r
+ nbody = toInt(args[i++]);\r
+ else\r
+ throw new Exception("-l requires the number of levels");\r
+ } else if (arg == "-s") {\r
+ if (i < args.length)\r
+ nsteps = toInt(args[i++]);\r
+ else\r
+ throw new Exception("-l requires the number of levels");\r
+ } else if (arg == "-m") {\r
+ printMsgs = true;\r
+ } else if (arg == "-p") {\r
+ printResults = true;\r
+ } else if (arg == "-h") {\r
+ usage();\r
+ }\r
+ }\r
+\r
+ if (nbody == 0)\r
+ usage();\r
+ }\r
+\r
+ /// The usage routine which describes the program options.\r
+ private static final void usage() {\r
+ fprintf(stderr, "usage: BH -b <size> [-s <steps>] [-p] [-m] [-h]\n");\r
+ fprintf(stderr, " -b the number of bodies\n");\r
+ fprintf(stderr, " -s the max. number of time steps (default=10)\n");\r
+ fprintf(stderr, " -p (print detailed results)\n");\r
+ fprintf(stderr, " -m (print information messages\n");\r
+ exit(0);\r
+ }\r
+}\r
+\r
+\r
+void main(char[][] args) {\r
+ BH.main(args);\r
+}\r
--- /dev/null
+/*
+A Java implementation of the bisort Olden benchmark.
+The Olden benchmark implements a Bitonic Sort as described in:
+
+G. Bilardi and A. Nicolau, "Adaptive Bitonic Sorting: An optimal parallel
+algorithm for shared-memory machines." SIAM J. Comput. 18(2):216-228, 1998.
+
+The benchmarks sorts N numbers where N is a power of 2. If the user provides
+an input value that is not a power of 2, then we use the nearest power of
+2 value that is less than the input value.
+
+Converted to D and obtimized by leonardo maffi, V.1.0, Oct 31 2009.
+
+Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.
+Downloaded from http://www.fantascienza.net/leonardo/js/
+ http://www.fantascienza.net/leonardo/js/dolden_bisort.zip
+*/
+
+version (Tango) {
+ import tango.stdc.stdio: printf, fprintf, stderr;
+ import tango.stdc.stdlib: exit;
+ import Integer = tango.text.convert.Integer;
+ alias Integer.parse toInt;
+ import tango.stdc.time: CLOCKS_PER_SEC, clock;
+} else {
+ import std.c.stdio: printf, fprintf, stderr;
+ import std.c.stdlib: exit;
+ import std.conv: toInt;
+ import std.c.time: CLOCKS_PER_SEC, clock;
+}
+
+
+double myclock() {
+ return clock() / cast(float)CLOCKS_PER_SEC;
+}
+
+
+/**
+ * A class that represents a value to be sorted by the <tt>BiSort</tt>
+ * algorithm. We represents a values as a node in a binary tree.
+ **/
+final class Value {
+ private int value;
+ private Value left, right;
+
+ const bool FORWARD = false;
+ const bool BACKWARD = true;
+
+ // These are used by the Olden benchmark random no. generator
+ private const int CONST_m1 = 10000;
+ private const int CONST_b = 31415821;
+ const int RANGE = 100;
+
+ /**
+ * Constructor for a node representing a value in the bitonic sort tree.
+ * @param v the integer value which is the sort key
+ **/
+ this(int v) {
+ value = v;
+ left = right = null;
+ }
+
+
+ /**
+ * Create a random tree of value to be sorted using the bitonic sorting algorithm.
+ *
+ * @param size the number of values to create.
+ * @param seed a random number generator seed value
+ * @return the root of the (sub) tree.
+ **/
+ static Value createTree(int size, int seed) {
+ if (size > 1) {
+ seed = random(seed);
+ int next_val = seed % RANGE;
+
+ Value retval = new Value(next_val);
+ retval.left = createTree(size/2, seed);
+ retval.right = createTree(size/2, skiprand(seed, size+1));
+ return retval;
+ } else {
+ return null;
+ }
+ }
+
+
+ /**
+ * Perform a bitonic sort based upon the Bilardi and Nicolau algorithm.
+ *
+ * @param spr_val the "spare" value in the algorithm.
+ * @param direction the direction of the sort (forward or backward)
+ * @return the new "spare" value.
+ **/
+ int bisort(int spr_val, bool direction) {
+ if (left is null) {
+ if ((value > spr_val) ^ direction) {
+ int tmpval = spr_val;
+ spr_val = value;
+ value = tmpval;
+ }
+ } else {
+ int val = value;
+ value = left.bisort(val, direction);
+ bool ndir = !direction;
+ spr_val = right.bisort(spr_val, ndir);
+ spr_val = bimerge(spr_val, direction);
+ }
+ return spr_val;
+ }
+
+
+ /**
+ * Perform the merge part of the bitonic sort. The merge part does
+ * the actualy sorting.
+ * @param spr_val the "spare" value in the algorithm.
+ * @param direction the direction of the sort (forward or backward)
+ * @return the new "spare" value
+ **/
+ int bimerge(int spr_val, bool direction) {
+ int rv = value;
+ Value pl = left;
+ Value pr = right;
+
+ bool rightexchange = (rv > spr_val) ^ direction;
+ if (rightexchange) {
+ value = spr_val;
+ spr_val = rv;
+ }
+
+ while (pl !is null) {
+ int lv = pl.value;
+ Value pll = pl.left;
+ Value plr = pl.right;
+ rv = pr.value;
+ Value prl = pr.left;
+ Value prr = pr.right;
+
+ bool elementexchange = (lv > rv) ^ direction;
+ if (rightexchange) {
+ if (elementexchange) {
+ pl.swapValRight(pr);
+ pl = pll;
+ pr = prl;
+ } else {
+ pl = plr;
+ pr = prr;
+ }
+ } else {
+ if (elementexchange) {
+ pl.swapValLeft(pr);
+ pl = plr;
+ pr = prr;
+ } else {
+ pl = pll;
+ pr = prl;
+ }
+ }
+ }
+
+ if (left !is null) {
+ value = left.bimerge(value, direction);
+ spr_val = right.bimerge(spr_val, direction);
+ }
+ return spr_val;
+ }
+
+
+ /**
+ * Swap the values and the right subtrees.
+ * @param n the other subtree involved in the swap.
+ **/
+ void swapValRight(Value n) {
+ int tmpv = n.value;
+ Value tmpr = n.right;
+
+ n.value = value;
+ n.right = right;
+
+ value = tmpv;
+ right = tmpr;
+ }
+
+
+ /**
+ * Swap the values and the left subtrees.
+ * @param n the other subtree involved in the swap.
+ **/
+ void swapValLeft(Value n) {
+ int tmpv = n.value;
+ Value tmpl = n.left;
+
+ n.value = value;
+ n.left = left;
+
+ value = tmpv;
+ left = tmpl;
+ }
+
+
+ /**
+ * Print out the nodes in the binary tree in infix order.
+ **/
+ void inOrder() {
+ if (left !is null)
+ left.inOrder();
+ printf("%d\n", value);
+ if (right !is null)
+ right.inOrder();
+ }
+
+
+ /**
+ * A random generator. The original Olden benchmark uses its
+ * own random generator. We use the same one in the Java version.
+ * @return the next random number in the sequence.
+ **/
+ private static int mult(int p, int q) {
+ int p1 = p / CONST_m1;
+ int p0 = p % CONST_m1;
+ int q1 = q / CONST_m1;
+ int q0 = q % CONST_m1;
+ return ((p0 * q1 + p1 * q0) % CONST_m1) * CONST_m1 + p0 * q0;
+ }
+
+
+ /**
+ * A routine to skip the next <i>n</i> random numbers.
+ * @param seed the current random no. seed
+ * @param n the number of numbers to skip
+ **/
+ private static int skiprand(int seed, int n) {
+ for (; n != 0; n--)
+ seed = random(seed);
+ return seed;
+ }
+
+
+ /**
+ * Return a random number based upon the seed value.
+ * @param seed the random number seed value
+ * @return a random number based upon the seed value.
+ **/
+ static int random(int seed) {
+ int tmp = mult(seed, CONST_b) + 1;
+ return tmp;
+ }
+}
+
+
+final public class BiSort {
+ private static int size = 0; // The number of values to sort.
+ private static bool printMsgs = false; // Print information messages
+ private static bool printResults = false; // Print the tree after each step
+
+
+ /**
+ * The main routine which creates a tree and sorts it a couple of times.
+ * @param args the command line arguments
+ **/
+ public static final void main(char[][] args) {
+ parseCmdLine(args);
+
+ if (printMsgs)
+ printf("Bisort with %d values\n", nextPow2(size+1) - 1);
+
+ auto start2 = myclock();
+ Value tree = Value.createTree(size, 12345768);
+ auto end2 = myclock();
+ int sval = Value.random(245867) % Value.RANGE;
+
+ if (printResults) {
+ tree.inOrder();
+ printf("%d\n", sval);
+ }
+
+ if (printMsgs)
+ printf("Beginning bitonic sort\n");
+
+ auto start0 = myclock();
+ sval = tree.bisort(sval, Value.FORWARD);
+ auto end0 = myclock();
+
+ if (printResults) {
+ tree.inOrder();
+ printf("%d\n", sval);
+ }
+
+ auto start1 = myclock();
+ sval = tree.bisort(sval, Value.BACKWARD);
+ auto end1 = myclock();
+
+ if (printResults) {
+ tree.inOrder();
+ printf("%d\n", sval);
+ }
+
+ if (printMsgs) {
+ printf("Creation time: %f\n", end2 - start2);
+ printf("Time to sort forward = %f\n", end0 - start0);
+ printf("Time to sort backward = %f\n", end1 - start1);
+ printf("Total: %f\n", end1 - start0);
+ printf("Done!\n");
+ }
+ }
+
+
+ /**
+ * Parse the command line options.
+ * @param args the command line options.
+ **/
+ private static final void parseCmdLine(char[][] args) {
+ int i = 1;
+ char[] arg;
+
+ while (i < args.length && args[i][0] == '-') {
+ arg = args[i++];
+
+ // check for options that require arguments
+ if (arg == "-s") {
+ if (i < args.length) {
+ size = toInt(args[i++]);
+ } else {
+ throw new Exception("-l requires the number of levels");
+ }
+ } else if (arg == "-m") {
+ printMsgs = true;
+ } else if (arg == "-p") {
+ printResults = true;
+ } else if (arg == "-h") {
+ usage();
+ }
+ }
+ if (size == 0)
+ usage();
+ }
+
+
+ /**
+ * The usage routine which describes the program options.
+ **/
+ private static final void usage() {
+ fprintf(stderr, "usage: bisort_d -s <size> [-p] [-i] [-h]\n");
+ fprintf(stderr, " -s the number of values to sort\n");
+ fprintf(stderr, " -m (print informative messages)\n");
+ fprintf(stderr, " -p (print the binary tree after each step)\n");
+ fprintf(stderr, " -h (print this message)\n");
+ exit(0);
+ }
+
+
+ private static /*unsigned*/ int nextPow2(/*unsigned*/ int x) {
+ if (x < 0)
+ throw new Exception("x must be >= 0");
+ x -= 1;
+ x |= x >> 1;
+ x |= x >> 2;
+ x |= x >> 4;
+ x |= x >> 8;
+ x |= x >> 16;
+ return x + 1;
+ }
+}
+
+
+void main(char[][] args) {
+ BiSort.main(args);
+}
--- /dev/null
+// Written by Leandro Lucarella
+//
+// The goal of this program is to do very CPU intensive work in threads
+
+import tango.core.Thread: Thread;
+import tango.core.sync.Atomic: atomicStore, atomicLoad, atomicAdd;
+import tango.io.device.File: File;
+import tango.util.digest.Sha512: Sha512;
+import tango.util.Convert: to;
+
+auto N = 100;
+auto NT = 2;
+ubyte[] BYTES;
+int running; // Atomic
+
+void main(char[][] args)
+{
+ auto fname = args[0];
+ if (args.length > 3)
+ fname = args[3];
+ if (args.length > 2)
+ NT = to!(int)(args[2]);
+ if (args.length > 1)
+ N = to!(int)(args[1]);
+ N /= NT;
+ atomicStore(running, NT);
+ BYTES = cast(ubyte[]) File.get(fname);
+ auto threads = new Thread[NT];
+ foreach(ref thread; threads) {
+ thread = new Thread(&doSha);
+ thread.start();
+ }
+ while (atomicLoad(running)) {
+ auto a = new void[](BYTES.length / 4);
+ a[] = cast(void[]) BYTES[];
+ Thread.yield();
+ }
+ foreach(thread; threads)
+ thread.join();
+}
+
+void doSha()
+{
+ for (size_t i = 0; i < N; i++) {
+ auto sha = new Sha512;
+ sha.update(BYTES);
+ }
+ atomicAdd(running, -1);
+}
+
--- /dev/null
+// Written by Leandro Lucarella
+//
+// The goal of this program is to do very CPU intensive work in threads
+
+import tango.core.Thread: Thread;
+import tango.core.sync.Atomic: atomicStore, atomicLoad, atomicAdd;
+import tango.io.device.File: File;
+import tango.util.digest.Sha512: Sha512;
+import tango.util.Convert: to;
+
+auto N = 100;
+auto NT = 2;
+ubyte[] BYTES;
+int running; // Atomic
+
+void main(char[][] args)
+{
+ auto fname = args[0];
+ if (args.length > 3)
+ fname = args[3];
+ if (args.length > 2)
+ NT = to!(int)(args[2]);
+ if (args.length > 1)
+ N = to!(int)(args[1]);
+ N /= NT;
+ atomicStore(running, NT);
+ BYTES = cast(ubyte[]) File.get(fname);
+ auto threads = new Thread[NT];
+ foreach(ref thread; threads) {
+ thread = new Thread(&doSha);
+ thread.start();
+ }
+ while (atomicLoad(running)) {
+ auto a = new void[](BYTES.length / 4);
+ a[] = cast(void[]) BYTES[];
+ Thread.yield();
+ }
+ foreach(thread; threads)
+ thread.join();
+}
+
+void doSha()
+{
+ auto sha = new Sha512;
+ for (size_t i = 0; i < N; i++) {
+ sha.update(BYTES);
+ }
+ atomicAdd(running, -1);
+}
+
--- /dev/null
+/**
+ * Java implementation of the <tt>em3d</tt> Olden benchmark. This Olden
+ * benchmark models the propagation of electromagnetic waves through
+ * objects in 3 dimensions. It is a simple computation on an irregular
+ * bipartite graph containing nodes representing electric and magnetic
+ * field values.
+ *
+ * <p><cite>
+ * D. Culler, A. Dusseau, S. Goldstein, A. Krishnamurthy, S. Lumetta, T. von
+ * Eicken and K. Yelick. "Parallel Programming in Split-C". Supercomputing
+ * 1993, pages 262-273.
+ * </cite>
+ *
+ * Java code converted to D by leonardo maffi, V.1.0, Oct 25 2009.
+ *
+ * Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.
+ * Downloaded from http://www.fantascienza.net/leonardo/js/
+ * http://www.fantascienza.net/leonardo/js/dolden_em3d.zip
+ **/
+
+
+version (Tango) {
+ import tango.stdc.stdio: printf, sprintf;
+ import tango.stdc.stdlib: exit;
+ import tango.stdc.time: CLOCKS_PER_SEC, clock;
+
+ import Integer = tango.text.convert.Integer;
+ alias Integer.parse toInt;
+} else {
+ import std.c.stdio: printf, sprintf;
+ import std.c.stdlib: exit;
+ import std.c.time: CLOCKS_PER_SEC, clock;
+
+ version (D_Version2) {
+ import std.conv: to;
+ alias to!(int, char[]) toInt;
+ } else {
+ import std.conv: toInt;
+ }
+}
+
+
+double myclock() {
+ return clock() / cast(double)CLOCKS_PER_SEC;
+}
+
+
+/**
+Basic uniform random generator: Minimal Standard in Park and
+Miller (1988): "Random Number Generators: Good Ones Are Hard to
+Find", Comm. of the ACM, 31, 1192-1201.
+Parameters: m = 2^31-1, a=48271.
+
+Adapted from Pascal code by Jesper Lund:
+http://www.gnu-pascal.de/crystal/gpc/en/mail1390.html
+*/
+final class Random {
+ const int m = int.max;
+ const int a = 48_271;
+ const int q = m / a;
+ const int r = m % a;
+ public int seed;
+
+ this(int the_seed) {
+ this.seed = the_seed;
+ }
+
+ public double nextDouble() {
+ int k = seed / q;
+ seed = a * (seed - k * q) - r * k;
+ if (seed < 1)
+ seed += m;
+ return cast(double)seed / m;
+ }
+
+ public int nextInt(int max) {
+ int n = max + 1;
+ int k = cast(int)(n * this.nextDouble());
+ return (k == n) ? k - 1 : k;
+ }
+}
+
+
+/**
+ * This class implements nodes (both E- and H-nodes) of the EM graph. Sets
+ * up random neighbors and propagates field values among neighbors.
+ */
+final class Node
+{
+ /**
+ * The value of the node.
+ **/
+ double value;
+ /**
+ * The next node in the list.
+ **/
+ private Node next;
+ /**
+ * Array of nodes to which we send our value.
+ **/
+ Node[] toNodes;
+ /**
+ * Array of nodes from which we receive values.
+ **/
+ Node[] fromNodes;
+ /**
+ * Coefficients on the fromNodes edges
+ **/
+ double[] coeffs;
+ /**
+ * The number of fromNodes edges
+ **/
+ int fromCount;
+ /**
+ * Used to create the fromEdges - keeps track of the number of edges that have
+ * been added
+ **/
+ int fromLength;
+
+ /**
+ * A random number generator.
+ **/
+ private static Random rand;
+
+ /**
+ * Initialize the random number generator
+ **/
+ public static void initSeed(long seed)
+ {
+ rand = new Random(seed);
+ }
+
+ /**
+ * Constructor for a node with given `degree'. The value of the
+ * node is initialized to a random value.
+ **/
+ this(int degree)
+ {
+ value = rand.nextDouble();
+ // create empty array for holding toNodes
+ toNodes = new Node[degree];
+
+ next = null;
+ fromNodes = null;
+ coeffs = null;
+ fromCount = 0;
+ fromLength = 0;
+ }
+
+ /**
+ * Create the linked list of E or H nodes. We create a table which is used
+ * later to create links among the nodes.
+ * @param size the no. of nodes to create
+ * @param degree the out degree of each node
+ * @return a table containing all the nodes.
+ **/
+ static Node[] fillTable(int size, int degree)
+ {
+ Node[] table = new Node[size];
+
+ Node prevNode = new Node(degree);
+ table[0] = prevNode;
+ for (int i = 1; i < size; i++) {
+ Node curNode = new Node(degree);
+ table[i] = curNode;
+ prevNode.next = curNode;
+ prevNode = curNode;
+ }
+ return table;
+ }
+
+ /**
+ * Create unique `degree' neighbors from the nodes given in nodeTable.
+ * We do this by selecting a random node from the give nodeTable to
+ * be neighbor. If this neighbor has been previously selected, then
+ * a different random neighbor is chosen.
+ * @param nodeTable the list of nodes to choose from.
+ **/
+ void makeUniqueNeighbors(Node[] nodeTable)
+ {
+ for (int filled = 0; filled < toNodes.length; filled++) {
+ int k;
+ Node otherNode;
+
+ do {
+ // generate a random number in the correct range
+ int index = rand.nextInt(nodeTable.length - 1);
+
+ // find a node with the random index in the given table
+ otherNode = nodeTable[index];
+
+ for (k = 0; k < filled; k++) {
+ if (otherNode == toNodes[k]) break; // fixed a bug of the original Java code
+ }
+ } while (k < filled);
+
+ // other node is definitely unique among "filled" toNodes
+ toNodes[filled] = otherNode;
+
+ // update fromCount for the other node
+ otherNode.fromCount++;
+ }
+ }
+
+ /**
+ * Allocate the right number of FromNodes for this node. This
+ * step can only happen once we know the right number of from nodes
+ * to allocate. Can be done after unique neighbors are created and known.
+ *
+ * It also initializes random coefficients on the edges.
+ **/
+ void makeFromNodes()
+ {
+ fromNodes = new Node[fromCount]; // nodes fill be filled in later
+ coeffs = new double[fromCount];
+ }
+
+ /**
+ * Fill in the fromNode field in "other" nodes which are pointed to
+ * by this node.
+ **/
+ void updateFromNodes()
+ {
+ for (int i = 0; i < toNodes.length; i++) {
+ Node otherNode = toNodes[i];
+ int count = otherNode.fromLength++;
+ otherNode.fromNodes[count] = this;
+ otherNode.coeffs[count] = rand.nextDouble();
+ }
+ }
+
+ /**
+ * Get the new value of the current node based on its neighboring
+ * from_nodes and coefficients.
+ **/
+ void computeNewValue()
+ {
+ for (int i = 0; i < fromCount; i++) {
+ value -= coeffs[i] * fromNodes[i].value;
+ }
+ }
+
+ int opApply(int delegate(ref Node) dg) {
+ int result;
+ for (Node current = this; current !is null; current = current.next) {
+ result = dg(current);
+ if (result)
+ break;
+ }
+ return result;
+ }
+
+ public char* toCString()
+ {
+ static char[256] repr;
+ sprintf(repr.ptr, "value %.17f, from_count %d", value, fromCount);
+ return repr.ptr;
+ }
+}
+
+
+/**
+ * A class that represents the irregular bipartite graph used in
+ * EM3D. The graph contains two linked structures that represent the
+ * E nodes and the N nodes in the application.
+ **/
+final class BiGraph
+{
+ /**
+ * Nodes that represent the electrical field.
+ **/
+ Node eNodes;
+ /**
+ * Nodes that representhe the magnetic field.
+ **/
+ Node hNodes;
+
+ /**
+ * Construct the bipartite graph.
+ * @param e the nodes representing the electric fields
+ * @param h the nodes representing the magnetic fields
+ **/
+ this(Node e, Node h)
+ {
+ eNodes = e;
+ hNodes = h;
+ }
+
+ /**
+ * Create the bi graph that contains the linked list of
+ * e and h nodes.
+ * @param numNodes the number of nodes to create
+ * @param numDegree the out-degree of each node
+ * @param verbose should we print out runtime messages
+ * @return the bi graph that we've created.
+ **/
+ static BiGraph create(int numNodes, int numDegree, bool verbose)
+ {
+ Node.initSeed(783);
+
+ // making nodes (we create a table)
+ if (verbose) printf("making nodes (tables in orig. version)\n");
+ Node[] hTable = Node.fillTable(numNodes, numDegree);
+ Node[] eTable = Node.fillTable(numNodes, numDegree);
+
+ // making neighbors
+ if (verbose) printf("updating from and coeffs\n");
+ foreach (n; hTable[0])
+ n.makeUniqueNeighbors(eTable);
+ foreach (n; eTable[0])
+ n.makeUniqueNeighbors(hTable);
+
+ // Create the fromNodes and coeff field
+ if (verbose) printf("filling from fields\n");
+ foreach (n; hTable[0])
+ n.makeFromNodes();
+ foreach (n; eTable[0])
+ n.makeFromNodes();
+
+ // Update the fromNodes
+ foreach (n; hTable[0])
+ n.updateFromNodes();
+ foreach (n; eTable[0])
+ n.updateFromNodes();
+
+ BiGraph g = new BiGraph(eTable[0], hTable[0]);
+ return g;
+ }
+
+ /**
+ * Update the field values of e-nodes based on the values of
+ * neighboring h-nodes and vice-versa.
+ **/
+ void compute() {
+ foreach (n; eNodes)
+ n.computeNewValue();
+ foreach (n; hNodes)
+ n.computeNewValue();
+ }
+
+ /**
+ * Print out the values of the e and h nodes.
+ **/
+ public void print() {
+ foreach (n; eNodes)
+ printf("E: %s\n", n.toCString());
+ foreach (n; hNodes)
+ printf("H: %s\n", n.toCString());
+ printf("\n");
+ }
+}
+
+public class Em3d1
+{
+ /**
+ * The number of nodes (E and H)
+ **/
+ private static int numNodes = 0;
+ /**
+ * The out-degree of each node.
+ **/
+ private static int numDegree = 0;
+ /**
+ * The number of compute iterations
+ **/
+ private static int numIter = 1;
+ /**
+ * Should we print the results and other runtime messages
+ **/
+ private static bool printResult = false;
+ /**
+ * Print information messages?
+ **/
+ private static bool printMsgs = false;
+
+ /**
+ * The main roitine that creates the irregular, linked data structure
+ * that represents the electric and magnetic fields and propagates the
+ * waves through the graph.
+ * @param args the command line arguments
+ **/
+ public static final void main(char[] args[])
+ {
+ parseCmdLine(args);
+
+ if (printMsgs)
+ printf("Initializing em3d random graph...\n");
+ auto start0 = myclock();
+ BiGraph graph = BiGraph.create(numNodes, numDegree, printResult);
+ auto end0 = myclock();
+
+ // compute a single iteration of electro-magnetic propagation
+ if (printMsgs)
+ printf("Propagating field values for %d iteration(s)...\n", numIter);
+ auto start1 = myclock();
+ for (int i = 0; i < numIter; i++) {
+ graph.compute();
+ }
+ auto end1 = myclock();
+
+ // print current field values
+ if (printResult)
+ graph.print();
+
+ if (printMsgs) {
+ printf("EM3D build time %.2f\n", end0 - start0);
+ printf("EM3D compute time %.2f\n", end1 - start1);
+ printf("EM3D total time %.2f\n", end1 - start0);
+ printf("Done!\n");
+ }
+ }
+
+
+ /**
+ * Parse the command line options.
+ * @param args the command line options.
+ **/
+ private static final void parseCmdLine(char[] args[])
+ {
+ int i = 1;
+ char[] arg;
+
+ while (i < args.length && args[i][0] == '-') {
+ arg = args[i++];
+
+ // check for options that require arguments
+ if (arg == "-n") {
+ if (i < args.length) {
+ numNodes = toInt(args[i++]);
+ } else throw new Exception("-n requires the number of nodes");
+ } else if (arg == "-d") {
+ if (i < args.length) {
+ numDegree = toInt(args[i++]);
+ } else throw new Exception("-d requires the out degree");
+ } else if (arg == "-i") {
+ if (i < args.length) {
+ numIter = toInt(args[i++]);
+ } else throw new Exception("-i requires the number of iterations");
+ } else if (arg == "-p") {
+ printResult = true;
+ } else if (arg == "-m") {
+ printMsgs = true;
+ } else if (arg == "-h") {
+ usage();
+ }
+ }
+ if (numNodes == 0 || numDegree == 0) usage();
+ }
+
+ /**
+ * The usage routine which describes the program options.
+ **/
+ private static final void usage() {
+ printf("usage: em3d -n <nodes> -d <degree> [-p] [-m] [-h]\n");
+ printf(" -n the number of nodes\n");
+ printf(" -d the out-degree of each node\n");
+ printf(" -i the number of iterations\n");
+ printf(" -p (print detailed results)\n");
+ printf(" -m (print informative messages)\n");
+ printf(" -h (this message)\n");
+ exit(0);
+ }
+}
+
+
+void main(char[][] args) {
+ Em3d1.main(args);
+}
+++ /dev/null
-// Written by Kevin Bealer <kevinbealer@gmail.com>
-// Found at http://www.digitalmars.com/webnews/newsgroups.php?art_group=digitalmars.D.announce&article_id=6978
-// Sightly modified by Leandro Lucarella <llucax@gmail.com>
-// (changed not to print anything and lower the total iterations; ported to
-// Tango)
-
-// Total residency should be ~160 MiB, but it usually increases a lot because
-// of false positives (probably in the static memory area)
-
-import tango.math.random.Random;
-
-const N = 2_000_000;
-const L = 20;
-const I = 50; // original: 200
-
-int main(char[][] args)
-{
- int[][] stuff;
-
- stuff.length = L;
-
- auto rand = new Random();
-
- for(int i = 0; i < I; i++) {
- int[] arr = new int[N];
-
- for(int j = 0; j < arr.length; j++) {
- rand(arr[j]);
- }
-
- int zig = i;
- if (zig >= stuff.length)
- zig = rand.uniform!(int) % stuff.length;
-
- stuff[zig] = arr;
- }
-
- return 0;
-}
-
import tango.math.random.Random;
-const IT = 50_000;
+const IT = 125; // number of iterations, each creates an object
+const BYTES = 1_000_000; // ~1MiB per object
+const N = 50; // ~50MiB of initial objects
+
+class C
+{
+ C c; // makes the compiler not set NO_SCAN
+ long[BYTES/long.sizeof] data;
+}
void main() {
- // The real memory use, ~55 KiB (original: ~20 MiB)
- uint[] data;
- data.length = 10_000; // original: 5_000_000
auto rand = new Random();
- foreach (ref x; data)
- rand(x);
+ C[] objs;
+ objs.length = N;
+ foreach (ref o; objs) {
+ o = new C;
+ foreach (ref x; o.data)
+ rand(x);
+ }
for (int i = 0; i < IT; ++i) {
- // simulate reading a few kb of data (14 KiB +/- 10 KiB)
- uint[] incoming;
- incoming.length = 1000 + rand.uniform!(uint) % 5000;
- foreach (ref x; incoming)
+ C o = new C;
+ foreach (ref x; o.data)
rand(x);
// do something with the data...
}
--- /dev/null
+// Written by Andrey Khropov <andkhropov_nosp@m_mtu-net.ru>
+// Found at http://www.digitalmars.com/webnews/newsgroups.php?art_group=digitalmars.D&article_id=43991
+// Modified by Leandro Lucarella
+// (ported to Tango, fixed some stylistic issues)
+
+import tango.util.Convert;
+
+alias char[] string;
+
+int main(string[] args)
+{
+ int N = args.length > 1 ? to!(int)(args[1]) : 1;
+
+ int minDepth = 4;
+ int maxDepth = (minDepth + 2) > N ? minDepth + 2 : N;
+ int stretchDepth = maxDepth + 1;
+
+ int check = TreeNode.BottomUpTree(0, stretchDepth).ItemCheck;
+ TreeNode longLivedTree = TreeNode.BottomUpTree(0, maxDepth);
+
+ for (int depth = minDepth; depth <= maxDepth; depth += 2) {
+ int iterations = 1 << (maxDepth - depth + minDepth);
+ check = 0;
+
+ for (int i = 1; i <= iterations; i++) {
+ check += TreeNode.BottomUpTree(i, depth).ItemCheck;
+ check += TreeNode.BottomUpTree(-i, depth).ItemCheck;
+ }
+
+ }
+
+ return 0;
+}
+
+class TreeNode
+{
+ TreeNode left, right;
+ int item;
+
+ this(int item, TreeNode left = null, TreeNode right = null)
+ {
+ this.item = item;
+ this.left = left;
+ this.right = right;
+ }
+
+ static TreeNode BottomUpTree(int item, int depth)
+ {
+ if (depth > 0)
+ return new TreeNode(item,
+ BottomUpTree(2 * item - 1, depth - 1),
+ BottomUpTree(2 * item, depth - 1));
+ return new TreeNode(item);
+ }
+
+ int ItemCheck()
+ {
+ if (left)
+ return item + left.ItemCheck() - right.ItemCheck();
+ return item;
+ }
+}
+
+++ /dev/null
-// Written by Andrey Khropov <andkhropov_nosp@m_mtu-net.ru>
-// Found at http://www.digitalmars.com/webnews/newsgroups.php?art_group=digitalmars.D&article_id=43991
-// Modified by Leandro Lucarella
-// (ported to Tango)
-
-import tango.util.Convert;
-
-alias char[] string;
-
-int main(string[] args)
-{
- int N = args.length > 1 ? to!(int)(args[1]) : 1;
-
- int minDepth = 4;
- int maxDepth = (minDepth + 2) > N ? minDepth + 2 : N;
- int stretchDepth = maxDepth + 1;
-
- TreeNode stretchTree = TreeNode.BottomUpTree(0, stretchDepth);
-
- TreeNode longLivedTree = TreeNode.BottomUpTree(0, maxDepth);
-
- int depth;
- for(depth = minDepth; depth <= maxDepth; depth += 2)
- {
- int check, iterations = 1 << (maxDepth - depth + minDepth);
-
- for(int i = 1; i <= iterations; i++)
- {
- auto tempTree = TreeNode.BottomUpTree(i, depth);
- check += tempTree.ItemCheck;
- //delete tempTree;
-
- tempTree = TreeNode.BottomUpTree(-i, depth);
- check += tempTree.ItemCheck;
- //delete tempTree;
- }
-
- }
-
- return 0;
-}
-
-class TreeNode
-{
-public:
- this(int item, TreeNode left = null, TreeNode right = null)
- {
- this.item = item;
- this.left = left;
- this.right = right;
- }
-
- static TreeNode BottomUpTree(int item, int depth)
- {
- if(depth > 0)
- return new TreeNode(item
- ,BottomUpTree(2 * item - 1, depth - 1)
- ,BottomUpTree(2 * item, depth - 1));
- return new TreeNode(item);
- }
-
- int ItemCheck()
- {
- if(left)
- return item + left.ItemCheck() - right.ItemCheck();
- return item;
- }
-private:
- TreeNode left, right;
- int item;
-}
-
+++ /dev/null
-// Written by Piotr Modzelewski <http://petermodzelewski.blogspot.com/>
-// Found at http://www.dsource.org/projects/tango/wiki/GCBenchmark
-
-class TreeNode {
- int item;
- TreeNode left, right;
-
- this(int item, TreeNode left=null, TreeNode right=null) {
- this.item = item;
- this.left = left;
- this.right = right;
- }
-
- int check() {
- return left is null ? item : item + left.check - right.check;
- }
-}
-
-TreeNode makeTree(int item, int depth) {
- if (depth > 0)
- return new TreeNode(item, makeTree(2*item-1, depth-1),
- makeTree(2*item, depth-1));
- else
- return new TreeNode(item);
-}
-
-void main(char[][] args) {
- const minDepth = 5;
- int n = 15;
- int maxDepth = (minDepth + 2) > n ? minDepth + 2 : n;
-
- int check = makeTree(0, maxDepth + 1).check;
-
- auto longLivedTree = makeTree(0, maxDepth);
-
- for (int depth = minDepth; depth <= maxDepth; depth += 2) {
- int iterations = 1 << (maxDepth - depth + minDepth);
- check = 0;
-
- for (int i = 1; i <= iterations; i++)
- check += (makeTree(i, depth)).check
- + (makeTree(-i, depth)).check;
- }
-
-}
-
--- /dev/null
+/**\r
+A D implementation of the "tsp" Olden benchmark, the traveling\r
+salesman problem.\r
+\r
+R. Karp, "Probabilistic analysis of partitioning algorithms for the\r
+traveling-salesman problem in the plane." Mathematics of Operations Research\r
+2(3):209-224, August 1977.\r
+\r
+Converted to D and optimized by leonardo maffi, V.1.0, Oct 29 2009.\r
+\r
+Removed output unless an option is passed by Leandro Lucarella, 2010-08-04.\r
+Downloaded from http://www.fantascienza.net/leonardo/js/\r
+ http://www.fantascienza.net/leonardo/js/dolden_tsp.zip\r
+*/\r
+\r
+version (Tango) {\r
+ import tango.stdc.stdio: printf, fprintf, stderr;\r
+ import tango.stdc.time: CLOCKS_PER_SEC, clock;\r
+ import tango.math.Math: sqrt, log;\r
+\r
+ import Integer = tango.text.convert.Integer;\r
+ alias Integer.parse toInt;\r
+} else {\r
+ import std.c.stdio: printf, fprintf, stderr;\r
+ import std.c.time: CLOCKS_PER_SEC, clock;\r
+ import std.math: sqrt, log;\r
+\r
+ version (D_Version2) {\r
+ import std.conv: to;\r
+ alias to!(int, char[]) toInt;\r
+ } else {\r
+ import std.conv: toInt;\r
+ }\r
+}\r
+\r
+\r
+double myclock() {\r
+ return clock() / cast(double)CLOCKS_PER_SEC;\r
+}\r
+\r
+\r
+/**\r
+Basic uniform random generator: Minimal Standard in Park and\r
+Miller (1988): "Random Number Generators: Good Ones Are Hard to\r
+Find", Comm. of the ACM, 31, 1192-1201.\r
+Parameters: m = 2^31-1, a=48271.\r
+\r
+Very simple portable random number generator adapted\r
+from Pascal code by Jesper Lund:\r
+http://www.gnu-pascal.de/crystal/gpc/en/mail1390.html\r
+*/\r
+final class Random {\r
+ const int m = int.max;\r
+ const int a = 48271;\r
+ const int q = m / a;\r
+ const int r = m % a;\r
+\r
+ public int seed;\r
+\r
+ this(int the_seed) {\r
+ this.seed = the_seed;\r
+ }\r
+\r
+ public double nextDouble() {\r
+ int k = seed / q;\r
+ seed = a * (seed - k * q) - r * k;\r
+ if (seed < 1)\r
+ seed += m;\r
+ return cast(double)seed / m;\r
+ }\r
+\r
+ public int nextInt(int max) {\r
+ int n = max + 1;\r
+ int k = cast(int)(n * this.nextDouble());\r
+ return (k == n) ? k - 1 : k;\r
+ }\r
+}\r
+\r
+\r
+\r
+/**\r
+* A class that represents a node in a binary tree. Each node represents\r
+* a city in the TSP benchmark.\r
+**/\r
+final class Tree {\r
+ /**\r
+ * The number of nodes (cities) in this subtree\r
+ **/\r
+ private int sz;\r
+\r
+ /**\r
+ * The coordinates that this node represents\r
+ **/\r
+ private double x, y;\r
+\r
+ /**\r
+ * Left and right child of tree\r
+ **/\r
+ private Tree left, right;\r
+\r
+ /**\r
+ * The next pointer in a linked list of nodes in the subtree. The list\r
+ * contains the order of the cities to visit.\r
+ **/\r
+ private Tree next;\r
+\r
+ /**\r
+ * The previous pointer in a linked list of nodes in the subtree. The list\r
+ * contains the order of the cities to visit.\r
+ **/\r
+ private Tree prev;\r
+\r
+ static Random rnd;\r
+\r
+ // used by the random number generator\r
+ private const double M_E2 = 7.3890560989306502274;\r
+ private const double M_E3 = 20.08553692318766774179;\r
+ private const double M_E6 = 403.42879349273512264299;\r
+ private const double M_E12 = 162754.79141900392083592475;\r
+\r
+ public static void initRnd(int seed) {\r
+ rnd = new Random(seed);\r
+ }\r
+\r
+ /**\r
+ * Construct a Tree node (a city) with the specified size\r
+ * @param size the number of nodes in the (sub)tree\r
+ * @param x the x coordinate of the city\r
+ * @param y the y coordinate of the city\r
+ * @param left the left subtree\r
+ * @param right the right subtree\r
+ **/\r
+ this(int size, double x, double y, Tree l, Tree r) {\r
+ sz = size;\r
+ this.x = x;\r
+ this.y = y;\r
+ left = l;\r
+ right = r;\r
+ next = null;\r
+ prev = null;\r
+ }\r
+\r
+ /**\r
+ * Find Euclidean distance between this node and the specified node.\r
+ * @param b the specified node\r
+ * @return the Euclidean distance between two tree nodes.\r
+ **/\r
+ double distance(Tree b) {\r
+ return sqrt((x - b.x) * (x - b.x) + (y - b.y) * (y - b.y));\r
+ }\r
+\r
+ /**\r
+ * Create a list of tree nodes. Requires root to be the tail of the list.\r
+ * Only fills in next field, not prev.\r
+ * @return the linked list of nodes\r
+ **/\r
+ Tree makeList() {\r
+ Tree myleft, myright, tleft, tright;\r
+ Tree retval = this;\r
+\r
+ // head of left list\r
+ if (left !is null)\r
+ myleft = left.makeList();\r
+ else\r
+ myleft = null;\r
+\r
+ // head of right list\r
+ if (right !is null)\r
+ myright = right.makeList();\r
+ else\r
+ myright = null;\r
+\r
+ if (myright !is null) {\r
+ retval = myright;\r
+ right.next = this;\r
+ }\r
+\r
+ if (myleft !is null) {\r
+ retval = myleft;\r
+ if (myright !is null)\r
+ left.next = myright;\r
+ else\r
+ left.next = this;\r
+ }\r
+ next = null;\r
+\r
+ return retval;\r
+ }\r
+\r
+ /**\r
+ * Reverse the linked list. Assumes that there is a dummy "prev"\r
+ * node at the beginning.\r
+ **/\r
+ void reverse() {\r
+ Tree prev = this.prev;\r
+ prev.next = null;\r
+ this.prev = null;\r
+ Tree back = this;\r
+ Tree tmp = this;\r
+\r
+ // reverse the list for the other nodes\r
+ Tree next;\r
+ for (Tree t = this.next; t !is null; back = t, t = next) {\r
+ next = t.next;\r
+ t.next = back;\r
+ back.prev = t;\r
+ }\r
+\r
+ // reverse the list for this node\r
+ tmp.next = prev;\r
+ prev.prev = tmp;\r
+ }\r
+\r
+ /**\r
+ * Use closest-point heuristic from Cormen, Leiserson, and Rivest.\r
+ * @return a\r
+ **/\r
+ Tree conquer() {\r
+ // create the list of nodes\r
+ Tree t = makeList();\r
+\r
+ // Create initial cycle\r
+ Tree cycle = t;\r
+ t = t.next;\r
+ cycle.next = cycle;\r
+ cycle.prev = cycle;\r
+\r
+ // Loop over remaining points\r
+ Tree donext;\r
+ for ( ; t !is null; t = donext) {\r
+ donext = t.next; /* value won't be around later */\r
+ Tree min = cycle;\r
+ double mindist = t.distance(cycle);\r
+ for (Tree tmp = cycle.next; tmp != cycle; tmp = tmp.next) {\r
+ double test = tmp.distance(t);\r
+ if (test < mindist) {\r
+ mindist = test;\r
+ min = tmp;\r
+ }\r
+ }\r
+\r
+ Tree next = min.next;\r
+ Tree prev = min.prev;\r
+\r
+ double mintonext = min.distance(next);\r
+ double mintoprev = min.distance(prev);\r
+ double ttonext = t.distance(next);\r
+ double ttoprev = t.distance(prev);\r
+\r
+ if ((ttoprev - mintoprev) < (ttonext - mintonext)) {\r
+ // insert between min and prev\r
+ prev.next = t;\r
+ t.next = min;\r
+ t.prev = prev;\r
+ min.prev = t;\r
+ } else {\r
+ next.prev = t;\r
+ t.next = next;\r
+ min.next = t;\r
+ t.prev = min;\r
+ }\r
+ }\r
+\r
+ return cycle;\r
+ }\r
+\r
+ /**\r
+ * Merge two cycles as per Karp.\r
+ * @param a a subtree with one cycle\r
+ * @param b a subtree with the other cycle\r
+ **/\r
+ Tree merge(Tree a, Tree b) {\r
+ // Compute location for first cycle\r
+ Tree min = a;\r
+ double mindist = distance(a);\r
+ Tree tmp = a;\r
+ for (a = a.next; a != tmp; a = a.next) {\r
+ double test = distance(a);\r
+ if (test < mindist) {\r
+ mindist = test;\r
+ min = a;\r
+ }\r
+ }\r
+\r
+ Tree next = min.next;\r
+ Tree prev = min.prev;\r
+ double mintonext = min.distance(next);\r
+ double mintoprev = min.distance(prev);\r
+ double ttonext = distance(next);\r
+ double ttoprev = distance(prev);\r
+\r
+ Tree p1, n1;\r
+ double tton1, ttop1;\r
+ if ((ttoprev - mintoprev) < (ttonext - mintonext)) {\r
+ // would insert between min and prev\r
+ p1 = prev;\r
+ n1 = min;\r
+ tton1 = mindist;\r
+ ttop1 = ttoprev;\r
+ } else {\r
+ // would insert between min and next\r
+ p1 = min;\r
+ n1 = next;\r
+ ttop1 = mindist;\r
+ tton1 = ttonext;\r
+ }\r
+\r
+ // Compute location for second cycle\r
+ min = b;\r
+ mindist = distance(b);\r
+ tmp = b;\r
+ for (b = b.next; b != tmp; b = b.next) {\r
+ double test = distance(b);\r
+ if (test < mindist) {\r
+ mindist = test;\r
+ min = b;\r
+ }\r
+ }\r
+\r
+ next = min.next;\r
+ prev = min.prev;\r
+ mintonext = min.distance(next);\r
+ mintoprev = min.distance(prev);\r
+ ttonext = this.distance(next);\r
+ ttoprev = this.distance(prev);\r
+\r
+ Tree p2, n2;\r
+ double tton2, ttop2;\r
+ if ((ttoprev - mintoprev) < (ttonext - mintonext)) {\r
+ // would insert between min and prev\r
+ p2 = prev;\r
+ n2 = min;\r
+ tton2 = mindist;\r
+ ttop2 = ttoprev;\r
+ } else {\r
+ // would insert between min andn ext\r
+ p2 = min;\r
+ n2 = next;\r
+ ttop2 = mindist;\r
+ tton2 = ttonext;\r
+ }\r
+\r
+ // Now we have 4 choices to complete:\r
+ // 1:t,p1 t,p2 n1,n2\r
+ // 2:t,p1 t,n2 n1,p2\r
+ // 3:t,n1 t,p2 p1,n2\r
+ // 4:t,n1 t,n2 p1,p2\r
+ double n1ton2 = n1.distance(n2);\r
+ double n1top2 = n1.distance(p2);\r
+ double p1ton2 = p1.distance(n2);\r
+ double p1top2 = p1.distance(p2);\r
+\r
+ mindist = ttop1 + ttop2 + n1ton2;\r
+ int choice = 1;\r
+\r
+ double test = ttop1 + tton2 + n1top2;\r
+ if (test < mindist) {\r
+ choice = 2;\r
+ mindist = test;\r
+ }\r
+\r
+ test = tton1 + ttop2 + p1ton2;\r
+ if (test < mindist) {\r
+ choice = 3;\r
+ mindist = test;\r
+ }\r
+\r
+ test = tton1 + tton2 + p1top2;\r
+ if (test < mindist)\r
+ choice = 4;\r
+\r
+ switch (choice) {\r
+ case 1:\r
+ // 1:p1,this this,p2 n2,n1 -- reverse 2!\r
+ n2.reverse();\r
+ p1.next = this;\r
+ this.prev = p1;\r
+ this.next = p2;\r
+ p2.prev = this;\r
+ n2.next = n1;\r
+ n1.prev = n2;\r
+ break;\r
+\r
+ case 2:\r
+ // 2:p1,this this,n2 p2,n1 -- OK\r
+ p1.next = this;\r
+ this.prev = p1;\r
+ this.next = n2;\r
+ n2.prev = this;\r
+ p2.next = n1;\r
+ n1.prev = p2;\r
+ break;\r
+\r
+ case 3:\r
+ // 3:p2,this this,n1 p1,n2 -- OK\r
+ p2.next = this;\r
+ this.prev = p2;\r
+ this.next = n1;\r
+ n1.prev = this;\r
+ p1.next = n2;\r
+ n2.prev = p1;\r
+ break;\r
+\r
+ default: // case 4:\r
+ // 4:n1,this this,n2 p2,p1 -- reverse 1!\r
+ n1.reverse();\r
+ n1.next = this;\r
+ this.prev = n1;\r
+ this.next = n2;\r
+ n2.prev = this;\r
+ p2.next = p1;\r
+ p1.prev = p2;\r
+ break;\r
+ }\r
+\r
+ return this;\r
+ } // end merge\r
+\r
+ /**\r
+ * Compute TSP for the tree t. Use conquer for problems <= sz\r
+ * @param sz the cutoff point for using conquer vs. merge\r
+ **/\r
+ Tree tsp(int sz) {\r
+ if (this.sz <= sz)\r
+ return conquer();\r
+ Tree leftval = left.tsp(sz);\r
+ Tree rightval = right.tsp(sz);\r
+ return merge(leftval, rightval);\r
+ }\r
+\r
+ /**\r
+ * Print the list of cities to visit from the current node. The\r
+ * list is the result of computing the TSP problem.\r
+ * The list for the root node (city) should contain every other node\r
+ * (city).\r
+ **/\r
+ void printVisitOrder() {\r
+ printf("x = %.15f y = %.15f\n", x, y);\r
+ for (Tree tmp = next; tmp != this; tmp = tmp.next)\r
+ printf("x = %.15f y = %.15f\n", tmp.x, tmp.y);\r
+ }\r
+\r
+ /**\r
+ * Computes the total length of the current tour.\r
+ **/\r
+ double tourLength() {\r
+ double total = 0.0;\r
+ Tree precedent = next;\r
+ Tree current = next.next;\r
+ if (current == this)\r
+ return total;\r
+\r
+ do {\r
+ total += current.distance(precedent);\r
+ precedent = current;\r
+ current = current.next;\r
+ } while (precedent != this);\r
+\r
+ total += current.distance(this);\r
+ return total;\r
+ }\r
+\r
+\r
+ // static methods ===============================================\r
+\r
+ /**\r
+ * Return an estimate of median of n values distributed in [min, max)\r
+ * @param min the minimum value\r
+ * @param max the maximum value\r
+ * @param n\r
+ * @return an estimate of median of n values distributed in [min, max)\r
+ **/\r
+ private static double median(double min, double max, int n) {\r
+ // get random value in [0.0, 1.0)\r
+ double t = rnd.nextDouble();\r
+\r
+ double retval;\r
+ if (t > 0.5)\r
+ retval = log(1.0 - (2.0 * (M_E12 - 1) * (t - 0.5) / M_E12)) / 12.0;\r
+ else\r
+ retval = -log(1.0 - (2.0 * (M_E12 - 1) * t / M_E12)) / 12.0;\r
+\r
+ // We now have something distributed on (-1.0, 1.0)\r
+ retval = (retval + 1.0) * (max - min) / 2.0;\r
+ return retval + min;\r
+ }\r
+\r
+ /**\r
+ * Get double uniformly distributed over [min,max)\r
+ * @return double uniformily distributed over [min,max)\r
+ **/\r
+ private static double uniform(double min, double max) {\r
+ // get random value between [0.0, 1.0)\r
+ double retval = rnd.nextDouble();\r
+ retval = retval * (max - min);\r
+ return retval + min;\r
+ }\r
+\r
+ /**\r
+ * Builds a 2D tree of n nodes in specified range with dir as primary\r
+ * axis (false for x, true for y)\r
+ *\r
+ * @param n the size of the subtree to create\r
+ * @param dir the primary axis\r
+ * @param min_x the minimum x coordinate\r
+ * @param max_x the maximum x coordinate\r
+ * @param min_y the minimum y coordinate\r
+ * @param max_y the maximum y coordinate\r
+ * @return a reference to the root of the subtree\r
+ **/\r
+ public static Tree buildTree(int n, bool dir, double min_x,\r
+ double max_x, double min_y, double max_y) {\r
+ if (n == 0)\r
+ return null;\r
+\r
+ Tree left, right;\r
+ double x, y;\r
+ if (dir) {\r
+ dir = !dir;\r
+ double med = median(min_x, max_x, n);\r
+ left = buildTree(n/2, dir, min_x, med, min_y, max_y);\r
+ right = buildTree(n/2, dir, med, max_x, min_y, max_y);\r
+ x = med;\r
+ y = uniform(min_y, max_y);\r
+ } else {\r
+ dir = !dir;\r
+ double med = median(min_y,max_y,n);\r
+ left = buildTree(n/2, dir, min_x, max_x, min_y, med);\r
+ right = buildTree(n/2, dir, min_x, max_x, med, max_y);\r
+ y = med;\r
+ x = uniform(min_x, max_x);\r
+ }\r
+\r
+ return new Tree(n, x, y, left, right);\r
+ }\r
+}\r
+\r
+\r
+public class TSP {\r
+ /**\r
+ * Number of cities in the problem.\r
+ **/\r
+ private static int cities;\r
+\r
+ /**\r
+ * Set to true if the result should be printed\r
+ **/\r
+ private static bool printResult = false;\r
+\r
+ /**\r
+ * Set to true to print informative messages\r
+ **/\r
+ private static bool printMsgs = false;\r
+\r
+ /**\r
+ * The main routine which creates a tree and traverses it.\r
+ * @param args the arguments to the program\r
+ **/\r
+ public static void main(char[] args[]) {\r
+ if (!parseCmdLine(args))\r
+ return;\r
+\r
+ if (printMsgs)\r
+ printf("Building tree of size: %d\n", nextPow2(cities+1) - 1);\r
+\r
+ Tree.initRnd(1);\r
+\r
+ auto t0 = myclock();\r
+ Tree t = Tree.buildTree(cities, false, 0.0, 1.0, 0.0, 1.0);\r
+\r
+ auto t1 = myclock();\r
+ t.tsp(150);\r
+ auto t2 = myclock();\r
+\r
+ if (printMsgs)\r
+ printf("Total tour length: %.15f\n", t.tourLength());\r
+ auto t3 = myclock();\r
+\r
+ if (printResult)\r
+ // if the user specifies, print the final result\r
+ t.printVisitOrder();\r
+\r
+ if (printMsgs) {\r
+ printf("Tsp build time: %.2f\n", t1 - t0);\r
+ printf("Tsp computing time: %.2f\n", t2 - t1);\r
+ printf("Tsp total time: %.2f\n", t3 - t0);\r
+ }\r
+ }\r
+\r
+ private static /*unsigned*/ int nextPow2(/*unsigned*/ int x) {\r
+ if (x < 0)\r
+ throw new Exception("x must be >= 0");\r
+ x = x - 1;\r
+ x = x | (x >> 1);\r
+ x = x | (x >> 2);\r
+ x = x | (x >> 4);\r
+ x = x | (x >> 8);\r
+ x = x | (x >> 16);\r
+ return x + 1;\r
+ }\r
+\r
+ /**\r
+ * Parse the command line options.\r
+ * @param args the command line options.\r
+ **/\r
+ private static final bool parseCmdLine(char[] args[]) {\r
+ int i = 1;\r
+\r
+ while (i < args.length && args[i][0] == '-') {\r
+ char[] arg = args[i++];\r
+\r
+ if (arg == "-c") {\r
+ if (i < args.length)\r
+ cities = toInt(args[i++]);\r
+ else\r
+ throw new Exception("-c requires the size of tree");\r
+ if (cities < 1)\r
+ throw new Exception("Number of cities must be > 0");\r
+ } else if (arg == "-p") {\r
+ printResult = true;\r
+ } else if (arg == "-m") {\r
+ printMsgs = true;\r
+ } else if (arg == "-h") {\r
+ return usage();\r
+ }\r
+ }\r
+\r
+ if (cities == 0)\r
+ return usage();\r
+\r
+ return true;\r
+ }\r
+\r
+ /**\r
+ * The usage routine which describes the program options.\r
+ **/\r
+ private static final bool usage() {\r
+ fprintf(stderr, "usage: tsp_d -c <num> [-p] [-m] [-h]\n");\r
+ fprintf(stderr, " -c number of cities (rounds up to the next power of 2 minus 1)\n");\r
+ fprintf(stderr, " -p (print the final result)\n");\r
+ fprintf(stderr, " -m (print informative messages)\n");\r
+ fprintf(stderr, " -h (print this message)\n");\r
+ return false;\r
+ }\r
+}\r
+\r
+\r
+void main(char[][] args) {\r
+ TSP.main(args);\r
+}\r
projects / software / dgc / dgcbench.git / commitdiff