games/galaxy: parallelize gravitational force calculations

Once the Barnes-Hut tree is constructed, the gravitational
force calculations can be done in parallel by dividing the
bodies up between a number of procs.

sys/man/1/galaxy

@@ -102,23 +102,11 @@ Exit the simulator.
 Certain aspects of the galaxy simulator are controlled by
 the following options:
 .TP
-.BI -t " throttle"
-Causes the process that calculates forces to relinquish
-the processor for
-.I throttle
-milliseconds after each calculation.
-.TP
 .BI -G " gravity"
 Sets the gravitational constant to
 .I gravity.
 The default value is 1.
 .TP
-.BI -ε " softening"
-Sets the
-.I softening
-factor to prevent gravitational singularities during
-collisions or near-collisions. The default value is 500.
-.TP
 .BI -f " file"
 Reads the galaxy file
 .I file
@@ -127,6 +115,25 @@ Reads the galaxy file
 .TP
 .B -i
 Reads a galaxy file from standard input.
+.TP
+.BI -p " procs"
+Specifies the number of extra processes to use in order
+to calculate the gravitational force on each body in
+parallel.
+The default value is
+.BR $NPROC-1 .
+.TP
+.BI -t " throttle"
+Causes the process that calculates forces to relinquish
+the processor for
+.I throttle
+milliseconds after each calculation.
+.TP
+.BI -ε " softening"
+Sets the
+.I softening
+factor to prevent gravitational singularities during
+collisions or near-collisions. The default value is 500.
 .SS Mkgalaxy
 .I Mkgalaxy
 is a utility to create galaxies for simulation.

sys/src/games/galaxy/body.c

@@ -11,7 +11,7 @@ glxyinit(void)
 	glxy.a = calloc(5, sizeof(Body));
 	if(glxy.a == nil)
 		sysfatal("could not allocate glxy: %r");
-	glxy.l = 0;
+	glxy.nb = 0;
 	glxy.max = 5;
 }
 
