# Runge-Kutta fourth order integration

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];
}
}
}

• How big is n? If it is small then the overhead from the threads will just make it slower. – Emily L. Jul 6 '17 at 8:58
• I know, but I have no better solution. Once I test creating the threads before the time loop, and using proper barriers at the end of each time step, It ruined the speed up, therefore creating the threads inside the time loop was better for large n about 1000. – Abolfazl Jul 6 '17 at 9:44
• You might have to try and balance the benefits of openmp vs using vectorization (intrinsics) your compiler might have auto vectorized some of the statements. Our you could use the Eigen math library that will perform vectorization on matrix and vector expressions. – Harald Scheirich Jul 6 '17 at 14:20
• Welcome to Code Review! You'll receive better reviews the more complete the code you show. For example, I recommend that you show the necessary #include lines, and a main() 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 for reduction(+:sum1), doesn't it? – Toby Speight Jul 6 '17 at 14:30
• @Abolfazl No, but it eliminates a possible error - sending the wrong n, so I thought it worth mentioning. – user1118321 Jul 8 '17 at 0:33

## 2 Answers

I would really suggest that you better structure your code. Runge Kutta methods are quite straightforward, so you should have a function that iterates a single step in parallel for every variable.

Furthermore, you can slightly deviate from the standard Runge Kutta scheme if you directly store $x + h_if_i$ in the $i$th component of that vector. You only have to adopt the summation at the end to account for that.

static void add_RK(std::vector<double>& var) {
var[0] = (-3*var[0] + 2*var[1] + 4*var[2] + 2*var[3] + var[4])/6;
}


A minimal improvement that I would do is changing this

for (int i=0; i<n; i++)
x_vec[0][i] = x[i];


To

x_vec[0] = x;


Where you utilize the underlying copy/move constructor, which might be more efficient depending on the compiler.

• It's very nice structured code, and I'm happy that is about neuroscience however It's still difficult for me to track the functions and classes. – Abolfazl Jul 6 '17 at 14:40

To follow the @miscco advice I wrote this version which uses valarray and has good performance. I only used one omp loop I did it here (link to gitlab repository)

void euler_integrator (StateVec &y, DerivFunc dydt, double dt) {
y += dydt (y) * dt;
}

/*------------------------------------------------------------*/
void runge_kutta4_integrator (StateVec &y, DerivFunc dydt, double dt) {
auto k1 = dydt (y);
auto k2 = dydt (y + dt*k1/2.);
auto k3 = dydt (y + dt*k2/2.);
auto k4 = dydt (y + dt*k3);
y += (k1 + 2.*k2 + 2.*k3 + k4) * dt/6;
}

/*------------------------------------------------------------*/
void integrate (Integrator integrator, DerivFunc dydt, OutputFunc output_func,
std::ofstream &output, StateVec y, int num_steps, double dt) {
for (auto step = 0; step < num_steps; ++step) {
output_func (step*dt, y, output);
integrator (y, dydt, dt);
}
}
/*------------------------------------------------------------*/
StateVec xprsys (const StateVec &x)
{
double sumx=0;
int N = PARN;
StateVec f;
f = StateVec(N);
#pragma omp parallel for reduction(+:sumx)
for (int i=0; i<N; i++){
sumx = 0;
for(int j=0; j<N; j++)
sumx += sin(x[j]-x[i]);
f[i] = PARw + PARka_over_N * sumx;
}
return f;
}