JetCracker

Life-time learner's blog

[MPI] Calculating integral in parallel

Problem: write a programme to calculate the integral with MPI:

Actually, the programme must be able to calculate any integral. But before writing a parallel programme, we should write a serial one to make sure it works and the algorithm is correct.

There are plenty of algorithms for calculating integrals. One of the most advanced (efficient), among non-parallel, is the method of local stack.

Serial programme: Local stack

/*
 * MIPT. Parallel programming.
 * Problem #3. Computing the integral of (sin(x)) / x
 * Serial algorithm: Local stack
 * Author: Anton Danshin <anton.danshin@frtk.ru>
 * Date: 15 May, 2012
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define eps 1e-10

// user function to integrate
double f(double x) {
	return sin(x)/x;
}

typedef struct {
	double a, b, fa, fb, s;
} Entry;

Entry stack[100000];
int sp = 0;

double abs (double x) {
    if(x>0) return x;
    else return -x;
}

inline int empty() {
	//puts("empty?");
	return sp <= 0;
}

inline void put(double a, double b, double fa, double fb, double s) {
	//puts("put");
	stack[sp].a = a;
	stack[sp].b = b;
	stack[sp].fa = fa;
	stack[sp].fb = fb;
	stack[sp].s = s;
	sp++;
}

inline void get(double &a, double &b, double &fa, double &fb, double &s) {
	//puts("get");
	if(!empty()) {
		sp--;
		a = stack[sp].a;
		b = stack[sp].b;
		fa = stack[sp].fa;
		fb = stack[sp].fb;
		s = stack[sp].s;
	}
}

double integrate(double a, double b) {
	double j = 0.0;
	double fa = f(a);
	double fb = f(b);
	double Sab = (fa+fb)*(b-a)*0.5;
	while(1) {
		//printf("%f %f\n", a, b);
		double c = (a+b)*0.5;
		double fc = f(c);
		double Sac = (fa+fc)*(c-a)*0.5;
		double Scb = (fb+fc)*(b-c)*0.5;
		double Sacb = Sac + Scb;
		//printf("%f %f\n", Sab, Sacb);
		if(abs(Sab-Sacb) >= eps*abs(Sacb)) {
			put(a, c, fa, fc, Sac);
			a = c;
			fa = fc;
			Sab = Scb;
		} else {
			j += Sacb;
			if(empty())
				break;
			get(a, b, fa, fb, Sab);
		}
	}
	return j;
}

int main() {
		printf("%.10f\n",integrate(0.01, 0.5));
}

Parallel programme: Collective solution method + Local stack

UPDATE: I’ve fixed some bugs!

/*
 * MIPT. Parallel programming.
 * Problem #3. Computing the integral of (sin(x)) / x
 * Parallel algorithm: Collective solution with local stack
 * Author: Anton Danshin <anton.danshin@frtk.ru>
 * Date: 16 May, 2012
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <mpi/mpi.h>

#define eps 1e-10
#define N 1000

// user function to integrate
double f(double x) {
	return sin(x)/x;
}

typedef struct {
	double a, b, fa, fb, s;
} Entry;

Entry stack[100000];
int sp = 0;

double abs (double x) {
    if(x>0) return x;
    else return -x;
}

inline bool empty() {
	//puts("empty?");
	return sp <= 0;
}

inline void put(double a, double b, double fa, double fb, double s) {
	//puts("put");
	stack[sp].a = a;
	stack[sp].b = b;
	stack[sp].fa = fa;
	stack[sp].fb = fb;
	stack[sp].s = s;
	sp++;
}

inline void get(double &a, double &b, double &fa, double &fb, double &s) {
	//puts("get");
	if(!empty()) {
		sp--;
		a = stack[sp].a;
		b = stack[sp].b;
		fa = stack[sp].fa;
		fb = stack[sp].fb;
		s = stack[sp].s;
	}
}

double integrate(double a, double b) {
	double j = 0.0;
	double fa = f(a);
	double fb = f(b);
	double Sab = (fa+fb)*(b-a)*0.5;
	while(1) {
		//printf("%f %f\n", a, b);
		double c = (a+b)*0.5;
		double fc = f(c);
		double Sac = (fa+fc)*(c-a)*0.5;
		double Scb = (fb+fc)*(b-c)*0.5;
		double Sacb = Sac + Scb;
		//printf("%f %f\n", Sab, Sacb);
		if(abs(Sab-Sacb) >= eps*abs(Sacb)) {
			put(a, c, fa, fc, Sac);
			a = c;
			fa = fc;
			Sab = Scb;
		} else {
			j += Sacb;
			if(empty())
				break;
			get(a, b, fa, fb, Sab);
		}
	}
	return j;
}

double * wtime = NULL;
int WORLD_SIZE = 1;

void send_task(int p, double a, double b) {
	MPI_Send(&a, 1, MPI_DOUBLE, p, 0, MPI_COMM_WORLD);
	MPI_Send(&b, 1, MPI_DOUBLE, p, 0, MPI_COMM_WORLD);
}

void print_stats() {
	int i;
	wtime[0] = wtime[1];
	puts("P\tTime\tAcc.\tEff.");
	for(i=1;i<=WORLD_SIZE;i++)
		printf("%i\t%.3f\t%.3f\t%.3f\n", i, wtime[i], wtime[0]/wtime[i], wtime[0]/wtime[i]/(double)i);
}

int main(int argc, char ** argv) {
	int world_size, world_rank;
	double A = 1e-5, B = 6, J = 0.0, s = 0, dx = (B-A)/N;
	double startwtime, endwtime;

	int i, k;

	// Initialize the MPI environment
	MPI_Init(&argc, &argv);
	// Find out rank, size
	MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &WORLD_SIZE);

	if(world_rank==0) {
		wtime = (double *)malloc(sizeof(double)*(WORLD_SIZE+1));
	}

	for(world_size=1;world_size<=WORLD_SIZE;world_size++) {
		// MAIN LOGIC BEGINS
		if(world_rank==0) { // master
			double a = A, b = B;
			J = 0; s = 0; k = 0; i=1;
			startwtime = MPI_Wtime();
			if(world_size==1) // we are the only machine
				J = integrate(a, b);
			else { // we have some slaves to do our job
				MPI_Status status;
				for(i=1;i<world_size;i++) {
					// send task
					send_task(i, a+k*dx, a+(k+1)*dx);
					k++; // the number of processed intervals
				}
				while(k<N) {
					MPI_Recv(&s, 1, MPI_DOUBLE, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);
					//printf("Received %f from %i; J=%f\n", s, status.MPI_SOURCE, J+s);
					J += s;
					send_task(status.MPI_SOURCE, a+k*dx, a+(k+1)*dx);
					k++;
				}
				for(i=1;i<world_size;i++) {
					MPI_Recv(&s, 1, MPI_DOUBLE, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
					J += s;
					send_task(i, 1, 0); // terminate slave #i
				}
			}
		} else { // slave
			while(1) {
				double a, b;
				MPI_Recv(&a, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
				MPI_Recv(&b, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
				if(b<a) break;
				s = integrate(a, b);
				MPI_Send(&s, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
			}
		}
		// MAIN LOGIC END

		if(world_rank==0) {
			endwtime = MPI_Wtime();
			wtime[world_size] = (endwtime-startwtime)*1000;
			printf("%i nodes. Result: %.10f\n",world_size, J);
		}

	}

	// Finalization
	if(world_rank == 0) {
		for(i=1;i<WORLD_SIZE;i++) {
			send_task(i, 1, 0); // terminate slave #i
		}
		print_stats();
		free(wtime);
	}

	MPI_Finalize();

	return 0;
}

Constants N and eps are predefined but you can modify the code to enter the values from console or file.

I’ve tested it on our training 16-core cluster.

1 nodes. Result: 1.4246775513
2 nodes. Result: 1.4246775513
3 nodes. Result: 1.4246775513
4 nodes. Result: 1.4246775513
5 nodes. Result: 1.4246775513
6 nodes. Result: 1.4246775513
7 nodes. Result: 1.4246775513
8 nodes. Result: 1.4246775513
9 nodes. Result: 1.4246775513
10 nodes. Result: 1.4246775513
11 nodes. Result: 1.4246775513
12 nodes. Result: 1.4246775513
13 nodes. Result: 1.4246775513
14 nodes. Result: 1.4246775513
15 nodes. Result: 1.4246775513
16 nodes. Result: 1.4246775513
P	Time	Acc.	Eff.
1	97.573	1.000	1.000
2	79.513	1.227	0.614
3	40.011	2.439	0.813
4	26.728	3.651	0.913
5	20.020	4.874	0.975
6	16.121	6.053	1.009
7	13.491	7.232	1.033
8	11.684	8.351	1.044
9	10.296	9.477	1.053
10	9.274	10.521	1.052
11	8.531	11.438	1.040
12	7.929	12.306	1.025
13	7.450	13.097	1.007
14	7.249	13.460	0.961
15	6.966	14.007	0.934
16	6.939	14.062	0.879

Good luck! :)

About these ads

2 responses to “[MPI] Calculating integral in parallel

  1. Pingback: Spring 2012: Overall stats « JetCracker

  2. USER_MPI January 17, 2013 at 18:45

    What command did you use to compile it??

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 245 other followers

%d bloggers like this: