The n-body problem is a classic example used in parallel programming. The straight forward solution takes O(N^2) time complexity, whereas the Barnes-Hut reduce it to O(NlogN) based on the assumption that, a group of bodies that is far enough from a body can be clustered in to a BIG body.

Based on that assumption, we can divide the space into sub 4 equal rectangles (2D) or 8 cubes (3D). Then we build a quad tree (2D) or oct tree (3D) whose root is the bounding box of all n bodies and divide each node into 4(2D)/8(3D) until each leaf contains no more than 1 body. Each node will need to store its bounding box, its mass (total mass of the children) and center of gravity, which is calculated by the formula:

COM = (C1*M1+C2*M2+C3*M3+C4*M4)/M (for quad tree)
COM = (C1*M1+…+C8*M8)/M (for quad tree)
where Ci is the position of Child ith, Mi is the respective mass and M is the total mass of all children

To compute the total force acting on body X, we just need to traverse down the tree, at each node Y, we calculate the distant R of the node to X. if R >= d*theta, where d is the largest size of the bounding box enclosing Y, theta is a constant (in my code I choose theta = 1.2), we treat Y as a body and calculate its force on X. Otherwise, we go down to Y’s children.

I have written 2 versions. One is Quadtree, and the other uses OpenMp. The quadtre version is no doubt faster but openmp is much shorter and easy to implement. And it considerably faster than naive O(N^2) version.

To compile the sourcode under Linux, you need to add “-lX11 -lm” as the program will rely on XWindows to draw
To compile OpenMp version, add -openmp

QuadTree version

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <sys/time.h>
#include <vector>

#define N 1000000
#define G 6.673e-11
#define T_MAX 100
#define dt 1
#define THETA 1.2
#define  X_RESN	800       /* x resolution */
#define  Y_RESN	800       /* y resolution */

using namespace std;


typedef struct complextype 
{
    float real, imag;
} Compl;

struct  QTreeNode {
	float x0,y0,x1,y1,mass;
	float COM[2]; // center of mass
	bool is_leaf;
	vector<int> id;	
	struct QTreeNode** children;
};



Window        win;  /* initialization for a window */
GC black_gc; // to draw the body
GC white_gc; // to clear the body
GC gc;
Display       *display;
unsigned int  screen; /* which screen */

float m[N]; // mass 
float f[N][2],v[N][2],p[N][2]; // forces, velocity, position
int n,t_max;

void init_display () {
    unsigned 
    int           width, height,  /* window size */
                  x, y, /* window position */
                  border_width, /*border width in pixels */
                  display_width, display_height;
                  /* size of screen */
                  char *window_name = "Mandelbrot Set", 
                  *display_name = NULL;
    unsigned
    long          valuemask = 0;
    XGCValues     values;
    XSizeHints    size_hints;
    Pixmap        bitmap;
    XPoint        points[800];
    FILE          *fp, *fopen ();
    char          str[100];

	XSetWindowAttributes attr[1];

    /* Mandlebrot variables */
    int           i, j, k;
    
    /* connect to Xserver */

    if ( (display = XOpenDisplay (display_name)) == NULL ) {
	  fprintf (stderr, 
                  "Cannot connect to X server %s\n",
                  XDisplayName (display_name) );
	  exit (-1);
    }
	
    /* get screen size */

    screen = DefaultScreen(display);
    display_width = DisplayWidth(display, screen);
    display_height = DisplayHeight(display, screen);

    /* set window size */

    width = X_RESN;
    height = Y_RESN;

    /* set window position */

    x = 0;
    y = 0;

    /* create opaque window */
	border_width = 4;
    win = XCreateSimpleWindow(display, 
                  RootWindow(display, screen),
                  x, y, width, height, border_width, 
                  BlackPixel(display, screen), 
                  WhitePixel(display, screen));

    size_hints.flags = USPosition|USSize;
    size_hints.x = x;
    size_hints.y = y;
    size_hints.width = width;
    size_hints.height = height;
    size_hints.min_width = 300;
    size_hints.min_height = 300;

    XSetNormalHints(display, win, &size_hints);
    XStoreName(display, win, window_name);
    /* create graphics context */
	/*
	values.foreground = WhitePixel (display, screen);
    white_gc = XCreateGC(display, win, valuemask, &values);
	values.foreground = BlackPixel (display, screen);
	black_gc = XCreateGC (display, win, valuemask, &values);
	*/
	gc = XCreateGC (display, win, valuemask, &values);
    XSetBackground(display, gc, 
                  WhitePixel (display, screen));
    XSetForeground(display, gc, 
                  BlackPixel (display, screen));
    XSetLineAttributes(display, gc, 1, 
                  LineSolid, CapRound, JoinRound);

    attr[0].backing_store = Always;
    attr[0].backing_planes = 1;
    attr[0].backing_pixel = BlackPixel(display, screen);

    XChangeWindowAttributes(display, win, 
                  CWBackingStore | CWBackingPlanes | CWBackingPixel, 
                  attr);

    XMapWindow(display, win);
    XSync(display, 0);

}


