This is the MPI solution of gravitation problem

Gravit.cpp

#include "pt4.h"
#include "mpi.h"
#include <ctime>
#include <cmath>

#define gravity 10 // гравитационная постоянная
#define dt 0.1 // шаг по времени
#define N 800 // количество частиц
#define fmax 1 // максимальное значение силы
#define Niter 100 // число итераций


struct Particle
{
	double x, y, vx, vy;
	void free() {
		x = y = vx = vy = 0;
	}

};
struct Force
{
	double x, y;
	void free() {
		x = y = 0;
	}
};

MPI_Datatype MPI_PARTICLE;
MPI_Datatype MPI_FORCE;
int maxBlockSize;
int k, r;
double * m;
int * block_size;
int * block_start;
Particle * p, *tp;
Force * f, *tf;

template<class T>
void free_(T * ar, int sz) {
	for (int i = 0; i < sz; ++i) {
		ar[i].free();
	}
}

void Init()
{
	for (int i = 0, j = block_start[r]; i < block_size[r]; ++i, ++j) {
		p[i].x = 20 * (j / 20 - 20) + 10;
		p[i].y = 20 * (j % 20 - 10) + 10;
		//Show(p[i].x); ShowLine(p[i].y);
		p[i].vx = p[i].y / 15;
		p[i].vy = -p[i].x / 50;
		f[i].x = 0;
		f[i].y = 0;
	}
	for (int i = 0; i < N; i++)
	{
		m[i] = 100 + i % 100;
	}
}

void check_output() {
	if (r == 0) {
		ShowLine();
		Show(p[0].x); Show(p[0].y);
	}
	else if (r == k - 1) {
		ShowLine();
		Show(p[block_size[r] - 1].x); Show(p[block_size[r] - 1].y);
	}
}

void CalcForces_inner()
{
	for (int i = 0; i < block_size[r] - 1; i++)
		for (int j = i + 1; j < block_size[r]; j++)
		{
			double dx = p[j].x - p[i].x, dy = p[j].y - p[i].y;
			double r_2 = 1 / (dx * dx + dy * dy);
			double r_1 = sqrt(r_2);
			double fabs = gravity * m[i + block_start[r]] * m[j + block_start[r]] * r_2;
			

			if (fabs > fmax) 
				fabs = fmax;
			f[i].x += dx = fabs * dx * r_1;
			f[i].y += dy = fabs * dy * r_1;
			f[j].x -= dx;
			f[j].y -= dy;
		}
}

void CalcForces_outer(int outer_p)
{
	for (int i = 0; i < block_size[r]; i++)
		for (int j = 0; j < block_size[outer_p]; j++)
		{
			double dx = tp[j].x - p[i].x, dy = tp[j].y - p[i].y;
			double r_2 = 1 / (dx * dx + dy * dy);
			double r_1 = sqrt(r_2);
			double fabs = gravity * m[i + block_start[r]] * m[j + block_start[outer_p]] * r_2;

			if (fabs > fmax)
				fabs = fmax;
			f[i].x += dx = fabs * dx * r_1;
			f[i].y += dy = fabs * dy * r_1;
			tf[j].x -= dx;
			tf[j].y -= dy;
		}
}

void MoveParticlesAndFreeForces()
{
	for (int i = 0; i < block_size[r]; i++)
	{
		double dvx = f[i].x * dt / m[i + block_start[r]],
			dvy = f[i].y * dt / m[i + block_start[r]];
		p[i].x += (p[i].vx + dvx / 2) * dt;
		p[i].y += (p[i].vy + dvy / 2) * dt;
		p[i].vx += dvx;
		p[i].vy += dvy;
		f[i].x = 0;
		f[i].y = 0;
	}
}



void DoIter() {
	free_(f, block_size[r]);
	free_(tf, maxBlockSize);
	free_(tp, maxBlockSize);
	
	MPI_Request req;
	MPI_Status state;
	for (int i = 0; i < r; ++i) {
		MPI_Isend(p, block_size[r], MPI_PARTICLE, i, 1, MPI_COMM_WORLD, &req);
	}
	
	CalcForces_inner();

	for (int i = r + 1; i < k; ++i) {
		
		MPI_Recv(tp, block_size[i], MPI_PARTICLE, i, 1, MPI_COMM_WORLD, &state);
		
		CalcForces_outer(i);
		MPI_Send(tf, block_size[i], MPI_FORCE, i, 1, MPI_COMM_WORLD);
		free_(tf, maxBlockSize);
	}
	
	for (int i = 0; i < r; ++i) {
		MPI_Recv(tf, block_size[r], MPI_FORCE, i, 1, MPI_COMM_WORLD, &state);
		for (int j = 0; j < block_size[r]; ++j) {
			f[j].x += tf[j].x;
			f[j].y += tf[j].y;
		}
	}

	MoveParticlesAndFreeForces();
}

void Solve()
{
    Task("MPIDebug8");
    int flag;
    MPI_Initialized(&flag);
    if (flag == 0)
        return;
    int rank, size;
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	HideTask();
	SetPrecision(12);

	k = size;
	r = rank;

	clock_t timer = clock();

	block_size = new int[k];
	block_start = new int[k];

	int piece_size = N*(N - 1) / 2 / k;
	int ctr = 0, j = 0;
	block_start[0] = 0;
	for (int i = 0; i < N; ++i) {
		ctr += N - i - 1;
		if(ctr > piece_size){
			ctr = N - i - 1;
			block_size[j] = i - block_start[j];
			if(j + 1 < k)
				block_start[j + 1] = i;
			++j;
		}
	}
	
	maxBlockSize = block_size[k - 1] = N - block_start[k - 1];

	m = new double[N];
	p = new Particle[block_size[r]];
	f = new Force[block_size[r]];
	tp = new Particle[maxBlockSize];
	tf = new Force[maxBlockSize];

	MPI_Type_vector(1, 4, 4, MPI_DOUBLE, &MPI_PARTICLE);
	Particle tprt;
	MPI_Type_struct(1, new int[1]{ 4 }, new int[1]{ 0 }, new MPI_Datatype[1]{ MPI_DOUBLE }, &MPI_PARTICLE);
	MPI_Type_commit(&MPI_PARTICLE);
	MPI_Type_struct(1, new int[1]{ 2 }, new int[1]{ 0 }, new MPI_Datatype[1]{ MPI_DOUBLE }, &MPI_FORCE);
	MPI_Type_commit(&MPI_FORCE);
	Init();

	for (int i = 0; i < Niter; ++i) {
		DoIter();
	}

	timer = clock() - timer;
	Show("Time: ", ((double)timer / CLOCKS_PER_SEC) * 1000);
	int pc = (N - block_start[r] + N - block_start[r] - block_size[r]) * block_size[r]/2;
	Show("PairsCount: ", pc);
	check_output();
};