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;
}
}
/Ox
, inspection of the generated assembly showsrkStep
being completely inlined, and MMX being used in various places (which works on 2 double calculations at a time). \$\endgroup\$ – 1201ProgramAlarm Mar 4 '18 at 19:50