# Runge-Kutta Fourth Order in C

I've found that the Runge-Kutta (4th order) calculations in some software I wrote are the bottleneck. Is there anything obvious I can do to improve efficiency here?

Note that Compiler optimizations are ON.

I am profiling without optimizations though (not sure if there's a good way around this).

## Profiling results:

Over a run of 2560000 calls to matsuoka_calc_nextVal_RK The three functions contributing the most are:

matsuoka_rkStep - 7107ms (function body, not inc called functions)

matsuoka_calc_deriv - 7550ms

matsuoka_euler_step - 6040ms

#define POSPART(X)  (X > 0.0 ? X : 0.0)
#define NEGPART(X)  (X < 0.0 ? X : 0.0)

#define STEP_DIVISION 16.0
#define DIVISIONS_PER_CYCLE 2.0

typedef struct _matsuoka_internals {
// equation internal variables
double x1;
double x2;
double v1;
double v2;
}   matsuoka_internals;

// calculates the next value via Runge Kutta, via all intermediate steps
double matsuoka_calc_nextVal_RK(double in, double t1, double t2,
double c1, double c2, double b, double g,
matsuoka_internals *internals)
{
double step = 1.0 / STEP_DIVISION;
double t1recip = 1.0 / t1;
double t2recip = 1.0 / t2;

matsuoka_enforceStability(&g, &c1, &c2, b, t1, t2);

int i = (int)STEP_DIVISION;
while (i--) {
matsuoka_rkStep(in, t1recip, t2recip, c1, c2, b, g, internals, step);
}
matsuoka_fixNAN(internals);
return POSPART(internals->x1) - POSPART(internals->x2);
}

// single step in Runge Kutta calculation
void matsuoka_rkStep(double in, double t1recip, double t2recip,
double c1, double c2, double b, double g,
matsuoka_internals *internals, double step) {
matsuoka_internals k1, k2, k3, k4;

k1 = matsuoka_calc_deriv(in, t1recip, t2recip, c1, c2, b, g,
*internals);
k2 = matsuoka_calc_deriv(in, t1recip, t2recip, c1, c2, b, g,
matsuoka_eulerStep(*internals, k1, step*0.5));
k3 = matsuoka_calc_deriv(in, t1recip, t2recip, c1, c2, b, g,
matsuoka_eulerStep(*internals, k2, step*0.5));
k4 = matsuoka_calc_deriv(in, t1recip, t2recip, c1, c2, b, g,
matsuoka_eulerStep(*internals, k3, step));

internals->x1 = internals->x1 + (k1.x1 + 2 * (k2.x1 + k3.x1) + k4.x1)*(1.0 / 6.0)*step;
internals->x2 = internals->x2 + (k1.x2 + 2 * (k2.x2 + k3.x2) + k4.x2)*(1.0 / 6.0)*step;
internals->v1 = internals->v1 + (k1.v1 + 2 * (k2.v1 + k3.v1) + k4.v1)*(1.0 / 6.0)*step;
internals->v2 = internals->v2 + (k1.v2 + 2 * (k2.v2 + k3.v2) + k4.v2)*(1.0 / 6.0)*step;

}

// calculates derivative internals for matsuoka oscillator
matsuoka_internals matsuoka_calc_deriv(double in, double t1recip, double t2recip,
double c1, double c2, double b, double g,
matsuoka_internals internals)
{
double posX1 = POSPART(internals.x1);
double posX2 = POSPART(internals.x2);
double posIn = POSPART(in);
double negIn = NEGPART(in);

matsuoka_internals deriv;
deriv.x1 = (c1 - internals.x1 - (b*(internals.v1)) - (g*posX2) - posIn) * t1recip;
deriv.x2 = (c2 - internals.x2 - (b*(internals.v2)) - (g*posX1) + negIn) * t1recip;
deriv.v1 = (posX1 - internals.v1) * t2recip;
deriv.v2 = (posX2 - internals.v2) * t2recip;

return deriv;
}

// single step in euler's method calculation
matsuoka_internals matsuoka_eulerStep(matsuoka_internals init, matsuoka_internals deriv, double step) {
matsuoka_internals newVal;
newVal.x1 = init.x1 + (deriv.x1*step);
newVal.x2 = init.x2 + (deriv.x2*step);
newVal.v1 = init.v1 + (deriv.v1*step);
newVal.v2 = init.v2 + (deriv.v2*step);
return newVal;
}

// enforces parameters within stable range (see K. Matsuoka, Biological cybernetics, 1985.)
void matsuoka_enforceStability(double *g, double *c1, double *c2,
double b, double t1, double t2) {
double gMin = 1.0 + (t1 / t2);
*g = *g < gMin ? gMin : *g;

double cMin = *g / (1.0 + b);
*c1 = *c1 < cMin ? cMin : *c1;
*c2 = *c2 < cMin ? cMin : *c2;
}

// protects against denormal
void matsuoka_fixNAN(matsuoka_internals *internals) {
if (!isnormal(internals->x1) || !isnormal(internals->x2)
|| !isnormal(internals->v1) || !isnormal(internals->v2)) {

internals->x1 = 0, internals->x2 = 0, internals->v1 = 0, internals->v2 = 0;
}
}

• Did you measure performance? If not, please do that first. There is a saying about optimization: "You do not care about the performance of code you don't measure." – Ben Steffan Mar 4 '18 at 17:38
• "I'm finding that the Runge Kutta calcs in some software I wrote are the bottleneck." I found that by profiling. – DanBennett Mar 4 '18 at 17:44
• @DanBenett I probably didn't express myself very well. I meant that you should measure and provide us with your measurements. Also, you are providing a lot of code for a bottleneck. Where is the bottleneck exactly? What function is slow? Having precise data would make this review a lot easier for us (and ultimately enable us to help you better). – Ben Steffan Mar 4 '18 at 17:47
• Ah, more likely it was me who misunderstood. Thanks. rkStep function body and its two called functions, calc_deriv and eulerStep are all (roughly) equally contributing (in the debug version with no compiler optimisations off - I'm a bit stuck when it comes to profiling meaningfully with optimisations on). I'll fill in some more details later tonight. – DanBennett Mar 4 '18 at 18:55
• Profiling with a debug build is just about useless. Using VS17 and /Ox, inspection of the generated assembly shows rkStep being completely inlined, and MMX being used in various places (which works on 2 double calculations at a time). – 1201ProgramAlarm Mar 4 '18 at 19:50

Not a thorough answer but some minor contributions to speed up:

1. Get rid of divisions: at each step there are four 1/6 divisions, these can be replaced with multiplication with a constant (one_sixth = 1/6.0).
2. Enforce stability function can better receive t2recip, to avoid another division.
3. 1/6 * step can be combined to become one_sixth_times_step, to avoid redoing the calculations 4 times. (Saves 3 multiplications at each step).

A few more multiplications could possibly be avoided by merging some definitions (step*0.5 can be replaced with half step).

Calling different functions can look better but uniting them might get rid of some overhead. Just a thought.

Finally, calling inline assembly (maybe as lambda functions) can be worth the effort.

• Thanks very much Gürkan, I'll try these suggestions - I should have spotted at least a couple of these myself really. I had assumed the functions will be inlined by the compiler (my compiler settings are prioritising speed over size and permitting any inlining) but I reach the limits of my understanding quite quickly on that point. – DanBennett Mar 4 '18 at 18:58
• Are you running optimization level 2 or 3? I'm not experienced with details of optimizations, maybe they do part of what I've suggested. – Gürkan Çetin Mar 4 '18 at 19:01
• This is in VisualStudio - Ox which I think is roughly equivalent to gcc 3, with "inline function expansion" turned on in addition. I have no idea how to tell if it's doing what you're suggesting either - I just assumed it would! But I'm not at all sure if I'm correct on that. If I still need improvements after trying your other suggestions I'll try uniting the functions as you suggest. – DanBennett Mar 4 '18 at 19:13
• Here’s a video I came across regarding compiler optimizations. youtu.be/w5Z4JlMJ1VQ – Gürkan Çetin Mar 4 '18 at 19:40
• Your first point is moot, as that constant computation will be done at compile time. – 1201ProgramAlarm Mar 4 '18 at 19:51