float abs (float x) {
	return (x>=0)?x:-x;
}

void build_qtree (struct QTreeNode* aNode) {
	if (aNode->id.size () < 2) {
		aNode->is_leaf = true;
		return;
	}
	aNode->children = new QTreeNode*[4];
	float x0,y0,x1,y1;
	x0 = aNode->x0; 
	x1 = aNode->x1;
	y0 = aNode->y0;
	y1 = aNode->y1;
	for (int i=0;i<4;++i) { 
		aNode->children[i] = new QTreeNode;
		aNode->is_leaf = false;
		if (i==0||i==2) {
			aNode->children[i]->x0 = x0;
			aNode->children[i]->x1 = (x0+x1)*0.5;
		}
		else {
			aNode->children[i]->x0 = (x0+x1)*0.5;
			aNode->children[i]->x1 = x1;
		}
		if (i<2) {
			aNode->children[i]->y0 = y0;
			aNode->children[i]->y1 = (y0+y1)*0.5;
		}
		else {
			aNode->children[i]->y0 = (y0+y1)*0.5;
			aNode->children[i]->y1 = y1;
		}
	}
	for (int i = 0;i<aNode->id.size ();++i) {
		// decide which child the body belongs to 
		int pid = aNode->id[i];
		float x = p[pid][0];
		float y = p[pid][1];
		for (int j=0;j<4;++j) {
			if (x>=aNode->children[j]->x0 && x<=aNode->children[j]->x1 && y>=aNode->children[j]->y0 && y<=aNode->children[j]->y1) {
				aNode->children[j]->id.push_back (pid);
				break;
			}
		}
	}
	for (int i=0;i<4;++i)
		build_qtree (aNode->children[i]);
}

void compute_mass (struct QTreeNode* aNode) {
	if (aNode->id.size () == 0) return;
	aNode->mass = 0;
	aNode->COM[0] = aNode->COM[1] = 0;
	if (aNode->is_leaf) {
		aNode->mass = m[aNode->id[0]];
		aNode->COM[0] = p[aNode->id[0]][0];
		aNode->COM[1] = p[aNode->id[0]][1];
	}
	else {
		for (int i=0;i<4;++i) {
			struct QTreeNode* child = aNode->children[i];
			compute_mass (child);
			aNode->mass+=child->mass;
			aNode->COM[0]+=child->mass*child->COM[0];
			aNode->COM[1]+=child->mass*child->COM[1];
		}
		aNode->COM[0]/= aNode->mass;
		aNode->COM[1]/= aNode->mass;
	}
}

void force (const float& x1, const float& y1, const float& x2, const float& y2, const float& m1, const float& m2, float& Fx, float& Fy) {
	float r = sqrt ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2));
	if (r>0) {
		r*=r*r;
		Fx +=G*m1*m2*(x2-x1) / r;
		Fy +=G*m1*m2*(y2-y1) / r;
	}
}
void cal_force (struct QTreeNode* aNode, const float& mm, const float& x, const float& y, float& Fx, float& Fy) {
	if (aNode->is_leaf) {
		if (aNode->id.size () == 0) return;
		int pid = aNode->id[0];
		force (x,y,p[pid][0],p[pid][1],mm,m[pid], Fx, Fy);
		return;
	} 
	float r = sqrt (pow (aNode->COM[0]-x,2) + pow (aNode->COM[1]-y,2));	
	float d = max (abs (aNode->x1-aNode->x0), abs (aNode->y1-aNode->y0));

	// if only cluster of bodies is far enough
	if (r>d*THETA) {
		force (x,y,aNode->COM[0],aNode->COM[1],mm,aNode->mass,Fx,Fy);
		return;
	}
	for (int i=0;i<4;++i) 
		cal_force (aNode->children[i], mm, x,y,Fx,Fy);
}

