RK4 w/ MPI

rk4_mpi.c

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <mpi.h>

#define formaster if(id==0)
#define forothers if(id!=0)
#define dim 2
#define USE_DEFAULT     1		// skip user input
#define DO_HARLEM_SHAKE 1		// show the result

double f(double t, int index)
{
	return 1.;							// here you should return whatever f(t) should be equal to
}

int main(int argc, char* argv[])
{
	int id;
	int p;
	double k1, k2, k3, k4;
	int n, i, j, k;
	double t, h = 0.0001, h_2, ti = 0., tf = 1.;
	double *y, *sub_y;
	double **A;
	int sub_dim;
	int *sendcounts;	// array describing how many elements to send to each process
	int *displs;		// array describing the displacements where each segment begins
	int rem;			// elements remaining after division among processes
	time_t tic, toc;						// measure time
	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &id);		// get the rank of this processor
	MPI_Comm_size(MPI_COMM_WORLD, &p);		// get the number of processors
	sub_dim = dim/p + (dim%p!=0);
	
	// allocate required memory
	A = calloc(dim, sizeof(*A));
	for(i=0; i<dim; ++i)
		A[i] = calloc(dim, sizeof(**A));
	y = calloc(dim, sizeof(*y));
	sub_y = calloc(sub_dim, sizeof(*sub_y));
	
	formaster
	{
		printf("%s", "Enter the initial value of t: ");
		#if USE_DEFAULT
			printf("%lf\n", ti);
		#else
			scanf("%lf", &ti);
			MPI_Bcast(&ti, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		#endif // USE_DEFAULT
		printf("%s", "Enter the final value of t: ");
		#if USE_DEFAULT
			printf("%lf\n", tf);
		#else
			scanf("%lf", &tf);
			MPI_Bcast(&tf, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		#endif // USE_DEFAULT
		printf("%s", "Enter the value of h (step): ");
		#if USE_DEFAULT
			printf("%lf\n", h);
		#else
			scanf("%lf", &h);
			MPI_Bcast(&h, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		#endif // USE_DEFAULT
		for(i=1; i<=dim; ++i)
		{
			printf("Enter the initial value of the %i%s component of y: ", i, ((((i%10)==1)&&((i%100)!=11))?"st":((((i%10)==2)&&((i%100)!=12))?"nd":"th")));
			#if USE_DEFAULT
				y[i-1] = (double)i;
				printf("%lf\n", y[i-1]);
			#else
				scanf("%lf", y+i-1);
			#endif // USE_DEFAULT
		}
		printf("\n");
		
		if(p==1)
			printf("single thread execution\n");
		else
			printf("multiple thread execution (%i procs.)\n", p);
	}

	n = (int)((tf-ti)/h);
	h_2 = h/2.;
	
	sendcounts = calloc(p, sizeof(*sendcounts));
	displs = calloc(p, sizeof(*displs));
	rem = dim%p;
	// calculate send counts and displacements
	for(i=0; i<p; ++i)
	{
		sendcounts[i] = dim/p;
		if(rem --> 0)
			++sendcounts[i];
		displs[i] = i?(displs[i-1]+sendcounts[i-1]):0;
	}
	
	// initialize A
	for(i=0; i<dim; ++i)
		for(j=0; j<dim; ++j)
			A[i][j]=(i==j)*(1-2*(i&1));

	#if DO_HARLEM_SHAKE
		formaster
		{
			printf("t");
			for(j=0; j<dim; ++j)
				printf("\ty(%i)", j+1);
			printf("\n");
		}
	#endif // DO_HARLEM_SHAKE
	i = 0;
	MPI_Barrier(MPI_COMM_WORLD);
	formaster
	{
		// tic
		tic = time(NULL);
	}
	while(i<=n)
	{
		t = ti + h*i;
		#if DO_HARLEM_SHAKE
			formaster
			{
				printf("%lf", t);
				for(j=0; j<dim; ++j)
					printf("\t%lf", y[j]);
				printf("\n");
			}
		#endif // DO_HARLEM_SHAKE
		// paralleled calculation
		MPI_Bcast(y, dim, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		MPI_Scatterv(y, sendcounts, displs, MPI_DOUBLE, sub_y, sub_dim, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		for(j=displs[id]; j<displs[id]+sendcounts[id]; ++j)
		{
			k1 = f(t    , j);
			for(k=0; k<dim; ++k)
				k1+=A[j][k]* y[k]        ;
			k1 *= h;
			k2 = f(t+h_2, j);
			for(k=0; k<dim; ++k)
				k2+=A[j][k]*(y[k]+0.5*k1);
			k2 *= h;
			k3 = f(t+h_2, j);
			for(k=0; k<dim; ++k)
				k3+=A[j][k]*(y[k]+0.5*k2);
			k3 *= h;
			k4 = f(t+h  , j);
			for(k=0; k<dim; ++k)
				k4+=A[j][k]*(y[k]+    k3);
			k4 *= h;
			sub_y[j-displs[id]] += (k1+k4)/6.+(k2+k3)/3.;
		}
		MPI_Gatherv(sub_y, sendcounts[id], MPI_DOUBLE, y, sendcounts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD);
		// wait till all threads are completed
		MPI_Barrier(MPI_COMM_WORLD);
		// and then go to a next step
		++i;
	}
	MPI_Barrier(MPI_COMM_WORLD);
	formaster
	{
		// toc
		toc = time(NULL);
		printf("the calculation took %.0lf seconds\n", difftime(toc, tic));
	}
	free(displs);
	free(sendcounts);
	free(sub_y);
	free(y);
	for(i=0; i<dim; ++i)
		free(A[i]);
	free(A);
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Finalize();
	formaster
		printf("done\n");
	return 0;
}