Skip to content

Instantly share code, notes, and snippets.

@StSav012
Created November 29, 2016 12:58
Show Gist options
  • Save StSav012/10ebc88f864b4ce1b562f16aa25dda9b to your computer and use it in GitHub Desktop.
Save StSav012/10ebc88f864b4ce1b562f16aa25dda9b to your computer and use it in GitHub Desktop.
RK4 w/ MPI
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment