I have written a simple code for Runge-Kutta fourth order integration to solve a system of ordinary differential equations and parallelized it using OpenMP. I don't know if it is the best we can do for maximum performance of the code with little effort.
I need all values of to be returned, so I kept values in all steps.
I also create threads in each time step and paralleled in position, i.e.
pragma opm parallel
is inside the loop over time.
Here is my try: (link to gitlab repository)
//RHS for a system of equations
void xprsys(const int n,const vector<double>& x, vector<double>& f)
{
/*
* n : number of equations
* x : value in each time step
* f : RHS of equtions dx/dt = f
*/
double sum1=0;
#pragma omp parallel for reduction(+:sum1)
for (int i=0; i<n; i++){
sum1 = 0;
for(int j=0; j<n; j++)
sum1 += sin(x[j]-x[i]);
f[i] = M_PI + 2.0 * sum1;
}
}
void SolveRK4(const int n, double h,vector<double> x,
vector<vector<double>>& x_vec,
int nstep, vector<double>& times)
{
/*
* times : vector contains the time 0 : t_final step dt
* x_vec : [nstep by N] 2 Dimensional vector
* dim1 : defined as typedef vector<double> dim1
*/
times[0] = 0.0;
dim1 y(n);
dim1 f1(n),f2(n),f3(n),f4(n);
double half_h = 0.5 * h;
double h_sixth = h/6.0;
// x_vec[nstepxN]
for (int i=0; i<n; i++)
x_vec[0][i] = x[i];
for (int k=1; k<nstep; k++){
times[k] = k*h;
xprsys(n,x,f1);
#pragma omp parallel for
for(int i=0; i<n; i++)
y[i] = x[i] + half_h * f1[i];
#pragma omp master
xprsys(n,y,f2);
#pragma omp barrier
#pragma omp for
for(int i=0; i<n; i++)
y[i] = x[i] + half_h * f2[i];
#pragma omp master
xprsys(n,y,f3);
#pragma omp barrier
#pragma omp for
for(int i=0; i<n; i++)
y[i] = x[i] + h * f3[i];
#pragma omp master
xprsys(n,y,f4);
#pragma omp barrier
#pragma omp for
for(int j=0; j<n; j++) {
x[j] = x[j] + h_sixth * (f1[j] + f4[j] + 2.0 * (f2[j] + f3[j]));
x_vec[k][j] = x[j];
}
}
}
n
? If it is small then the overhead from the threads will just make it slower. \$\endgroup\$ – Emily L. Jul 6 '17 at 8:58n
about 1000. \$\endgroup\$ – Abolfazl Jul 6 '17 at 9:44#include
lines, and amain()
that shows how to call your function. It's not mandatory, but it really helps! Also, you might want to indicate which OpenMP version your code requires - I think it needs to be v4 or later forreduction(+:sum1)
, doesn't it? \$\endgroup\$ – Toby Speight Jul 6 '17 at 14:30n
, so I thought it worth mentioning. \$\endgroup\$ – user1118321 Jul 8 '17 at 0:33