Runge-Kutta Fourth Order in C

Question

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

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.