barnes.cpp

一个用OpenMP实现的并行Barnes Hut算法。有schedule和chunk size的设置功能。运行环境：vs2005

#include "stdafx.h"
#include <iostream>
#include <string.h>
#include <math.h>
#include <omp.h>
#include <fstream>
#include <time.h>
using namespace std;

#define NoOfBody
#ifdef NoOfBody
#define InitNo 1000
#define InitBody 1000
char* FileName="1000bodies.txt";
#else
#define InitNo 20
#define InitBody 20
char* FileName="20bodies.txt";
#endif

#define MaxX 200
#define MaxY 200
#define MaxMass 10e9
#define TimeStep 1
#define G 6.673e-11
#define Theta 1.0
#define Scale 4
#define StepCount 500
//#define StepCount 1

const int num_of_threads = 4; //number of threads created
int ForceCalculations=0;

typedef struct{
	float mass;
	float dx;  //velocity in the x dimension
	float dy;  //velocity in the y dimension
	float x;   //x dimension offset from (0,0) at lower left
	float y;   //y dimension offset from (0,0) at lower left
} body;

typedef struct{
	float x;
	float y;
} point;

typedef struct{
	float x;
	float y;
} force;

body** bodies;

//class TreeNode
class TreeNode{
public:
	float minX;
	float minY;
	float maxX; 
	float maxY;
	float edgeLength;
	float containedMass;
	int childNo; //the no of children
	int bodyNo; //the no of contained bodies
	point centerOfMass;
	body** containedBodies;
	TreeNode** children;

	TreeNode(float, float, float, float, body***, int);
	~TreeNode();
	void addBodies(body***, int);
	void buildQuadTree();
    bool inBounds(float, float);
	void computeCenterOfMass();
	bool notEmpty();
	point centerPoint();	
};

TreeNode::TreeNode(float x1, float y1, float x2, float y2, body*** bodies, int count){
//The 5th parameter is the length of list bodies
	this->minX=x1;
	this->minY=y1;
	this->maxX=x2;
	this->maxY=y2;
	this->edgeLength=x2-x1;
	this->bodyNo=0;
	this->childNo=0;
	this->containedMass=0;
	this->addBodies(bodies,count);
	buildQuadTree();
	if(this->notEmpty()){
		this->computeCenterOfMass();
	}
}

bool TreeNode::inBounds(float x, float y){
	return (this->minX<=x)&&(x<this->maxX)&&(this->minY<=y)&&(y<this->maxY);
}

void TreeNode::addBodies(body*** bodies, int count){
	
	body* body_tmp=new body [InitNo];
	int i,j=0;
    for(i=0;i<count;i++){
		if(inBounds((*bodies)[i]->x,(*bodies)[i]->y)){
			body_tmp[j++]=(*(*bodies)[i]);
			this->bodyNo++;
		}
	}
	if(this->bodyNo==0){
		delete body_tmp;
		this->containedBodies=NULL;
		return;
	}
	this->containedBodies=new body* [this->bodyNo];
	for(i=0;i<this->bodyNo;i++){
		this->containedBodies[i]=new body;
		*(this->containedBodies[i])=body_tmp[i];
	}
	delete body_tmp;
}

point TreeNode::centerPoint(){
	point p;
	p.x=this->minX+(this->maxX-this->minX)/2.0;
	p.y=this->minY+(this->maxY-this->minY)/2.0;
	return p;
}

bool TreeNode::notEmpty(){
	return (!(this->bodyNo==0));
}

void TreeNode::buildQuadTree(){
	if(this->bodyNo<=1){
		this->childNo=0;this->children=NULL;
		return;
	}
	TreeNode** children_tmp=new TreeNode*[4];
	float cx=this->centerPoint().x;
	float cy=this->centerPoint().y;
	children_tmp[0]=new TreeNode(this->minX, this->minY, cx, cy, &(this->containedBodies), this->bodyNo);
	children_tmp[1]=new TreeNode(cx, this->minY, this->maxX, cy, &(this->containedBodies), this->bodyNo);
	children_tmp[2]=new TreeNode(this->minX, cy, cx, this->maxY, &(this->containedBodies), this->bodyNo);
	children_tmp[3]=new TreeNode(cx, cy, this->maxX, this->maxY, &(this->containedBodies), this->bodyNo);
	this->childNo=0;
	for(int i=0;i<4;i++){
		if(children_tmp[i]->bodyNo>0)
			this->childNo++;
	}
	this->children=new TreeNode*[childNo];
	if(this->childNo!=0){
		int i=0;
		for(int j=0;j<4;j++){
			if (children_tmp[j]->bodyNo>0){
				children[i]=children_tmp[j];
				i++;}//if
			else{
				delete children_tmp[j];
			}
		}//for
	}//if	
	else
		this->children=NULL;
	delete [] children_tmp;
}

void TreeNode::computeCenterOfMass(){
	if(this->bodyNo==1){
		this->containedMass=this->containedBodies[0]->mass;
		(this->centerOfMass).x=this->containedBodies[0]->x;
		(this->centerOfMass).y=this->containedBodies[0]->y;
	}
	else{
		float xSum=0.0;
		float ySum=0.0;				
		for (int i= 0; i<this->childNo; i++){
			xSum+=(this->children[i]->containedMass)*((this->children[i]->centerOfMass).x);
			ySum+=(this->children[i]->containedMass)*((this->children[i]->centerOfMass).y);
			this->containedMass+=(this->children[i]->containedMass);
		}
		(this->centerOfMass).x=xSum/(this->containedMass);
		(this->centerOfMass).y=ySum/(this->containedMass);
	}
}

TreeNode::~TreeNode(){
	   for(int i=0;i<this->childNo;i++){
			delete this->children[i];
		}
	   delete [] children;
	   children=NULL;
	   for(int i=0;i<this->bodyNo;i++){
		   delete this->containedBodies[i];
	   }
       delete [] this->containedBodies;
	   containedBodies=NULL;
}

float ComputeRange(point a, point b){
	float xd=b.x-a.x;
	float yd=b.y-a.y;
	return sqrt(xd*xd+yd*yd);
}

bool PointInList(point p, body** list, int count)
{
	for (int i=0; i<count; i++){
		if ((p.x==list[i]->x)&&(p.y==list[i]->y))
			return 1;
	}
	return 0;
}

force ComputeForce(body* b, TreeNode* tree){
	point position;
	force f;
	position.x=b->x;
	position.y=b->y;
	float range=ComputeRange(position,tree->centerOfMass);
	if(range<0.2){
		f.x=0;
		f.y=0;
		return f;
	}
	if(range>=tree->edgeLength/Theta||tree->bodyNo==1){
		float rcubed=range*range*range;
		float invariant=G*(b->mass)*(tree->containedMass)/rcubed;
		float fx=invariant*((tree->centerOfMass).x-position.x);
		float fy=invariant*((tree->centerOfMass).y-position.y);
		ForceCalculations++;
		f.x=fx;
		f.y=fy;
		return f;
	}
	else{
		f.x=0;
		f.y=0;
		force childForce;
		for(int i=0;i<tree->childNo;i++){
			childForce=ComputeForce(b,(tree->children)[i]);
			f.x+=childForce.x;
			f.y+=childForce.y;
		}//for
		return f;
	}//else
}