@@ -20,13 +20,13 @@ body(void)
 {
 	Body *b;
 
-	if(glxy.l == glxy.max) {
+	if(glxy.nb == glxy.max) {
 		glxy.max *= 2;
 		glxy.a = realloc(glxy.a, sizeof(Body) * glxy.max);
 		if(glxy.a == nil)
 			sysfatal("could not realloc glxy: %r");
 	}
-	b = glxy.a + glxy.l++;
+	b = glxy.a + glxy.nb++;
 	*b = ZB;
 	return b;
 }
@@ -64,11 +64,11 @@ center(void)
 	Vector gc, gcv;
 	double mass;
 
-	if(glxy.l == 0)
+	if(glxy.nb == 0)
 		return (Vector){0, 0};
 
 	gc.x = gc.y = gcv.x = gcv.y = mass = 0;
-	for(b = glxy.a; b < glxy.a+glxy.l; b++) {
+	for(b = glxy.a; b < glxy.a+glxy.nb; b++) {
 		gc.x += b->x * b->mass;
 		gc.y += b->y * b->mass;
 		gcv.x += b->v.x * b->mass;
@@ -79,7 +79,7 @@ center(void)
 	gc.y /= mass;
 	gcv.x /= mass;
 	gcv.y /= mass;
-	for(b = glxy.a; b < glxy.a+glxy.l; b++) {
+	for(b = glxy.a; b < glxy.a+glxy.nb; b++) {
 		b->x -= gc.x;
 		b->y -= gc.y;
 		b->v.x -= gcv.x;
@@ -145,7 +145,7 @@ readglxy(int fd)
 	int cmd, len;
 	Body *b;
 
-	glxy.l = 0;
+	glxy.nb = 0;
 	Binit(&bin, fd, OREAD);
 	for(;;) {
 		line = Brdline(&bin, ' ');
@@ -212,7 +212,7 @@ writeglxy(int fd)
 	Bprint(&bout, "DT %g\n", dt);
 	Bprint(&bout, "GRAV %g\n", G);
 
-	for(b = glxy.a; b < glxy.a + glxy.l; b++)
+	for(b = glxy.a; b < glxy.a + glxy.nb; b++)
 		Bprint(&bout, "%B\n", b);
 
 	Bterm(&bout);

sys/src/games/galaxy/galaxy.c

@@ -46,7 +46,6 @@ Cursor pausecursor={
 };
 
 enum {
-	STK = 8192,
 	DOBODY = 0,
 	SPEED,
 	GRAV,
@@ -68,7 +67,7 @@ double
 	LIM = 10,
 	dt²;
 char *file;
-int showv, showa, throttle, paused;
+int showv, showa, paused;
 
 char *menustr[] = {
 	[DOBODY]	"new body",
@@ -149,7 +148,7 @@ drawstats(void)
 	Point p;
 	static char buf[1024];
 
-	snprint(buf, sizeof(buf), "Number of bodies: %d", glxy.l);
+	snprint(buf, sizeof(buf), "Number of bodies: %d", glxy.nb);
 	p = addpt(screen->r.min, (Point){5, 3});
 	string(screen, p, display->white, ZP, font, buf);
 
@@ -170,7 +169,7 @@ drawglxy(void)
 	int s;
 
 	draw(screen, screen->r, display->black, 0, ZP);
-	for(b = glxy.a; b < glxy.a + glxy.l; b++) {
+	for(b = glxy.a; b < glxy.a + glxy.nb; b++) {
 		pos = topoint(b->Vector);
 		s = b->size/scale;
 		fillellipse(screen, pos, s, s, b->col, ZP);
@@ -529,57 +528,6 @@ kbdthread(void*)
 	}
 }
 
-/* verlet barnes-hut */
-void
-simulate(void*)
-{
-	Body *b;
-	double f;
-
-	threadsetname("simulate");
-
-	for(;;) {
-		qlock(&glxy);
-
-		if(throttle)
-			sleep(throttle);
-
-		drawglxy();
-
-Again:
-		space.t = EMPTY;
-		quads.l = 0;
-		STATS(quaddepth = 0;)
-		for(b = glxy.a; b < glxy.a + glxy.l; b++) {
-			if(quadins(b, LIM) == -1) {
-				growquads();
-				goto Again;
-			}
-		}
-
-		STATS(avgcalcs = 0;)
-		for(b = glxy.a; b < glxy.a + glxy.l; b++) {
-			b->a.x = b->newa.x;
-			b->a.y = b->newa.y;
-			b->newa.x = b->newa.y = 0;
-			STATS(calcs = 0;)
-			quadcalc(b, space, LIM);
-			STATS(avgcalcs += calcs;)
-		}
-		STATS(avgcalcs /= glxy.l;)
-
-		for(b = glxy.a; b < glxy.a + glxy.l; b++) {
-			b->x += dt*b->v.x + dt²*b->a.x/2;
-			b->y += dt*b->v.y + dt²*b->a.y/2;
-			b->v.x += dt*(b->a.x + b->newa.x)/2;
-			b->v.y += dt*(b->a.y + b->newa.y)/2;
-			CHECKLIM(b, f);
-		}
-
-		qunlock(&glxy);
-	}
-}
-
 Vector
 tovector(Point p)
 {
@@ -600,6 +548,7 @@ usage(void)
 void
 threadmain(int argc, char **argv)
 {
+	char* nproc;
 	int doload;
 
 	doload = 0;
@@ -607,6 +556,9 @@ threadmain(int argc, char **argv)
 	default:
 		usage();
 		break;
+	case 'p':
+		extraproc = strtol(EARGF(usage()), nil, 0);
+		break;
 	case 't':
 		throttle = strtol(EARGF(usage()), nil, 0);
 		break;
@@ -637,6 +589,16 @@ threadmain(int argc, char **argv)
 			sysfatal("threadmain: could not open file: %r");
 	}
 
+	if(extraproc < 0) {
+		nproc = getenv("NPROC");
+		if(nproc == nil)
+			extraproc = 0;
+		else
+			extraproc = strtol(nproc, nil, 10) - 1;
+		if(extraproc < 0)
+			extraproc = 0;
+	}
+
 	if(initdraw(nil, nil, "Galaxy") < 0)
 		sysfatal("initdraw failed: %r");
 	if(mc = initmouse(nil, screen), mc == nil)

sys/src/games/galaxy/galaxy.h

@@ -33,7 +33,7 @@ struct Quad {
 struct {
 	QLock;
 	Body *a;
-	int l, max;
+	int nb, max;
 } glxy;
 
 struct {
@@ -49,11 +49,6 @@ enum {
 	BODY,
 };
 
-Point orig;
-double G, θ, scale, ε, LIM, dt, dt²;
-Body ZB;
-QB space;
-
 Image *randcol(void);
 Point topoint(Vector);
 Vector tovector(Point);
@@ -65,16 +60,24 @@ void glxyinit(void);
 int Bfmt(Fmt*);
 void readglxy(int);
 void writeglxy(int);
+void drawglxy(void);
+
+void simulate(void*);
 
 void quadcalc(Body*, QB, double);
 int quadins(Body*, double);
 void growquads(void);
 void quadsinit(void);
 
-int stats;
+int stats, quaddepth, calcs, extraproc, throttle;
-int quaddepth;
+double G, θ, scale, ε, LIM, dt, dt², avgcalcs;
-int calcs;
+Point orig;
-double avgcalcs;
+QB space;
+Body ZB;
+
+enum {
+	STK = 8192,
+};
 
 #define STATS(e) if(stats) {e}
 

sys/src/games/galaxy/mkfile

@@ -2,7 +2,7 @@
 
 TARG=galaxy mkgalaxy
 BIN=/$objtype/bin/games
-OGALAXY=galaxy.$O quad.$O body.$O
+OGALAXY=galaxy.$O quad.$O body.$O simulate.$O
 OMKGALAXY=mkgalaxy.$O body.$O
 
 </sys/src/cmd/mkmany

sys/src/games/galaxy/mkgalaxy.c

@@ -152,7 +152,7 @@ main(int argc, char **argv)
 	if(argc != 1)
 		usage();
 
-	new = glxy.l;
+	new = glxy.nb;
 	lim = strtod(*argv, nil);
 	mkbodies(lim);
 
@@ -160,7 +160,7 @@ main(int argc, char **argv)
 	for(b = glxy.a; b < glxy.a + new; b++)
 		Bprint(&bout, "%B\n", b);
 
-	for(b = glxy.a+new; b < glxy.a+glxy.l; b++) {
+	for(b = glxy.a+new; b < glxy.a+glxy.nb; b++) {
 		b->x += o.x;
 		b->y += o.y;
 		Bprint(&bout, "%B\n", b);

sys/src/games/galaxy/simulate.c

@@ -0,0 +1,142 @@
+#include <u.h>
+#include <libc.h>
+#include <draw.h>
+#include <thread.h>
+#include "galaxy.h"
+
+int extraproc = -1, throttle;
+
+static QLock* golock;
+static Rendez* gorend;
+static int* go;
+
+static QLock runninglock;
+static Rendez runningrend;
+static int running;
+
+static void
+calcproc(void *v)
+{
+	Body *b, *e;
+	int nbody;
+	int pid;
+
+	pid = (uintptr)v;
+
+	for(;;) {
+		qlock(golock+pid);
+		while(!go[pid])
+			rsleep(gorend+pid);
+		go[pid] = 0;
+		qunlock(golock+pid);
+
+		nbody = glxy.nb / (extraproc+1);
+		b = glxy.a + nbody * pid;
+		e = b + nbody;
+		while(b < e) {
+			b->a.x = b->newa.x;
+			b->a.y = b->newa.y;
+			b->newa.x = b->newa.y = 0;
+			quadcalc(b++, space, LIM);
+		}	
+
+		qlock(&runninglock);
+		if(--running == 0)
+			rwakeup(&runningrend);
+		qunlock(&runninglock);
+	}
+}
+
+static void
+startprocs(void)
+{
+	int pid;
+
+	golock = calloc(extraproc, sizeof(*golock));
+	if(golock == nil)
+		sysfatal("Could not create go locks: %r\n");
+
+	gorend = calloc(extraproc, sizeof(*gorend));
+	if(gorend == nil)
+		sysfatal("Could not create go rendez: %r\n");
+
+	go = calloc(extraproc, sizeof(*go));
+	if(go == nil)
+		sysfatal("Could not create go flags: %r\n");
+
+	for(pid = 0; pid < extraproc; pid++)
+		gorend[pid].l = golock+pid;
+	runningrend.l = &runninglock;
+
+	for(pid = 0; pid < extraproc; pid++) 
+		proccreate(calcproc, (void*)pid, STK);
+}
+
+/* verlet barnes-hut */
+void
+simulate(void*)
+{
+	Body *b, *s;
+	int nbody, pid;
+	double f;
+
+	threadsetname("simulate");
+
+	startprocs();
+
+	for(;;) {
+		qlock(&glxy);
+
+		if(throttle)
+			sleep(throttle);
+
+		drawglxy();
+
+Again:
+		space.t = EMPTY;
+		quads.l = 0;
+		STATS(quaddepth = 0;)
+		for(b = glxy.a; b < glxy.a + glxy.nb; b++) {
+			if(quadins(b, LIM) == -1) {
+				growquads();
+				goto Again;
+			}
+		}
+
+		running = extraproc;
+		for(pid = 0; pid < extraproc; pid++) {
+			qlock(golock+pid);
+			go[pid] = 1;
+			rwakeup(gorend+pid);
+			qunlock(golock+pid);
+		}
+
+		STATS(avgcalcs = 0;)
+		nbody = glxy.nb / (extraproc+1);
+		s = glxy.a + nbody * (extraproc);
+		for(b = s; b < glxy.a + glxy.nb; b++) {
+			b->a.x = b->newa.x;
+			b->a.y = b->newa.y;
+			b->newa.x = b->newa.y = 0;
+			STATS(calcs = 0;)
+			quadcalc(b, space, LIM);
+			STATS(avgcalcs += calcs;)
+		}
+		STATS(avgcalcs /= glxy.nb;)
+
+		qlock(&runninglock);
+		while(running > 0)
+			rsleep(&runningrend);
+		qunlock(&runninglock);
+
+		for(b = glxy.a; b < glxy.a + glxy.nb; b++) {
+			b->x += dt*b->v.x + dt²*b->a.x/2;
+			b->y += dt*b->v.y + dt²*b->a.y/2;
+			b->v.x += dt*(b->a.x + b->newa.x)/2;
+			b->v.y += dt*(b->a.y + b->newa.y)/2;
+			CHECKLIM(b, f);
+		}
+
+		qunlock(&glxy);
+	}
+}