From 226a056a64171b7464faaf97ed5d331ffc56478a Mon Sep 17 00:00:00 2001 From: Leandro Lucarella Date: Sun, 14 Nov 2010 00:39:16 -0300 Subject: [PATCH] Update micro benchmarks Remove some micro benchmarks that provided no added value, rename others to have shorter names and add a few; some Olden benchmarks and a couple of benchmarks to exercise concurrency. --- Makefile | 114 ++-- micro/bh.d | 1045 ++++++++++++++++++++++++++++++ micro/{big_arrays.d => bigarr.d} | 0 micro/bisort.d | 366 +++++++++++ micro/conalloc.d | 50 ++ micro/concpu.d | 50 ++ micro/em3d.d | 468 +++++++++++++ micro/{multicore.d => mcore.d} | 0 micro/rnd_data_2.d | 40 -- micro/{rnd_data.d => rnddata.d} | 28 +- micro/sbtree.d | 63 ++ micro/shootout_binarytrees.d | 72 -- micro/tree.d | 46 -- micro/tsp.d | 650 +++++++++++++++++++ 14 files changed, 2776 insertions(+), 216 deletions(-) create mode 100644 micro/bh.d rename micro/{big_arrays.d => bigarr.d} (100%) create mode 100644 micro/bisort.d create mode 100644 micro/conalloc.d create mode 100644 micro/concpu.d create mode 100644 micro/em3d.d rename micro/{multicore.d => mcore.d} (100%) delete mode 100644 micro/rnd_data_2.d rename micro/{rnd_data.d => rnddata.d} (52%) create mode 100644 micro/sbtree.d delete mode 100644 micro/shootout_binarytrees.d delete mode 100644 micro/tree.d create mode 100644 micro/tsp.d diff --git a/Makefile b/Makefile index 5ff8b62..f32edb6 100644 --- a/Makefile +++ b/Makefile @@ -1,10 +1,11 @@ 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 @@ -13,7 +14,7 @@ ifndef V 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'; @@ -24,46 +25,81 @@ endif # 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 @@ -73,17 +109,14 @@ $O/time/voronoi.t.csv $O/time/voronoi.s.csv \ 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 $( $@ ($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/% @@ -140,10 +165,3 @@ $O/stats/%.eps: $O/stats/%.c.csv $O/stats/%.a.csv $O/stats/%.h.csv \ 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 - diff --git a/micro/bh.d b/micro/bh.d new file mode 100644 index 0000000..a1c092a --- /dev/null +++ b/micro/bh.d @@ -0,0 +1,1045 @@ +/** +A D implementation of the _bh_ Olden benchmark. +The Olden benchmark implements the Barnes-Hut benchmark +that is decribed in: + +J. Barnes and P. Hut, "A hierarchical o(N log N) force-calculation algorithm", +Nature, 324:446-449, Dec. 1986 + +The original code in the Olden benchmark suite is derived from the +ftp://hubble.ifa.hawaii.edu/pub/barnes/treecode +source distributed by Barnes. + +Java code converted to D by leonardo maffi, V.1.0, Dec 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_bh.zip +*/ + + +version (Tango) { + import tango.stdc.stdio: printf, sprintf, fprintf, stderr; + import tango.stdc.stdlib: exit, malloc, free; + import tango.stdc.time: CLOCKS_PER_SEC, clock; + import tango.math.Math: sqrt, floor, PI, pow; + + import Integer = tango.text.convert.Integer; + alias Integer.parse toInt; +} else { + import std.c.stdio: printf, sprintf, fprintf, stderr; + import std.c.stdlib: exit, malloc, free; + import std.c.time: CLOCKS_PER_SEC, clock; + import std.math: sqrt, floor, PI, pow; + + 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 +*/ +struct Random { + const int m = int.max; + const int a = 48_271; + const int q = m / a; + const int r = m % a; + int seed; + + public double uniform(double min, double max) { + int k = seed / q; + seed = a * (seed - k * q) - r * k; + if (seed < 1) + seed += m; + double r = cast(double)seed / m; + return r * (max - min) + min; + } +} + + +interface Enumeration { + bool hasMoreElements(); + Object nextElement(); +} + + +/** +A class representing a three dimensional vector that implements +several math operations. To improve speed we implement the +vector as an array of doubles rather than use the exising +code in the java.util.Vec3 class. +*/ +final class Vec3 { + /// The number of dimensions in the vector + // Can't be changed because of crossProduct() + const int NDIM = 3; + + /// An array containing the values in the vector. + private double data[]; + + /// Construct an empty 3 dimensional vector for use in Barnes-Hut algorithm. + this() { + data.length = NDIM; + data[] = 0.0; + } + + /// Create a copy of the vector. Return a clone of the math vector + public Vec3 clone() { + Vec3 v = new Vec3; + v.data.length = NDIM; + v.data[] = this.data[]; + return v; + } + + /** + Return the value at the i'th index of the vector. + @param i the vector index + @return the value at the i'th index of the vector. + */ + final double value(int i) { + return data[i]; + } + + /** + Set the value of the i'th index of the vector. + @param i the vector index + @param v the value to store + */ + final void value(int i, double v) { + data[i] = v; + } + + /** + Set one of the dimensions of the vector to 1.0 + param j the dimension to set. + */ + final void unit(int j) { + for (int i = 0; i < NDIM; i++) + data[i] = (i == j) ? 1.0 : 0.0; + } + + /** + Add two vectors and the result is placed in this vector. + @param u the other operand of the addition + */ + final void addition(Vec3 u) { + for (int i = 0; i < NDIM; i++) + data[i] += u.data[i]; + } + + /** + * Subtract two vectors and the result is placed in this vector. + * This vector contain the first operand. + * @param u the other operand of the subtraction. + **/ + final void subtraction(Vec3 u) { + for (int i = 0; i < NDIM; i++) + data[i] -= u.data[i]; + } + + /** + Subtract two vectors and the result is placed in this vector. + @param u the first operand of the subtraction. + @param v the second opernd of the subtraction + */ + final void subtraction(Vec3 u, Vec3 v) { + for (int i = 0; i < NDIM; i++) + data[i] = u.data[i] - v.data[i]; + } + + /** + Multiply the vector times a scalar. + @param s the scalar value + **/ + final void multScalar(double s) { + for (int i = 0; i < NDIM; i++) + data[i] *= s; + } + + /** + * Multiply the vector times a scalar and place the result in this vector. + * @param u the vector + * @param s the scalar value + **/ + final void multScalar(Vec3 u, double s) { + for (int i = 0; i < NDIM; i++) + data[i] = u.data[i] * s; + } + + /** + * Divide each element of the vector by a scalar value. + * @param s the scalar value. + **/ + final void divScalar(double s) { + for (int i = 0; i < NDIM; i++) + data[i] /= s; + } + + /** + * Return the dot product of a vector. + * @return the dot product of a vector. + **/ + final double dotProduct() { + double s = 0.0; + for (int i = 0; i < NDIM; i++) + s += data[i] * data[i]; + return s; + } + + final double absolute() { + double tmp = 0.0; + for (int i = 0; i < NDIM; i++) + tmp += data[i] * data[i]; + return sqrt(tmp); + } + + final double distance(Vec3 v) { + double tmp = 0.0; + for (int i = 0; i < NDIM; i++) + tmp += ((data[i] - v.data[i]) * (data[i] - v.data[i])); + return sqrt(tmp); + } + + final void crossProduct(Vec3 u, Vec3 w) { + data[0] = u.data[1] * w.data[2] - u.data[2] * w.data[1]; + data[1] = u.data[2] * w.data[0] - u.data[0] * w.data[2]; + data[2] = u.data[0] * w.data[1] - u.data[1] * w.data[0]; + } + + final void incrementalAdd(Vec3 u) { + for (int i = 0; i < NDIM; i++) + data[i] += u.data[i]; + } + + final void incrementalSub(Vec3 u) { + for (int i = 0; i < NDIM; i++) + data[i] -= u.data[i]; + } + + final void incrementalMultScalar(double s) { + for (int i = 0; i < NDIM; i++) + data[i] *= s; + } + + final void incrementalDivScalar(double s) { + for (int i = 0; i < NDIM; i++) + data[i] /= s; + } + + final void addScalar(Vec3 u, double s) { + for (int i = 0; i < NDIM; i++) + data[i] = u.data[i] + s; + } + + public char* toCString() { + // this is not generic code at all + char* result = cast(char*)malloc(100); + if (result == null) exit(1); + sprintf(result, "%.17f %.17f %.17f ", data[0], data[1], data[2]); + return result; + } +} + + +/// A class that represents the common fields of a cell or body data structure. +abstract class Node { + // mass of the node + double mass; + + // Position of the node + Vec3 pos; + + // highest bit of int coord + const int IMAX = 1073741824; + + // potential softening parameter + const double EPS = 0.05; + + /// Construct an empty node + protected this() { + mass = 0.0; + pos = new Vec3(); + } + + abstract Node loadTree(Body p, Vec3 xpic, int l, Tree root); + + abstract double hackcofm(); + + abstract HG walkSubTree(double dsq, HG hg); + + static final int oldSubindex(Vec3 ic, int l) { + int i = 0; + for (int k = 0; k < Vec3.NDIM; k++) + if ((cast(int)ic.value(k) & l) != 0) + i += Cell.NSUB >> (k + 1); + return i; + } + + public char* toCString() { + char* result = cast(char*)malloc(130); + if (result == null) exit(1); + char* pos_str = pos.toCString(); + sprintf(result, "%f : %s", mass, pos_str); + free(pos_str); + return result; + } + + /// Compute a single body-body or body-cell interaction + final HG gravSub(HG hg) { + Vec3 dr = new Vec3(); + dr.subtraction(pos, hg.pos0); + + double drsq = dr.dotProduct() + (EPS * EPS); + double drabs = sqrt(drsq); + + double phii = mass / drabs; + hg.phi0 -= phii; + double mor3 = phii / drsq; + dr.multScalar(mor3); + hg.acc0.addition(dr); + return hg; + } + + /** + * A sub class which is used to compute and save information during the + * gravity computation phase. + **/ + static protected final class HG { + /// Body to skip in force evaluation + Body pskip; + + /// Point at which to evaluate field + Vec3 pos0; + + /// Computed potential at pos0 + + double phi0; + + /// computed acceleration at pos0 + Vec3 acc0; + + /** + * Create a HG object. + * @param b the body object + * @param p a vector that represents the body + **/ + this(Body b, Vec3 p) { + pskip = b; + pos0 = p.clone(); + phi0 = 0.0; + acc0 = new Vec3(); + } + } +} + + +/// A class used to representing particles in the N-body simulation. +final class Body : Node { + Vec3 vel, acc, newAcc; + double phi = 0.0; + private Body next, procNext; + + /// Create an empty body. + this() { + vel = new Vec3(); + acc = new Vec3(); + newAcc = new Vec3(); + } + + /** + * Set the next body in the list. + * @param n the body + **/ + final void setNext(Body n) { + next = n; + } + + /** + * Get the next body in the list. + * @return the next body + **/ + final Body getNext() { + return next; + } + + /** + * Set the next body in the list. + * @param n the body + **/ + final void setProcNext(Body n) { + procNext = n; + } + + /** + * Get the next body in the list. + * @return the next body + **/ + final Body getProcNext() { + return procNext; + } + + /** + * Enlarge cubical "box", salvaging existing tree structure. + * @param tree the root of the tree. + * @param nsteps the current time step + **/ + final void expandBox(Tree tree, int nsteps) { + Vec3 rmid = new Vec3(); + + bool inbox = icTest(tree); + while (!inbox) { + double rsize = tree.rsize; + rmid.addScalar(tree.rmin, 0.5 * rsize); + + for (int k = 0; k < Vec3.NDIM; k++) + if (pos.value(k) < rmid.value(k)) { + double rmin = tree.rmin.value(k); + tree.rmin.value(k, rmin - rsize); + } + + tree.rsize = 2.0 * rsize; + if (tree.root !is null) { + Vec3 ic = tree.intcoord(rmid); + if (ic is null) + throw new Exception("Value is out of bounds"); + int k = oldSubindex(ic, IMAX >> 1); + Cell newt = new Cell(); + newt.subp[k] = tree.root; + tree.root = newt; + inbox = icTest(tree); + } + } + } + + /** + * Check the bounds of the body and return true if it isn't in the + * correct bounds. + **/ + final bool icTest(Tree tree) { + double pos0 = pos.value(0); + double pos1 = pos.value(1); + double pos2 = pos.value(2); + + // by default, it is in bounds + bool result = true; + + double xsc = (pos0 - tree.rmin.value(0)) / tree.rsize; + if (!(0.0 < xsc && xsc < 1.0)) + result = false; + + xsc = (pos1 - tree.rmin.value(1)) / tree.rsize; + if (!(0.0 < xsc && xsc < 1.0)) + result = false; + + xsc = (pos2 - tree.rmin.value(2)) / tree.rsize; + if (!(0.0 < xsc && xsc < 1.0)) + result = false; + + return result; + } + + /** + * Descend Tree and insert particle. We're at a body so we need to + * create a cell and attach this body to the cell. + * @param p the body to insert + * @param xpic + * @param l + * @param tree the root of the data structure + * @return the subtree with the new body inserted + **/ + final Node loadTree(Body p, Vec3 xpic, int l, Tree tree) { + // create a Cell + Cell retval = new Cell(); + int si = subindex(tree, l); + // attach this Body node to the cell + retval.subp[si] = this; + + // move down one level + si = oldSubindex(xpic, l); + Node rt = retval.subp[si]; + if (rt !is null) + retval.subp[si] = rt.loadTree(p, xpic, l >> 1, tree); + else + retval.subp[si] = p; + return retval; + } + + /** + * Descend tree finding center of mass coordinates + * @return the mass of this node + **/ + final double hackcofm() { + return mass; + } + + /// iteration of the bodies + int opApply(int delegate(ref Body) dg) { + int result; + Body current = this; + while (current !is null) { + result = dg(current); + current = current.next; + if (result) + break; + } + return result; + } + + /// inverse iteration of the bodies + int opApplyReverse(int delegate(ref Body) dg) { + int result; + Body current = this; + while (current !is null) { + result = dg(current); + current = current.procNext; + if (result) + break; + } + return result; + } + + + /** + * Return an enumeration of the bodies + * @return an enumeration of the bodies + **/ + final Enumeration elements() { + // a local class that implements the enumerator + static final class Enumerate : Enumeration { + private Body current; + public this(Body outer) { + //this.current = Body.this; + this.current = outer; + } + public bool hasMoreElements() { + return current !is null; + } + public Object nextElement() { + Object retval = current; + current = current.next; + return retval; + } + } + + return new Enumerate(this); + } + + final Enumeration elementsRev() { + // a local class that implements the enumerator + static class Enumerate : Enumeration { + private Body current; + public this(Body outer) { + //this.current = Body.this; + this.current = outer; + } + public bool hasMoreElements() { + return current !is null; + } + public Object nextElement() { + Object retval = current; + current = current.procNext; + return retval; + } + } + + return new Enumerate(this); + } + + /** + * Determine which subcell to select. + * Combination of intcoord and oldSubindex. + * @param t the root of the tree + **/ + final int subindex(Tree tree, int l) { + Vec3 xp = new Vec3(); + + double xsc = (pos.value(0) - tree.rmin.value(0)) / tree.rsize; + xp.value(0, floor(IMAX * xsc)); + + xsc = (pos.value(1) - tree.rmin.value(1)) / tree.rsize; + xp.value(1, floor(IMAX * xsc)); + + xsc = (pos.value(2) - tree.rmin.value(2)) / tree.rsize; + xp.value(2, floor(IMAX * xsc)); + + int i = 0; + for (int k = 0; k < Vec3.NDIM; k++) + if ((cast(int)xp.value(k) & l) != 0) + i += Cell.NSUB >> (k + 1); + return i; + } + + /** + * Evaluate gravitational field on the body. + * The original olden version calls a routine named "walkscan", + * but we use the same name that is in the Barnes code. + **/ + final void hackGravity(double rsize, Node root) { + Vec3 pos0 = pos.clone(); + HG hg = new HG(this, pos); + hg = root.walkSubTree(rsize * rsize, hg); + phi = hg.phi0; + newAcc = hg.acc0; + } + + /// Recursively walk the tree to do hackwalk calculation + final HG walkSubTree(double dsq, HG hg) { + if (this != hg.pskip) + hg = gravSub(hg); + return hg; + } + + /** + * Return a string represenation of a body. + * @return a string represenation of a body. + **/ + public char* toCString() { + char* result = cast(char*)malloc(140); + if (result == null) exit(1); + char* super_str = super.toCString(); + sprintf(result, "Body %s", super_str); + free(super_str); + return result; + } +} + + + +/// A class used to represent internal nodes in the tree +final class Cell : Node { + // subcells per cell + const int NSUB = 1 << Vec3.NDIM; + + /** + * The children of this cell node. Each entry may contain either + * another cell or a body. + **/ + Node[] subp; + Cell next; + + this() { + subp.length = NSUB; + } + + /** + * Descend Tree and insert particle. We're at a cell so + * we need to move down the tree. + * @param p the body to insert into the tree + * @param xpic + * @param l + * @param tree the root of the tree + * @return the subtree with the new body inserted + **/ + Node loadTree(Body p, Vec3 xpic, int l, Tree tree) { + // move down one level + int si = oldSubindex(xpic, l); + Node rt = subp[si]; + if (rt !is null) + subp[si] = rt.loadTree(p, xpic, l >> 1, tree); + else + subp[si] = p; + return this; + } + + /** + * Descend tree finding center of mass coordinates + * @return the mass of this node + **/ + double hackcofm() { + double mq = 0.0; + Vec3 tmpPos = new Vec3(); + Vec3 tmpv = new Vec3(); + for (int i = 0; i < NSUB; i++) { + Node r = subp[i]; + if (r !is null) { + double mr = r.hackcofm(); + mq = mr + mq; + tmpv.multScalar(r.pos, mr); + tmpPos.addition(tmpv); + } + } + mass = mq; + pos = tmpPos; + pos.divScalar(mass); + + return mq; + } + + /// Recursively walk the tree to do hackwalk calculation + HG walkSubTree(double dsq, HG hg) { + if (subdivp(dsq, hg)) { + for (int k = 0; k < Cell.NSUB; k++) { + Node r = subp[k]; + if (r !is null) + hg = r.walkSubTree(dsq / 4.0, hg); + } + } else { + hg = gravSub(hg); + } + return hg; + } + + /** + * Decide if the cell is too close to accept as a single term. + * @return true if the cell is too close. + **/ + bool subdivp(double dsq, HG hg) { + Vec3 dr = new Vec3(); + dr.subtraction(pos, hg.pos0); + double drsq = dr.dotProduct(); + + // in the original olden version drsp is multiplied by 1.0 + return drsq < dsq; + } + + /** + * Return a string represenation of a cell. + * @return a string represenation of a cell. + **/ + public char* toCString() { + char* result = cast(char*)malloc(140); + if (result == null) exit(1); + char* super_str = super.toCString(); + sprintf(result, "Cell %s", super_str); + free(super_str); + return result; + } +} + + +/** +* A class that represents the root of the data structure used +* to represent the N-bodies in the Barnes-Hut algorithm. +**/ +final class Tree { + Vec3 rmin; + double rsize; + + /// A reference to the root node. + Node root; + + /// The complete list of bodies that have been created. + private Body bodyTab; + + /// The complete list of bodies that have been created - in reverse. + private Body bodyTabRev; + + /// Construct the root of the data structure that represents the N-bodies. + this() { + rmin = new Vec3(); + rsize = -2.0 * -2.0; + root = null; + bodyTab = null; + bodyTabRev = null; + rmin.value(0, -2.0); + rmin.value(1, -2.0); + rmin.value(2, -2.0); + } + + /** + * Return an enumeration of the bodies. + * @return an enumeration of the bodies. + **/ + final Enumeration bodies() { + return bodyTab.elements(); + } + + /** + * Return an enumeration of the bodies - in reverse. + * @return an enumeration of the bodies - in reverse. + **/ + final Enumeration bodiesRev() { + return bodyTabRev.elementsRev(); + } + + /** + * Create the testdata used in the benchmark. + * @param nbody the number of bodies to create + **/ + final void createTestData(int nbody) { + Vec3 cmr = new Vec3(); + Vec3 cmv = new Vec3(); + Body head = new Body(); + Body prev = head; + + double rsc = 3.0 * PI / 16.0; + double vsc = sqrt(1.0 / rsc); + int seed = 123; + Random rnd = Random(seed); + + for (int i = 0; i < nbody; i++) { + Body p = new Body(); + + prev.setNext(p); + prev = p; + p.mass = 1.0 / cast(double)nbody; + + double t1 = rnd.uniform(0.0, 0.999); + t1 = pow(t1, (-2.0 / 3.0)) - 1.0; + double r = 1.0 / sqrt(t1); + + double coeff = 4.0; + for (int k = 0; k < Vec3.NDIM; k++) { + r = rnd.uniform(0.0, 0.999); + p.pos.value(k, coeff * r); + } + + cmr.addition(p.pos); + + double x, y; + do { + x = rnd.uniform(0.0, 1.0); + y = rnd.uniform(0.0, 0.1); + } while (y > (x * x * pow(1.0 - x * x, 3.5))); + + double v = sqrt(2.0) * x / pow(1 + r * r, 0.25); + + double rad = vsc * v; + double rsq; + do { + for (int k = 0; k < Vec3.NDIM; k++) + p.vel.value(k, rnd.uniform(-1.0, 1.0)); + rsq = p.vel.dotProduct(); + } while (rsq > 1.0); + double rsc1 = rad / sqrt(rsq); + p.vel.multScalar(rsc1); + cmv.addition(p.vel); + } + + // mark end of list + prev.setNext(null); + + // toss the dummy node at the beginning and set a reference to the first element + bodyTab = head.getNext(); + + cmr.divScalar(cast(double)nbody); + cmv.divScalar(cast(double)nbody); + + prev = null; + + for (Enumeration e = bodyTab.elements(); e.hasMoreElements();) { + Body b = cast(Body)e.nextElement(); + b.pos.subtraction(cmr); + b.vel.subtraction(cmv); + b.setProcNext(prev); + prev = b; + } + + // set the reference to the last element + bodyTabRev = prev; + } + + + /** + * Advance the N-body system one time-step. + * @param nstep the current time step + **/ + void stepSystem(int nstep) { + // free the tree + root = null; + + makeTree(nstep); + + // compute the gravity for all the particles + for (Enumeration e = bodyTabRev.elementsRev(); e.hasMoreElements();) { + Body b = cast(Body)e.nextElement(); + b.hackGravity(rsize, root); + } + vp(bodyTabRev, nstep); + } + + /** + * Initialize the tree structure for hack force calculation. + * @param nsteps the current time step + **/ + private void makeTree(int nstep) { + for (Enumeration e = bodiesRev(); e.hasMoreElements();) { + Body q = cast(Body)e.nextElement(); + if (q.mass != 0.0) { + q.expandBox(this, nstep); + Vec3 xqic = intcoord(q.pos); + if (root is null) { + root = q; + } else { + root = root.loadTree(q, xqic, Node.IMAX >> 1, this); + } + } + } + root.hackcofm(); + } + + /** + * Compute integerized coordinates. + * @return the coordinates or null if rp is out of bounds + **/ + final Vec3 intcoord(Vec3 vp) { + Vec3 xp = new Vec3(); + + double xsc = (vp.value(0) - rmin.value(0)) / rsize; + if (0.0 <= xsc && xsc < 1.0) + xp.value(0, floor(Node.IMAX * xsc)); + else + return null; + + xsc = (vp.value(1) - rmin.value(1)) / rsize; + if (0.0 <= xsc && xsc < 1.0) + xp.value(1, floor(Node.IMAX * xsc)); + else + return null; + + xsc = (vp.value(2) - rmin.value(2)) / rsize; + if (0.0 <= xsc && xsc < 1.0) + xp.value(2, floor(Node.IMAX * xsc)); + else + return null; + + return xp; + } + + static final private void vp(Body p, int nstep) { + Vec3 dacc = new Vec3(); + Vec3 dvel = new Vec3(); + double dthf = 0.5 * BH.DTIME; + + for (Enumeration e = p.elementsRev(); e.hasMoreElements();) { + Body b = cast(Body)e.nextElement(); + Vec3 acc1 = b.newAcc.clone(); + if (nstep > 0) { + dacc.subtraction(acc1, b.acc); + dvel.multScalar(dacc, dthf); + dvel.addition(b.vel); + b.vel = dvel.clone(); + } + + b.acc = acc1.clone(); + dvel.multScalar(b.acc, dthf); + + Vec3 vel1 = b.vel.clone(); + vel1.addition(dvel); + Vec3 dpos = vel1.clone(); + dpos.multScalar(BH.DTIME); + dpos.addition(b.pos); + b.pos = dpos.clone(); + vel1.addition(dvel); + b.vel = vel1.clone(); + } + } +} + + +final class BH { + /// The user specified number of bodies to create. + private static int nbody = 0; + + /// The maximum number of time steps to take in the simulation + private static int nsteps = 10; + + /// Should we print information messsages + private static bool printMsgs = false; + + /// Should we print detailed results + private static bool printResults = false; + + const double DTIME = 0.0125; + private const double TSTOP = 2.0; + + public static final void main(char[][] args) { + parseCmdLine(args); + + if (printMsgs) + printf("nbody = %d\n", nbody); + + auto start0 = myclock(); + Tree root = new Tree(); + root.createTestData(nbody); + auto end0 = myclock(); + if (printMsgs) + printf("Bodies created\n"); + + auto start1 = myclock(); + double tnow = 0.0; + int i = 0; + while ((tnow < TSTOP + 0.1 * DTIME) && (i < nsteps)) { + root.stepSystem(i++); + tnow += DTIME; + } + auto end1 = myclock(); + + if (printResults) { + int j = 0; + for (Enumeration e = root.bodies(); e.hasMoreElements();) { + Body b = cast(Body)e.nextElement(); + char* str_ptr = b.pos.toCString(); + printf("body %d: %s\n", j++, str_ptr); + free(str_ptr); + } + } + + if (printMsgs) { + printf("Build time %.2f\n", end0 - start0); + printf("Compute Time %.2f\n", end1 - start1); + printf("Total Time %.2f\n", end1 - start0); + printf("Done!\n"); + } + } + + private static final void parseCmdLine(char[][] args) { + int i = 1; + while (i < args.length && args[i][0] == '-') { + char[] arg = args[i++]; + + // check for options that require arguments + if (arg == "-b") { + if (i < args.length) + nbody = toInt(args[i++]); + else + throw new Exception("-l requires the number of levels"); + } else if (arg == "-s") { + if (i < args.length) + nsteps = 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 (nbody == 0) + usage(); + } + + /// The usage routine which describes the program options. + private static final void usage() { + fprintf(stderr, "usage: BH -b [-s ] [-p] [-m] [-h]\n"); + fprintf(stderr, " -b the number of bodies\n"); + fprintf(stderr, " -s the max. number of time steps (default=10)\n"); + fprintf(stderr, " -p (print detailed results)\n"); + fprintf(stderr, " -m (print information messages\n"); + exit(0); + } +} + + +void main(char[][] args) { + BH.main(args); +} diff --git a/micro/big_arrays.d b/micro/bigarr.d similarity index 100% rename from micro/big_arrays.d rename to micro/bigarr.d diff --git a/micro/bisort.d b/micro/bisort.d new file mode 100644 index 0000000..6425535 --- /dev/null +++ b/micro/bisort.d @@ -0,0 +1,366 @@ +/* +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 BiSort + * 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 n 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 [-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); +} diff --git a/micro/conalloc.d b/micro/conalloc.d new file mode 100644 index 0000000..fba0f20 --- /dev/null +++ b/micro/conalloc.d @@ -0,0 +1,50 @@ +// 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); +} + diff --git a/micro/concpu.d b/micro/concpu.d new file mode 100644 index 0000000..6a123b2 --- /dev/null +++ b/micro/concpu.d @@ -0,0 +1,50 @@ +// 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); +} + diff --git a/micro/em3d.d b/micro/em3d.d new file mode 100644 index 0000000..349f6a7 --- /dev/null +++ b/micro/em3d.d @@ -0,0 +1,468 @@ +/** + * Java implementation of the em3d 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. + * + *

+ * 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. + * + * + * 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 -d [-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); +} diff --git a/micro/multicore.d b/micro/mcore.d similarity index 100% rename from micro/multicore.d rename to micro/mcore.d diff --git a/micro/rnd_data_2.d b/micro/rnd_data_2.d deleted file mode 100644 index 545d588..0000000 --- a/micro/rnd_data_2.d +++ /dev/null @@ -1,40 +0,0 @@ -// Written by Kevin Bealer -// Found at http://www.digitalmars.com/webnews/newsgroups.php?art_group=digitalmars.D.announce&article_id=6978 -// Sightly modified by Leandro Lucarella -// (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; -} - diff --git a/micro/rnd_data.d b/micro/rnddata.d similarity index 52% rename from micro/rnd_data.d rename to micro/rnddata.d index aa71355..c0d2294 100644 --- a/micro/rnd_data.d +++ b/micro/rnddata.d @@ -5,20 +5,28 @@ 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... } diff --git a/micro/sbtree.d b/micro/sbtree.d new file mode 100644 index 0000000..3f8c904 --- /dev/null +++ b/micro/sbtree.d @@ -0,0 +1,63 @@ +// Written by Andrey Khropov +// 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; + } +} + diff --git a/micro/shootout_binarytrees.d b/micro/shootout_binarytrees.d deleted file mode 100644 index 7663458..0000000 --- a/micro/shootout_binarytrees.d +++ /dev/null @@ -1,72 +0,0 @@ -// Written by Andrey Khropov -// 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; -} - diff --git a/micro/tree.d b/micro/tree.d deleted file mode 100644 index a5179a3..0000000 --- a/micro/tree.d +++ /dev/null @@ -1,46 +0,0 @@ -// Written by Piotr Modzelewski -// 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; - } - -} - diff --git a/micro/tsp.d b/micro/tsp.d new file mode 100644 index 0000000..26f0376 --- /dev/null +++ b/micro/tsp.d @@ -0,0 +1,650 @@ +/** +A D implementation of the "tsp" Olden benchmark, the traveling +salesman problem. + +R. Karp, "Probabilistic analysis of partitioning algorithms for the +traveling-salesman problem in the plane." Mathematics of Operations Research +2(3):209-224, August 1977. + +Converted to D and optimized by leonardo maffi, V.1.0, Oct 29 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_tsp.zip +*/ + +version (Tango) { + import tango.stdc.stdio: printf, fprintf, stderr; + import tango.stdc.time: CLOCKS_PER_SEC, clock; + import tango.math.Math: sqrt, log; + + import Integer = tango.text.convert.Integer; + alias Integer.parse toInt; +} else { + import std.c.stdio: printf, fprintf, stderr; + import std.c.time: CLOCKS_PER_SEC, clock; + import std.math: sqrt, log; + + 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. + +Very simple portable random number generator 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 = 48271; + 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; + } +} + + + +/** +* A class that represents a node in a binary tree. Each node represents +* a city in the TSP benchmark. +**/ +final class Tree { + /** + * The number of nodes (cities) in this subtree + **/ + private int sz; + + /** + * The coordinates that this node represents + **/ + private double x, y; + + /** + * Left and right child of tree + **/ + private Tree left, right; + + /** + * The next pointer in a linked list of nodes in the subtree. The list + * contains the order of the cities to visit. + **/ + private Tree next; + + /** + * The previous pointer in a linked list of nodes in the subtree. The list + * contains the order of the cities to visit. + **/ + private Tree prev; + + static Random rnd; + + // used by the random number generator + private const double M_E2 = 7.3890560989306502274; + private const double M_E3 = 20.08553692318766774179; + private const double M_E6 = 403.42879349273512264299; + private const double M_E12 = 162754.79141900392083592475; + + public static void initRnd(int seed) { + rnd = new Random(seed); + } + + /** + * Construct a Tree node (a city) with the specified size + * @param size the number of nodes in the (sub)tree + * @param x the x coordinate of the city + * @param y the y coordinate of the city + * @param left the left subtree + * @param right the right subtree + **/ + this(int size, double x, double y, Tree l, Tree r) { + sz = size; + this.x = x; + this.y = y; + left = l; + right = r; + next = null; + prev = null; + } + + /** + * Find Euclidean distance between this node and the specified node. + * @param b the specified node + * @return the Euclidean distance between two tree nodes. + **/ + double distance(Tree b) { + return sqrt((x - b.x) * (x - b.x) + (y - b.y) * (y - b.y)); + } + + /** + * Create a list of tree nodes. Requires root to be the tail of the list. + * Only fills in next field, not prev. + * @return the linked list of nodes + **/ + Tree makeList() { + Tree myleft, myright, tleft, tright; + Tree retval = this; + + // head of left list + if (left !is null) + myleft = left.makeList(); + else + myleft = null; + + // head of right list + if (right !is null) + myright = right.makeList(); + else + myright = null; + + if (myright !is null) { + retval = myright; + right.next = this; + } + + if (myleft !is null) { + retval = myleft; + if (myright !is null) + left.next = myright; + else + left.next = this; + } + next = null; + + return retval; + } + + /** + * Reverse the linked list. Assumes that there is a dummy "prev" + * node at the beginning. + **/ + void reverse() { + Tree prev = this.prev; + prev.next = null; + this.prev = null; + Tree back = this; + Tree tmp = this; + + // reverse the list for the other nodes + Tree next; + for (Tree t = this.next; t !is null; back = t, t = next) { + next = t.next; + t.next = back; + back.prev = t; + } + + // reverse the list for this node + tmp.next = prev; + prev.prev = tmp; + } + + /** + * Use closest-point heuristic from Cormen, Leiserson, and Rivest. + * @return a + **/ + Tree conquer() { + // create the list of nodes + Tree t = makeList(); + + // Create initial cycle + Tree cycle = t; + t = t.next; + cycle.next = cycle; + cycle.prev = cycle; + + // Loop over remaining points + Tree donext; + for ( ; t !is null; t = donext) { + donext = t.next; /* value won't be around later */ + Tree min = cycle; + double mindist = t.distance(cycle); + for (Tree tmp = cycle.next; tmp != cycle; tmp = tmp.next) { + double test = tmp.distance(t); + if (test < mindist) { + mindist = test; + min = tmp; + } + } + + Tree next = min.next; + Tree prev = min.prev; + + double mintonext = min.distance(next); + double mintoprev = min.distance(prev); + double ttonext = t.distance(next); + double ttoprev = t.distance(prev); + + if ((ttoprev - mintoprev) < (ttonext - mintonext)) { + // insert between min and prev + prev.next = t; + t.next = min; + t.prev = prev; + min.prev = t; + } else { + next.prev = t; + t.next = next; + min.next = t; + t.prev = min; + } + } + + return cycle; + } + + /** + * Merge two cycles as per Karp. + * @param a a subtree with one cycle + * @param b a subtree with the other cycle + **/ + Tree merge(Tree a, Tree b) { + // Compute location for first cycle + Tree min = a; + double mindist = distance(a); + Tree tmp = a; + for (a = a.next; a != tmp; a = a.next) { + double test = distance(a); + if (test < mindist) { + mindist = test; + min = a; + } + } + + Tree next = min.next; + Tree prev = min.prev; + double mintonext = min.distance(next); + double mintoprev = min.distance(prev); + double ttonext = distance(next); + double ttoprev = distance(prev); + + Tree p1, n1; + double tton1, ttop1; + if ((ttoprev - mintoprev) < (ttonext - mintonext)) { + // would insert between min and prev + p1 = prev; + n1 = min; + tton1 = mindist; + ttop1 = ttoprev; + } else { + // would insert between min and next + p1 = min; + n1 = next; + ttop1 = mindist; + tton1 = ttonext; + } + + // Compute location for second cycle + min = b; + mindist = distance(b); + tmp = b; + for (b = b.next; b != tmp; b = b.next) { + double test = distance(b); + if (test < mindist) { + mindist = test; + min = b; + } + } + + next = min.next; + prev = min.prev; + mintonext = min.distance(next); + mintoprev = min.distance(prev); + ttonext = this.distance(next); + ttoprev = this.distance(prev); + + Tree p2, n2; + double tton2, ttop2; + if ((ttoprev - mintoprev) < (ttonext - mintonext)) { + // would insert between min and prev + p2 = prev; + n2 = min; + tton2 = mindist; + ttop2 = ttoprev; + } else { + // would insert between min andn ext + p2 = min; + n2 = next; + ttop2 = mindist; + tton2 = ttonext; + } + + // Now we have 4 choices to complete: + // 1:t,p1 t,p2 n1,n2 + // 2:t,p1 t,n2 n1,p2 + // 3:t,n1 t,p2 p1,n2 + // 4:t,n1 t,n2 p1,p2 + double n1ton2 = n1.distance(n2); + double n1top2 = n1.distance(p2); + double p1ton2 = p1.distance(n2); + double p1top2 = p1.distance(p2); + + mindist = ttop1 + ttop2 + n1ton2; + int choice = 1; + + double test = ttop1 + tton2 + n1top2; + if (test < mindist) { + choice = 2; + mindist = test; + } + + test = tton1 + ttop2 + p1ton2; + if (test < mindist) { + choice = 3; + mindist = test; + } + + test = tton1 + tton2 + p1top2; + if (test < mindist) + choice = 4; + + switch (choice) { + case 1: + // 1:p1,this this,p2 n2,n1 -- reverse 2! + n2.reverse(); + p1.next = this; + this.prev = p1; + this.next = p2; + p2.prev = this; + n2.next = n1; + n1.prev = n2; + break; + + case 2: + // 2:p1,this this,n2 p2,n1 -- OK + p1.next = this; + this.prev = p1; + this.next = n2; + n2.prev = this; + p2.next = n1; + n1.prev = p2; + break; + + case 3: + // 3:p2,this this,n1 p1,n2 -- OK + p2.next = this; + this.prev = p2; + this.next = n1; + n1.prev = this; + p1.next = n2; + n2.prev = p1; + break; + + default: // case 4: + // 4:n1,this this,n2 p2,p1 -- reverse 1! + n1.reverse(); + n1.next = this; + this.prev = n1; + this.next = n2; + n2.prev = this; + p2.next = p1; + p1.prev = p2; + break; + } + + return this; + } // end merge + + /** + * Compute TSP for the tree t. Use conquer for problems <= sz + * @param sz the cutoff point for using conquer vs. merge + **/ + Tree tsp(int sz) { + if (this.sz <= sz) + return conquer(); + Tree leftval = left.tsp(sz); + Tree rightval = right.tsp(sz); + return merge(leftval, rightval); + } + + /** + * Print the list of cities to visit from the current node. The + * list is the result of computing the TSP problem. + * The list for the root node (city) should contain every other node + * (city). + **/ + void printVisitOrder() { + printf("x = %.15f y = %.15f\n", x, y); + for (Tree tmp = next; tmp != this; tmp = tmp.next) + printf("x = %.15f y = %.15f\n", tmp.x, tmp.y); + } + + /** + * Computes the total length of the current tour. + **/ + double tourLength() { + double total = 0.0; + Tree precedent = next; + Tree current = next.next; + if (current == this) + return total; + + do { + total += current.distance(precedent); + precedent = current; + current = current.next; + } while (precedent != this); + + total += current.distance(this); + return total; + } + + + // static methods =============================================== + + /** + * Return an estimate of median of n values distributed in [min, max) + * @param min the minimum value + * @param max the maximum value + * @param n + * @return an estimate of median of n values distributed in [min, max) + **/ + private static double median(double min, double max, int n) { + // get random value in [0.0, 1.0) + double t = rnd.nextDouble(); + + double retval; + if (t > 0.5) + retval = log(1.0 - (2.0 * (M_E12 - 1) * (t - 0.5) / M_E12)) / 12.0; + else + retval = -log(1.0 - (2.0 * (M_E12 - 1) * t / M_E12)) / 12.0; + + // We now have something distributed on (-1.0, 1.0) + retval = (retval + 1.0) * (max - min) / 2.0; + return retval + min; + } + + /** + * Get double uniformly distributed over [min,max) + * @return double uniformily distributed over [min,max) + **/ + private static double uniform(double min, double max) { + // get random value between [0.0, 1.0) + double retval = rnd.nextDouble(); + retval = retval * (max - min); + return retval + min; + } + + /** + * Builds a 2D tree of n nodes in specified range with dir as primary + * axis (false for x, true for y) + * + * @param n the size of the subtree to create + * @param dir the primary axis + * @param min_x the minimum x coordinate + * @param max_x the maximum x coordinate + * @param min_y the minimum y coordinate + * @param max_y the maximum y coordinate + * @return a reference to the root of the subtree + **/ + public static Tree buildTree(int n, bool dir, double min_x, + double max_x, double min_y, double max_y) { + if (n == 0) + return null; + + Tree left, right; + double x, y; + if (dir) { + dir = !dir; + double med = median(min_x, max_x, n); + left = buildTree(n/2, dir, min_x, med, min_y, max_y); + right = buildTree(n/2, dir, med, max_x, min_y, max_y); + x = med; + y = uniform(min_y, max_y); + } else { + dir = !dir; + double med = median(min_y,max_y,n); + left = buildTree(n/2, dir, min_x, max_x, min_y, med); + right = buildTree(n/2, dir, min_x, max_x, med, max_y); + y = med; + x = uniform(min_x, max_x); + } + + return new Tree(n, x, y, left, right); + } +} + + +public class TSP { + /** + * Number of cities in the problem. + **/ + private static int cities; + + /** + * Set to true if the result should be printed + **/ + private static bool printResult = false; + + /** + * Set to true to print informative messages + **/ + private static bool printMsgs = false; + + /** + * The main routine which creates a tree and traverses it. + * @param args the arguments to the program + **/ + public static void main(char[] args[]) { + if (!parseCmdLine(args)) + return; + + if (printMsgs) + printf("Building tree of size: %d\n", nextPow2(cities+1) - 1); + + Tree.initRnd(1); + + auto t0 = myclock(); + Tree t = Tree.buildTree(cities, false, 0.0, 1.0, 0.0, 1.0); + + auto t1 = myclock(); + t.tsp(150); + auto t2 = myclock(); + + if (printMsgs) + printf("Total tour length: %.15f\n", t.tourLength()); + auto t3 = myclock(); + + if (printResult) + // if the user specifies, print the final result + t.printVisitOrder(); + + if (printMsgs) { + printf("Tsp build time: %.2f\n", t1 - t0); + printf("Tsp computing time: %.2f\n", t2 - t1); + printf("Tsp total time: %.2f\n", t3 - t0); + } + } + + private static /*unsigned*/ int nextPow2(/*unsigned*/ int x) { + if (x < 0) + throw new Exception("x must be >= 0"); + x = x - 1; + x = x | (x >> 1); + x = x | (x >> 2); + x = x | (x >> 4); + x = x | (x >> 8); + x = x | (x >> 16); + return x + 1; + } + + /** + * Parse the command line options. + * @param args the command line options. + **/ + private static final bool parseCmdLine(char[] args[]) { + int i = 1; + + while (i < args.length && args[i][0] == '-') { + char[] arg = args[i++]; + + if (arg == "-c") { + if (i < args.length) + cities = toInt(args[i++]); + else + throw new Exception("-c requires the size of tree"); + if (cities < 1) + throw new Exception("Number of cities must be > 0"); + } else if (arg == "-p") { + printResult = true; + } else if (arg == "-m") { + printMsgs = true; + } else if (arg == "-h") { + return usage(); + } + } + + if (cities == 0) + return usage(); + + return true; + } + + /** + * The usage routine which describes the program options. + **/ + private static final bool usage() { + fprintf(stderr, "usage: tsp_d -c [-p] [-m] [-h]\n"); + fprintf(stderr, " -c number of cities (rounds up to the next power of 2 minus 1)\n"); + fprintf(stderr, " -p (print the final result)\n"); + fprintf(stderr, " -m (print informative messages)\n"); + fprintf(stderr, " -h (print this message)\n"); + return false; + } +} + + +void main(char[][] args) { + TSP.main(args); +} -- 2.43.0