<![CDATA[#include <stdio.h>
"Gravitational N-Body Problem" task.
Parallel implementation.
---------------------------------------
---------------------------------------
---------------------------------------
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include "mpi.h"
int len,pid,pnr;
int ln,rn,un,dn;
MPI_Status status;
int AnzahlKoerper=3;
int KR=2;
const double G=2.959122083e-04;
double AlleR[100000][3]; //body's koordinats
double AlleV[100000][3]; //velocity
double AlleA[100000][3];
double Massen[100000];
double T,dt;
int driver,mode;
FILE *super_file;
void Linearebewegung(double *R, double *V, double *A, double dt, double *R1, double *V1);
double SQR(double x);
void A_Berechen();
void Integrationsschritt(double dt);
void initialisierung();
void Zeichne();
float supertime[4] = {0,0,0,0};
void Linearebewegung(double *R, double *V, double *A, double dt, double *R1, double *V1)
{
int i;
if ((pid>0)&&(pid<=pnr))
for (i=0; i<3; i++)
{
R1[i]=R[i] + V[i]*dt + 0.5*dt*dt*A[i];
V1[i]=V[i] + A[i]*dt;
}
}
double SQR(double x)
{
return x*x;
}
void A_Berechen()
{
int i,k,kk,kkk,j;
double d;
double Massen_new;
double AlleR_new[3];
double s_buffer[500000];
double r_buffer[500000];
float time;
for (i=0;i<KR;i++)
for (j=0;j<3;j++)
AlleA[i][j]=0;
i=0;
MPI_Barrier(MPI_COMM_WORLD);
if (pid==0) time = MPI_Wtime();
for (k=0;k<KR;k++)
{
s_buffer[i] = AlleR[k][0];
s_buffer[i+1] = AlleR[k][1];
s_buffer[i+2] = AlleR[k][2];
s_buffer[i+3] = Massen[k];
i+=4;
}
MPI_Allgather( s_buffer , KR*4 , MPI_DOUBLE, r_buffer, KR*4, MPI_DOUBLE, MPI_COMM_WORLD);
if (pid==0) supertime[0] = (MPI_Wtime() - time);
if (pid==0) time = MPI_Wtime();
for (k=1 ; k < pnr ; k++)
{
if ( k != pid )
{
for (i=0;i<KR;i++)
{
AlleR_new[0]=r_buffer[k*KR*4 + i*4];
AlleR_new[1]=r_buffer[k*KR*4 + i*4 + 1];
AlleR_new[2]=r_buffer[k*KR*4 + i*4 + 2];
Massen_new =r_buffer[k*KR*4 + i*4 + 3];
for (kk=0;kk<KR;kk++)
{
d = sqrt( SQR( AlleR[kk][0]-AlleR_new[0] ) + SQR( AlleR[kk][1]-AlleR_new[1] ) + SQR( AlleR[kk][2]-AlleR_new[2] ) );
for (kkk=0; kkk < 3; kkk++)
AlleA[kk][kkk]+= G*Massen_new * (AlleR_new[kkk] - AlleR[kk][kkk] ) / (d*d*d);
}
}
}
}
for (i=0; i<KR; i++)
for (j=0; j<KR; j++)
{
if (i!=j)
{
d = sqrt( SQR( AlleR[i][0]-AlleR[j][0] ) + SQR( AlleR[i][1]-AlleR[j][1] ) + SQR( AlleR[i][2]-AlleR[j][2] ) );
for (k=0; k < 3; k++)
AlleA[i][k]+= G*Massen[j] * (AlleR[j][k] - AlleR[i][k] ) / (d*d*d);
}
}
if (pid==0) supertime[1] = (MPI_Wtime() - time);
}
void Integrationsschritt(double dt)
{
int i,j;
float time;
A_Berechen();
if (pid==0) time = MPI_Wtime();
for (i=0; i<KR; i++)
Linearebewegung(&AlleR[i][0],&AlleV[i][0],&AlleA[i][0], dt,&AlleR[i][0],&AlleV[i][0]);
if (pid==0) supertime[2] = (MPI_Wtime() - time);
}
void initialisierung()
{
int i,k,j;
float f;
for (k=0; k<100000;k++)
for (i=0; i<3;i++)
{
AlleR[k][i]=0;
AlleV[k][i]=0;
AlleA[k][i]=0;
}
super_file = fopen("NBODY.data","r");
fscanf(super_file,"%d",&AnzahlKoerper);
KR=AnzahlKoerper/(pnr-1);
j=1;
for (i=0; i<AnzahlKoerper; i+=KR)
{
for (k=0;k<KR;k++)
{
fscanf(super_file,"%f",&f);
if (j==pid)
{
Massen[k]=f;
}
}
j++;
}
j=1;
for (i=0; i<AnzahlKoerper ;i+=KR)
{
for (k=0;k<KR;k++)
{
fscanf(super_file,"%f",&f);
if (j==pid)
AlleR[k][0]=f;
}
j++;
}
j=1;
for (i=0; i<AnzahlKoerper; i+=KR)
{
for (k=0;k<KR;k++)
{
fscanf(super_file,"%f",&f);
if (j==pid)
AlleV[k][1]=f;
}
j++;
}
dt=0.01;
T=0;
if (pid==0)
{
openwindow();
line(0,500,1400,500);
}
}
void Zeichne()
{
int i,x,y,k;
int buffer_xy[100000];
int send_buffer[1000000]; //x,y
i=0;
if (pid != 0)
for (k=0;k<KR;k++)
{
x = (int)((200*AlleR[k][0]) + 700);
y = (int)((200*AlleR[k][1]) + 500);
if (x >= 1390) AlleV[k][0] = -AlleV[k][0];
if (y >= 990) AlleV[k][1] = -AlleV[k][1];
if (x <= 10) AlleV[k][0] = -AlleV[k][0];
if (y <= 10) AlleV[k][1] = -AlleV[k][1];
send_buffer[i]=x;
i++;
send_buffer[i]=y;
i++;
}
MPI_Gather(send_buffer, KR*2, MPI_INT, buffer_xy, KR*2, MPI_INT, 0, MPI_COMM_WORLD);
if (pid==0)
for (i=KR*2 ; i <= (AnzahlKoerper+1)*2 ; i+=2)
putpixel(buffer_xy[i],buffer_xy[i+1]);
}
int main(int argc, char **argv)
{
int i=0,j=1;
float time;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&pnr);
MPI_Comm_rank(MPI_COMM_WORLD,&pid);
initialisierung();
// MPI_Barrier(MPI_COMM_WORLD);
// if (pid == 0) time = MPI_Wtime();
for ( i=0;i<j;i++ )
{
Integrationsschritt(dt);
// Zeichne();
}
// if (pid==0)
// printf("\n Time of computing = %f sec", MPI_Wtime() - time);
if (pid==0)
{
printf("\n Time of transmission = %f sec", supertime[0]/j);
printf("\n Time of computing = %f sec", supertime[1]+supertime[2]);
}
MPI_Finalize();
}
]]>