nbody.c

一个在linux环境下的基于mpi＋openmp的混合编程

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <mpi.h>
#include<omp.h>
#include <sys/types.h>
#include <stdarg.h>
#include <string.h>
#include <errno.h>
#include <fcntl.h>

#define NUM_THREADS 4
typedef struct {
    double x;
    double y;
} coords;

struct body_s {
    coords p;
    coords v;
    coords f;
    double m;
    struct body_s *next;
    struct body_s *prev;
};
typedef struct body_s body;

typedef struct {
  double px, py;
  double vx, vy;
  double fx, fy;
  double m;
} Body;
MPI_Datatype MPI_BODY;
MPI_Status status;

int my_rank;
int size;
Body* bodies;
int* myStart;
int* myEnd;
  
static int N;
static double width, height,G,DT,C;

static body *blist_head = NULL;
static body *blist_tail = NULL;

static body **barray;

/* 插入body.
 */
int blist_insert(body *new)
{
    /*
     * Insert node at the tail of the list.
     */
    new->next = NULL;
    new->prev = blist_tail;
    if (blist_tail != NULL)
	blist_tail->next = new;
    if (blist_head == NULL) 
	blist_head = new;
    blist_tail = new;
	
    return 0;
}

/*读入每一行*/
int readline(FILE *fp, void *vptr, size_t maxlen) {
    int n, rc;
    char c, *ptr;

    ptr = vptr;
    /* n = 1, because we need space to store '\0' */
    for (n = 1; n < maxlen; n++) {
	if ((rc = fread(&c, 1, 1, fp)) == 1) { /* one byte read */
	    *ptr++ = c;

	    if (c == '\n') /* end of line */
		break;
	}
	else if (rc == 0) {
	    if (n == 1)
		/* EOF, no data read */
		return 0;
	    else
		/* EOF, data read */
		break;
	}
    }
    
    *ptr = '\0';
    return n;
}

#define LINE_SIZE 250   /* Maximum size of a line in the data file */
#define BODY_ENTRIES 7  /* Number of entries per body (as described belov) */

void read_input(const char *filename) {
    FILE *fp = fopen(filename, "r");
    char line[LINE_SIZE], *p, *p_new;
    int current_line = 0;
    body *new;
    int i;
    body *bp;
    
    /******************************************/
    /*     Read base information              */
    /******************************************/
    
    /*read N*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%d", &N);
	printf("N:%d\n",N);
    /*read G*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf", &G);
	printf("G:%lf\n",G);
    /*read Size*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf %lf", &width, &height);
	printf("Size:%lf %lf\n",width,height);
    /*read DT*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf", &DT);
	printf("DT:%lf\n",DT);
    /*read c*/
    if (readline(fp, (void *)line, LINE_SIZE) == 0)
	      printf("Invalid input file (first line should contain universe size)");
    sscanf(line, "%lf", &C);
	printf("C:%lf\n",C);
    /* Number of bodies */

    
    while (readline(fp, (void *)line, LINE_SIZE) != 0) {

	/* Check line for correctnes (number of entries) */
	current_line++;
//	if (check_line(line) == -1)
//	    printf("Error at line %d: %s\n", current_line, line);

	new = malloc(sizeof(body));

	/* Parse line */
	p = line;
	/* Goto next non-blank character */
	while (isspace(*p) && *p != '\0') {
	    *p = '\0';
	    p++;
	    printf(".");
	}

	/* Set entries */
	new->m = strtod(p, &p_new);
	p = p_new;
	new->p.x = strtod(p, &p_new);
	p = p_new;
	new->p.y = strtod(p, &p_new);
	p = p_new;
	new->v.x = strtod(p, &p_new);
	p = p_new;
	new->v.y = strtod(p, &p_new);
		
    
	/* Insert in list */
	blist_insert(new);
    }	
    
    //N = current_line; /* Number of bodies */
    bodies = (Body *) malloc(sizeof(Body) * N);
    for (i = 0, bp = blist_head; i < N; i++, bp = bp->next) {
      bodies[i].px = bp->p.x;
      bodies[i].py = bp->p.y;
      bodies[i].vx = bp->v.x;
      bodies[i].vy = bp->v.y;
      bodies[i].fx = 0.0;
      bodies[i].fy = 0.0;
      bodies[i].m = bp->m;
    }

    /* Initialize barray */
    barray = malloc(sizeof(body *) * N);
    for (i = 0, bp = blist_head; i < N; i++, bp = bp->next)
	barray[i] = bp;
}

/* 计算力 */
void calculate_forces(void) {
  int i, j;
  coords direction;
  double distance,distance2, magnitude;
    omp_set_num_threads(NUM_THREADS);
#pragma omp parallel for default(shared) private(i,j,direction,distance,distance2,magnitude)
  for(i = myStart[my_rank]; i < myEnd[my_rank]; i++) {
    for (j = 0; j < N; j++) {
      if (i == j) continue;
	  distance2 = pow((bodies[i].px - bodies[j].px), 2.0) + 
		       pow((bodies[i].py - bodies[j].py), 2.0);
      distance = sqrt( distance2 );
      magnitude = (G * bodies[i].m * bodies[j].m) / (distance2+pow(C, 2.0));
      direction.x = bodies[j].px - bodies[i].px;
      direction.y = bodies[j].py - bodies[i].py;

      bodies[i].fx = bodies[i].fx + magnitude * direction.x / distance;
      bodies[i].fy = bodies[i].fy + magnitude * direction.y / distance;
    }
  }
}
/*计算物体的位置*/
void move_bodies(void) {
  int i;
  coords deltav;
  coords deltap;
  double tempPointx,tempPointy;
    omp_set_num_threads(NUM_THREADS);
  #pragma omp parallel for default(shared) private(i,deltav,deltap,tempPointx,tempPointy)
  for (i = myStart[my_rank]; i < myEnd[my_rank]; i++) {
    deltav.x = bodies[i].fx / bodies[i].m * DT;
    deltav.y = bodies[i].fy / bodies[i].m * DT;
    deltap.x = (bodies[i].vx + deltav.x / 2) * DT;
    deltap.y = (bodies[i].vy + deltav.y / 2) * DT;

    tempPointx = bodies[i].px + deltap.x;
    tempPointy = bodies[i].py + deltap.y;
	//运行轨迹超出了空间范围
	if(tempPointx < 0 || tempPointx > width ||
	tempPointy < 0 || tempPointy > height){
	bodies[i].vx = -bodies[i].vx;
    bodies[i].vy = -bodies[i].vy;
    bodies[i].px = bodies[i].px;
    bodies[i].py = bodies[i].py;
	}
	else
	{
    bodies[i].vx = bodies[i].vx + deltav.x;
    bodies[i].vy = bodies[i].vy + deltav.y;
    bodies[i].px = tempPointx;
    bodies[i].py = tempPointy;
	}




    /* Reset force vectors */
    bodies[i].fx = 0.0;
    bodies[i].fy = 0.0;

  }
}

/* Run simulation for 'steps' steps. Results are writen to res_filename. */
void run(int steps, char *res_filename) {
  int step, j;
  FILE *fp;

  if (my_rank == 0) {
    printf("Running %i bodies, %i steps\n", N, steps);
    fp = fopen(res_filename, "w");
  }

  for (step = 0; step < steps; step++) {
    MPI_Bcast(bodies, N, MPI_BODY, 0, MPI_COMM_WORLD);


    /* Calculate */
    calculate_forces();
    move_bodies();

    if (my_rank == 0) {
      for (j = 1; j < size; j++)
		  MPI_Recv(&bodies[myStart[j]],
					  myEnd[j] - myStart[j],
					  MPI_BODY, j, 66, MPI_COMM_WORLD, &status);
      
    } else {
      MPI_Send(&bodies[myStart[my_rank]],
	       myEnd[my_rank] - myStart[my_rank],
	       MPI_BODY, 0, 66, MPI_COMM_WORLD);
    }
  }
  if (my_rank == 0) {
	  //输出结果;
      fprintf(fp, "%d\n",N);
      fprintf(fp, "%lf\n",G);
      fprintf(fp, "%lf %lf\n",width,height);
      fprintf(fp, "%lf\n",DT);
      fprintf(fp, "%lf\n",C);
      for (j = 0; j < N; j++)
				fprintf(fp, "%lf %lf %lf %lf %lf\n",bodies[j].m, bodies[j].px, bodies[j].py, bodies[j].vx, bodies[j].vy);
    fclose(fp);
  }
}

int main(int argc, char *argv[]) {
  int i, z;
  double s_time, e_time;


  if (argc < 4) {
    printf("Usage: nbody INPUT_FILE STEPS RESULT_FILE");
    return -1;
  }
  
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  MPI_Comm_size(MPI_COMM_WORLD, &size);
  MPI_Type_contiguous(7, MPI_DOUBLE, &MPI_BODY);
  MPI_Type_commit(&MPI_BODY);

  if (my_rank == 0) {
    read_input(argv[1]);

    for (z = 1; z < size; z++) 
	{
		MPI_Send(&N, 1, MPI_INT, z, 66, MPI_COMM_WORLD);
		MPI_Send(&G, 1, MPI_DOUBLE, z, 67, MPI_COMM_WORLD);
		MPI_Send(&width, 1, MPI_DOUBLE, z, 68, MPI_COMM_WORLD);
		MPI_Send(&height, 1, MPI_DOUBLE, z, 69, MPI_COMM_WORLD);
		MPI_Send(&DT, 1, MPI_DOUBLE, z, 70, MPI_COMM_WORLD);
		MPI_Send(&C, 1, MPI_DOUBLE, z, 71, MPI_COMM_WORLD);
	}
  } else {
    MPI_Recv(&N, 1, MPI_INT, 0, 66, MPI_COMM_WORLD, &status);
	MPI_Recv(&G, 1, MPI_DOUBLE, 0, 67, MPI_COMM_WORLD, &status);
	MPI_Recv(&width, 1, MPI_DOUBLE, 0, 68, MPI_COMM_WORLD, &status);
	MPI_Recv(&height, 1, MPI_DOUBLE, 0, 69, MPI_COMM_WORLD, &status);
	MPI_Recv(&DT, 1, MPI_DOUBLE, 0, 70, MPI_COMM_WORLD, &status);
	MPI_Recv(&C, 1, MPI_DOUBLE, 0, 71, MPI_COMM_WORLD, &status);

    bodies = (Body *) malloc(sizeof(Body) * N);
    
    memset(bodies, 0, sizeof(Body) * N);
  }

  myStart = (int *) malloc(size * sizeof(int));
  myEnd = (int *) malloc(size * sizeof(int));

  myStart[0] = 0;
  myEnd[0] = N / size;
  for (i = 1; i < size; i++) {
    myStart[i] = myEnd[i - 1];
    myEnd[i] = myStart[i] + N / size;
  }
  myEnd[size - 1] = N;

  if (my_rank == 0) {
    for (i = 0; i < size; i++) {
      fprintf(stderr, "Process %i owns particle %i through %i\n", i, myStart[i], myEnd[i]);
    }
  }


  s_time = MPI_Wtime();
  run(atoi(argv[2]), argv[3]);
  e_time = MPI_Wtime();

  if (my_rank == 0) printf("\nProcessing took %f seconds\n", e_time - s_time);
  MPI_Finalize();
  return 0;
}