void compute_force (struct QTreeNode* aNode) {
	for (int i=0;i<n;++i) {
		cal_force (aNode, m[i], p[i][0], p[i][1], f[i][0],f[i][1]);	
	}
}
/*
void traverse (struct QTreeNode* aNode) {
	if (aNode->id.size ()==0) cout<<"nil"<<endl;
	else if (aNode->is_leaf) cout<<"L: "<<aNode->id[0]<<endl;
	else {
		for (int i = 0;i<aNode->id.size();++i) 
			cout<<aNode->id[i]<<" ";
		cout<<endl;
		for (int i=0;i<4;++i)
			traverse (aNode->children[i]);
	}
}
*/

main (int argc, char* argv[])
{
	init_display ();     	 
	FILE* fin = fopen ("nbody.txt","r");
    /* Calculate and draw points */
    Compl         z, c;
    float         lengthsq, temp;
	fscanf (fin,"%d %d",&n,&t_max);

		
	for (int i = 0;i<n;++i) {
		m[i] = 1;
		f[i][0] = f[i][1] = v[i][0] = v[i][1] = 0;
		fscanf (fin, "%f %f",&p[i][0],&p[i][1]);
	}

	float p_new[n][2],v_new[n][2];
	cout<<n<<endl;

	struct timeval tv1,tv2;
	struct timezone tz1,tz2;
	gettimeofday (&tv1,&tz1);
	for (int t=0;t<t_max;++t) {
		// build quad tree
		struct QTreeNode* root = new QTreeNode;
		root->x0 = root->y0 = 0;
		root->x1 = root->y1 = 1;
		root->is_leaf = false;
		for (int i=0;i<n;++i) 
			root->id.push_back (i);
		build_qtree (root);
		//traverse (root);
		compute_mass (root);
		compute_force (root);
		
		for (int i=0;i<n;++i) {
			v_new[i][0] = v[i][0] + f[i][0]*dt/m[i];
			v_new[i][1] = v[i][1] + f[i][1]*dt/m[i];
			p_new[i][0] = p[i][0]+v_new[i][0]*dt/m[i];
			p_new[i][1] = p[i][1]+v_new[i][1]*dt/m[i];
		}

		for (int i=0;i<n;++i) {
			int x,y;
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			//XSetForeground (display, gc, WhitePixel (display,screen));
			//XDrawPoint (display, win, gc, x,y);
			v[i][0] = v_new[i][0];
			v[i][1] = v_new[i][1];
			p[i][0] = p_new[i][0];
			p[i][1] = p_new[i][1];
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			//cout<<x<<" "<<y<<endl;
			//XSetForeground (display, gc, BlackPixel (display,screen));
			//XDrawPoint (display, win, gc, x, y);
		}
		//sleep (1);
		
	}
	gettimeofday (&tv2,&tz2);
	printf ("Eslaped time = %ld\n", (tv2.tv_sec-tv1.tv_sec)*1000000 + tv2.tv_usec-tv1.tv_usec);
    XFlush(display);

    /* Program Finished */

}

OpenPM version

#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#include <sys/time.h>
#include <omp.h>

#define N 10   				// number of bodies
#define G 6.673e-11
#define T_MAX 100
#define N_THREADS 7			// default number of threads
#define dt 1
#define  X_RESN	800       /* x resolution */
#define  Y_RESN	800       /* y resolution */

using namespace std;


typedef struct complextype 
{
    float real, imag;
} Compl;


Window        win;  /* initialization for a window */
GC black_gc; // to draw the body
GC white_gc; // to clear the body
GC gc;
Display       *display;
unsigned int  screen; /* which screen */

