Runge-Kutta fourth order integration

Question

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

Answer 1

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.

I did so here https://github.com/miscco/Neural_Network_TC/blob/master/Iterate_ODE.h

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.

Answer 2

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