// parallel.cpp : Defines the entry point for the console application.
//
#include "stdafx.h"
#include "mpi.h"
#include "time.h"
#include < stdlib.h >
#include < iostream >
#include < math.h >
#define TIME 0.0001 //дискрет изменения времени
#define G 6.67e-11 //гравитационная постоянная
using namespace std;
struct body{
double weight;//масса тела
long double x,y,z;//координаты точки
long double vx,vy,vz;//компоненты скорости
};
void ProcessInizialization(int N,int ProcRank,int ProcNum,body* &psystem,int* &pSendNum,int *&pSendIndex);
void MoveSystem(body* system,int N,int* pSendNum,int* pSendIndex,body* psystem,int ProcNum,int ProcRank);
void PrintSystem(body* system,int N,double Duration);
long double Length(body b1,body b2);
void Motion(body &b,long double Fx,long double Fy,long double Fz);
void parallelalgorithm();
int _tmain(int argc, char** argv)
{
MPI_Init(&argc,&argv);
double time;//начальное время
int N;//число тел
double T;//длительность времени
body* system;//система тел
body* psystem;//часть тел для отдельного процессора
int ProcRank,ProcNum;//ранг процессора и число процессоров
double Start, Finish, Duration;
MPI_Barrier(MPI_COMM_WORLD);
Start=MPI_Wtime();
MPI_Comm_size(MPI_COMM_WORLD,&ProcNum);//вычисление числа процессоров
MPI_Comm_rank(MPI_COMM_WORLD,&ProcRank);//ранк процессора
if(ProcRank==0){
FILE* in;
in=fopen("inputdata.txt","r");
fscanf(in,"%d\t%lf\n",&N,&T);
system=new body[N];
int count=0;
while(!feof(in)){
system[count].vx=system[count].vy=system[count].vz=0;
system[count].weight=1e+20;
fscanf(in,"%lf\t%lf\t%lf\n",&system[count].x,&system[count].y,&system[count].z);
count++;
}
fclose(in);
}
time=0;
MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&T,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
int* pSendNum;
int* pSendIndex;
if(ProcRank!=0)system=new body[N];
ProcessInizialization(N,ProcRank,ProcNum,psystem,pSendNum,pSendIndex);//инициализация исходных данных
/////////////////////////////////////////
#define N_body 7
/* структуры данных для создания MPI-типа данных,
соответствующий С-типу body*/
MPI_Datatype body_types[N_body]={MPI_DOUBLE,
MPI_LONG_DOUBLE,MPI_LONG_DOUBLE,MPI_LONG_DOUBLE,
MPI_LONG_DOUBLE,MPI_LONG_DOUBLE,MPI_LONG_DOUBLE};
//количество элементов этого типа в каждом поле
int body_counts[N_body]={1,1,1,1,1,1,1};
//смещение каждого поля относительно начала структуры
MPI_Aint body_disp[N_body];
//переменные используемые для вычисления смещений
MPI_Aint start_address, address;
//новый MPI-тип
MPI_Datatype MPI_body;
//создаем MPI-тип для body
MPI_Address(&system[0],&start_address);
MPI_Address(&system[0].weight,&address);
body_disp[0]=address-start_address;//смещение первого поля
MPI_Address(&system[0].x,&address);
body_disp[1]=address-start_address;
MPI_Address(&system[0].y,&address);
body_disp[2]=address-start_address;
MPI_Address(&system[0].z,&address);
body_disp[3]=address-start_address;
MPI_Address(&system[0].vx,&address);
body_disp[4]=address-start_address;
MPI_Address(&system[0].vy,&address);
body_disp[5]=address-start_address;
MPI_Address(&system[0].vz,&address);
body_disp[6]=address-start_address;
MPI_Type_struct(N_body,body_counts,body_disp,body_types,&MPI_body);
//делаем созданный MPI-тип доступным при обмене сообщениями
MPI_Type_commit(&MPI_body);
////////////////////////////////////
//обобщенная рассылка данных от 0 процессора всем
MPI_Scatterv(system,pSendNum,pSendIndex,MPI_body,psystem,pSendNum[ProcRank],MPI_body,0,MPI_COMM_WORLD);
while(time < T){
//обобщенный сбор данных от всех процессоров одному
//информация об остальных телах
MPI_Allgatherv(psystem,pSendNum[ProcRank],MPI_body,system,pSendNum,pSendIndex,MPI_body,MPI_COMM_WORLD);
MoveSystem(system,N,pSendNum,pSendIndex,psystem,ProcNum,ProcRank);//движение части тел системы для данного процессора
time+=TIME;
}
MPI_Gatherv(psystem,pSendNum[ProcRank],MPI_body,system,pSendNum,pSendIndex,MPI_body,0,MPI_COMM_WORLD);
MPI_Type_free(&MPI_body);
MPI_Barrier(MPI_COMM_WORLD);
Finish=MPI_Wtime();
Duration=Finish-Start;
if(ProcRank==0)PrintSystem(system,N,Duration);
delete[] system;
delete[] psystem;
delete[] pSendIndex;
delete[] pSendNum;
MPI_Finalize();
return 0;
};
void ProcessInizialization(int N,int ProcRank,int ProcNum,body* &psystem,int* &pSendNum,int* &pSendIndex)
{
pSendNum=new int[ProcNum];
pSendIndex=new int[ProcNum];
pSendIndex[0]=0;
pSendNum[0]=N/ProcNum;
N=N-pSendNum[0];
for(int i=1;i < ProcNum;i++){
pSendIndex[i]=pSendIndex[i-1]+pSendNum[i-1];
pSendNum[i]=N/(ProcNum-i);
N=N-pSendNum[i];
}
psystem=new body[pSendNum[ProcRank]];//иницивлизация массива частиц для процессора с рангом ProcRank
};
void MoveSystem(body* system,int N,int* pSendNum,int* pSendIndex,body* psystem,int ProcNum,int ProcRank)
{
long double Fx,Fy,Fz;
for(int i=0;i < pSendNum[ProcRank];i++){
system[pSendIndex[ProcRank]+i];//текущий объект
Fx=0,Fy=0,Fz=0;//компоненты результирующей силы
for(int j=0;j < N;j++){
if(j!=pSendIndex[ProcRank]+i){
Fx+=G*system[pSendIndex[ProcRank]+i].weight*system[j].weight/pow(Length(system[pSendIndex[ProcRank]+i],system[j]),3)*(system[j].x-system[pSendIndex[ProcRank]+i].x);
Fy+=G*system[pSendIndex[ProcRank]+i].weight*system[j].weight/pow(Length(system[pSendIndex[ProcRank]+i],system[j]),3)*(system[j].y-system[pSendIndex[ProcRank]+i].y);
Fz+=G*system[pSendIndex[ProcRank]+i].weight*system[j].weight/pow(Length(system[pSendIndex[ProcRank]+i],system[j]),3)*(system[j].z-system[pSendIndex[ProcRank]+i].z);
}
}
Motion(psystem[i],Fx,Fy,Fz);
}
};
void PrintSystem(body* system,int N,double Duration)
{
FILE *out;
out=fopen("outputdataparallel.txt","w");
for(int i=0;i < N;i++)
fprintf(out,"%lf\t%lf\t%lf\n",system[i].x,system[i].y,system[i].z);
fprintf(out,"%s %lf","Duration of parallel calculating:",Duration);
fclose(out);
};
long double Length(body b1,body b2)
{
return sqrt((b1.x-b2.x)*(b1.x-b2.x)+(b1.y-b2.y)*(b1.y-b2.y)+(b1.z-b2.z)*(b1.z-b2.z));
};
void Motion(body &b,long double Fx,long double Fy,long double Fz){
long double accelx,accely,accelz;
accelx=Fx/b.weight;
accely=Fy/b.weight;
accelz=Fz/b.weight;
b.x=b.x+b.vx*TIME+accelx*TIME*TIME/2;
b.y=b.y+b.vy*TIME+accely*TIME*TIME/2;
b.z=b.z+b.vz*TIME+accelz*TIME*TIME/2;
b.vx=b.vx+accelx*TIME;
b.vy=b.vy+accely*TIME;
b.vz=b.vz+accelz*TIME;
};