void init_display () {
    unsigned 
    int           width, height,  /* window size */
                  x, y, /* window position */
                  border_width, /*border width in pixels */
                  display_width, display_height;
                  /* size of screen */
                  char *window_name = "Mandelbrot Set", 
                  *display_name = NULL;
    unsigned
    long          valuemask = 0;
    XGCValues     values;
    XSizeHints    size_hints;
    Pixmap        bitmap;
    XPoint        points[800];
    FILE          *fp, *fopen ();
    char          str[100];

	XSetWindowAttributes attr[1];

    /* Mandlebrot variables */
    int           i, j, k;
    
    /* connect to Xserver */

    if ( (display = XOpenDisplay (display_name)) == NULL ) {
	  fprintf (stderr, 
                  "Cannot connect to X server %s\n",
                  XDisplayName (display_name) );
	  exit (-1);
    }
	
    /* get screen size */

    screen = DefaultScreen(display);
    display_width = DisplayWidth(display, screen);
    display_height = DisplayHeight(display, screen);

    /* set window size */

    width = X_RESN;
    height = Y_RESN;

    /* set window position */

    x = 0;
    y = 0;

    /* create opaque window */

    border_width = 4;
    win = XCreateSimpleWindow(display, 
                  RootWindow(display, screen),
                  x, y, width, height, border_width, 
                  BlackPixel(display, screen), 
                  WhitePixel(display, screen));

    size_hints.flags = USPosition|USSize;
    size_hints.x = x;
    size_hints.y = y;
    size_hints.width = width;
    size_hints.height = height;
    size_hints.min_width = 300;
    size_hints.min_height = 300;

    XSetNormalHints(display, win, &size_hints);
    XStoreName(display, win, window_name);
    /* create graphics context */
	/*
	values.foreground = WhitePixel (display, screen);
    white_gc = XCreateGC(display, win, valuemask, &values);
	values.foreground = BlackPixel (display, screen);
	black_gc = XCreateGC (display, win, valuemask, &values);
	*/
	gc = XCreateGC (display, win, valuemask, &values);
    XSetBackground(display, gc, 
                  WhitePixel (display, screen));
    XSetForeground(display, gc, 
                  BlackPixel (display, screen));
    XSetLineAttributes(display, gc, 1, 
                  LineSolid, CapRound, JoinRound);

	//cout<<"here"<<endl;
    attr[0].backing_store = Always;
    attr[0].backing_planes = 1;
    attr[0].backing_pixel = BlackPixel(display, screen);

    XChangeWindowAttributes(display, win, 
                  CWBackingStore | CWBackingPlanes | CWBackingPixel, 
                  attr);

    XMapWindow(display, win);
    XSync(display, 0);
}

main (int argc, char* argv[])
{
	init_display ();     	 
    /* Calculate and draw points */
    Compl         z, c;
    float         lengthsq, temp;
	int n,t_max;
	FILE* fin = fopen ("nbody.txt","r");
	fscanf (fin, "%d %d",&n,&t_max);
	int n_threads = (argc==2)?atoi(argv[3]):N_THREADS;
	double m[n];
	double f[n][2],v[n][2],p[n][2];
	
	#pragma omp parallel for 
	for (int i = 0;i<n;++i) {
		m[i] = 1;
		f[i][0] = f[i][1] = v[i][0] = v[i][1] = 0;
		fscanf (fin,"%lf %lf",&p[i][0],&p[i][1]);
	}
	double p_new[n][2],v_new[n][2];
	omp_set_num_threads (n_threads);
	struct timeval tv1,tv2;
	struct timezone tz1,tz2;
	gettimeofday (&tv1,&tz1);
	for (int t=0;t<t_max;++t) {
		//cout<<t<<endl;
		#pragma omp parallel for 
		for (int i=0;i<n;++i) {
			for (int j=0;j<n;++j) {
				double r = sqrt (pow (p[j][0]-p[i][0],2) + pow(p[j][1]-p[i][1],2));
				r*=r*r;
				if (r>0) {
					f[i][0]+=G*m[i]*m[j]*(p[j][0]-p[i][0]) / r;
					f[i][1]+=G*m[i]*m[j]*(p[j][1]-p[i][1]) / r;
				}
			}
			v_new[i][0] = v[i][0] + f[i][0]*dt/m[i];
			v_new[i][1] = v[i][1] + f[i][1]*dt/m[i];
			p_new[i][0] = p[i][0]+v_new[i][0]*dt/m[i];
			p_new[i][1] = p[i][1]+v_new[i][1]*dt/m[i];

		}
		#pragma omp barrier

		#pragma omp parallel for 
		for (int i=0;i<n;++i) {
			int x,y;
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			#pragma omp critical 
			{
				XSetForeground (display, gc, WhitePixel (display,screen));
				XDrawPoint (display, win, gc, x,y);
			}
			v[i][0] = v_new[i][0];
			v[i][1] = v_new[i][1];
			p[i][0] = p_new[i][0];
			p[i][1] = p_new[i][1];
			x = floor (p[i][0]*X_RESN);
			y = floor (p[i][1]*Y_RESN);
			#pragma omp critical 
			{
				XSetForeground (display, gc, BlackPixel (display,screen));
				XDrawPoint (display, win, gc, x, y);
			}
		}
		//sleep(1);
	}
	gettimeofday (&tv2,&tz2);
	printf ("Eslaped time = %ld\n", (tv2.tv_sec-tv1.tv_sec)*1000000 + tv2.tv_usec-tv1.tv_usec);
    XFlush(display);

    /* Program Finished */

}

Advertisements