int ComputeForcesAndUpdateBodies(body*** list, TreeNode* tree, int count){
	float vx, vy, x, y;
	int countNew=0;
	
	//newlist=new body*[count];
    int i=0;
	force *f=new force [count];
	body** newlist=new body* [count];

/**********************************
//The following works in sequential
	for(i=0;i<count;i++){
			f[i]=ComputeForce((*list)[i], tree);
			vx=(*list)[i]->dx+f[i].x*TimeStep/(*list)[i]->mass;
			vy=(*list)[i]->dy+f[i].y*TimeStep/(*list)[i]->mass;
			x=(*list)[i]->x+vx*TimeStep;
			y=(*list)[i]->y+vy*TimeStep;
			if((x>=0)&&(x<=tree->maxX)&&(y>=0)&&(y<=tree->maxY)){
			newlist[countNew]=new body;
			newlist[countNew]->mass=(*list)[i]->mass;
			newlist[countNew]->dx=vx;
			newlist[countNew]->dy=vy;
			newlist[countNew]->x=x;
			newlist[countNew]->y=y;
			countNew++;
			}//if
		}//for
//***********************************/
	
	omp_set_num_threads(num_of_threads);
	#pragma omp parallel private(i) shared(list,tree) num_threads(num_of_threads)		
	{
		#pragma omp for schedule(dynamic,1)		//specify for loop schedule (dynamic,1)
		for(i=0;i<count;i++){
			f[i]=ComputeForce((*list)[i], tree);
		}
	}
	for(i=0;i<count;i++){
			vx=(*list)[i]->dx+f[i].x*TimeStep/(*list)[i]->mass;
			vy=(*list)[i]->dy+f[i].y*TimeStep/(*list)[i]->mass;
			x=(*list)[i]->x+vx*TimeStep;
			y=(*list)[i]->y+vy*TimeStep;
			if((x>=0)&&(x<=tree->maxX)&&(y>=0)&&(y<=tree->maxY)){
			newlist[countNew]=new body;
			newlist[countNew]->mass=(*list)[i]->mass;
			newlist[countNew]->dx=vx;
			newlist[countNew]->dy=vy;
			newlist[countNew]->x=x;
			newlist[countNew]->y=y;
			countNew++;
			}//if
		}//for

	for(i=0;i<count;i++){
		delete (*list)[i];
	}
	delete [] (*list);
    *list=newlist; 
	delete f;
	return countNew;
}


void traverse(TreeNode* root){
	if(root->bodyNo==0){
		return;
	}
	if(root->bodyNo==1){
		cout<<"bodyNo= childNo="<<root->bodyNo<<" "<<root->childNo<<endl;
		cout<<root->containedBodies[0]->x<<" "<<root->containedBodies[0]->y<<endl;
		return;
	}//
	cout<<"bodyNo= childNo="<<root->bodyNo<<" "<<root->childNo<<endl;
	for(int i=0;i<root->bodyNo;i++){
		cout<<root->containedBodies[i]->x<<" "<<root->containedBodies[i]->y<<endl;
	}

	for(int i=0;i<root->childNo;i++){
		traverse(root->children[i]);
	}
}

void main(){
    int count=0,i;
	//char fileName[256];
	float Mass,vx,vy,x,y;
	//argv[1]="20bodies.txt";
	clock_t start,finish;
	TreeNode* root;
	//strcpy_s(fileName,argv[1]);
    body** bodies_tmp=new body*[InitBody];
	//body** newlist=NULL;
    ifstream fin(FileName);

    while(fin>>Mass>>vx>>vy>>x>>y){
		bodies_tmp[count]=new body;
		bodies_tmp[count]->mass=Mass;
		bodies_tmp[count]->dx=vx;
		bodies_tmp[count]->dy=vy;
		bodies_tmp[count]->x=x;
		bodies_tmp[count]->y=y;
		count++;
	}
	if(count!=InitBody){
		bodies=new body*[count];
		for(i=0;i<count;i++){
			bodies[i]=bodies_tmp[i];
		}
		delete [] bodies_tmp;
	}
	else
		bodies=bodies_tmp;
    start=clock();
	for(i=0;i<StepCount;i++){
		root=new TreeNode(0,0,MaxX,MaxY,&bodies,count);
		ForceCalculations=0;
		count=ComputeForcesAndUpdateBodies(&bodies,root,count);
//		cout<<"Iteration "<<i<<" - Force calculations: "<<ForceCalculations<<"  count="<<count<<endl;
		delete root;
	}
    finish=clock();
	for (i=0;i<count;i++){		
		cout<<bodies[i]->mass<<" "<<bodies[i]->dx<<" "<<bodies[i]->dy<<" "<<bodies[i]->x<<" "<<bodies[i]->y<<endl;
	}
	cout<<"The duration time is: "<<(double)(finish-start)/CLOCKS_PER_SEC<<endl;
	for(i=0;i<count;i++){
		delete bodies[i];
	}
	delete [] bodies;
}