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

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);
    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,
    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;
  • \$\begingroup\$ 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." \$\endgroup\$ Mar 4, 2018 at 17:38
  • \$\begingroup\$ "I'm finding that the Runge Kutta calcs in some software I wrote are the bottleneck." I found that by profiling. \$\endgroup\$
    – DanBennett
    Mar 4, 2018 at 17:44
  • 1
    \$\begingroup\$ @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). \$\endgroup\$ Mar 4, 2018 at 17:47
  • \$\begingroup\$ 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. \$\endgroup\$
    – DanBennett
    Mar 4, 2018 at 18:55
  • 2
    \$\begingroup\$ 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). \$\endgroup\$ Mar 4, 2018 at 19:50

1 Answer 1


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.

  • \$\begingroup\$ 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. \$\endgroup\$
    – DanBennett
    Mar 4, 2018 at 18:58
  • \$\begingroup\$ 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. \$\endgroup\$ Mar 4, 2018 at 19:01
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – DanBennett
    Mar 4, 2018 at 19:13
  • \$\begingroup\$ Here’s a video I came across regarding compiler optimizations. youtu.be/w5Z4JlMJ1VQ \$\endgroup\$ Mar 4, 2018 at 19:40
  • 1
    \$\begingroup\$ Your first point is moot, as that constant computation will be done at compile time. \$\endgroup\$ Mar 4, 2018 at 19:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.