Efficiently simulating a control sequence

Question

I have a control problem that I am simulating in C. The premise is that I have some continuous-time signal s(t) and I am given some length-150 array c = [c(0), c(dt), c(2dt), ...] of controls. The array is a discretized version of a continuous-time control c(t), but I am not actually given access to the full c(t), only to the array c with a fixed value of dt.

Given c and dt, I create an array s = [s(0), s(dt), s(2dt), ...]. I then need to compute the following function:

#include <math.h> #include <complex.h> #define SEQUENCE_LENGTH 150 double f(double s[], double c[], double dt) { // s and c are both length SEQUENCE_LENGTH arrays double complex alpha = 1. / sqrt(2) + 0.0 * I; double complex beta = 1. / sqrt(2) + 0.0 * I; double complex alphanew; double complex betanew; double x, y, norm, co, si; for (int j = 0; j < SEQUENCE_LENGTH; j++) { x = c[j]; y = s[j]; norm = sqrt(x*x + y*y); co = cos(norm * dt); si = sin(norm * dt) / norm; alphanew = co * alpha - I * si * (y * alpha + x * beta); betanew = co * beta - I * si * (x * alpha - y * beta); alpha = alphanew; beta = betanew; } return (1 / 2.) * pow(cabs(alpha + beta), 2); } 

This function is called on the order of a hundred million times for various different c and s and it is the largest bottleneck in my simulation, so I really need to make it as fast as possible. Does anybody have any suggestions/tips/tricks to speed up this code?

Various thoughts. Currently, I compile with gcc with the flags -Wall -Wextra -Werror -O3 -ffast-math.

• Because I know SEQUENCE_LENGTH is 150, manually unrolling part of the loop would be straightforward, but I'm sure the compiler already does this.
• I could perhaps vectorize some of the calculations in the loop, though again I assume the compiler is already doing this. Nonetheless I would still be interested in seeing an explicit vectorization of my code (for academic purposes), but I am not very familiar with vectorization.
• Is there any way to parallelize this code, eg with #pragma amp parallel? I don't see an obvious way.

Answer 1

I really need to make it as fast as possible.

Does anybody have any suggestions/tips/tricks to speed up this code?

Make a test harness

Form test code that rates the performance and precision correctness.

There simply are too many factors for a general fastest. We are currently left with faster suggestions. Yet with a use case specific test harness, we can reliably improve things (simpler trig/power functions, etc.)

float vs. double

Often float functions are faster.

Use float objects and float functions like sqrtf(), sinf(), powf(), ...

Do you really need double precision?

What is the precision goal? Optimization level

In requesting speed, one needs to identify the required precision as various speed optimizations can trade off correctness for time simply by enabling higher floating point compile-time optimizations levels.

Was ffast-math really OK to use with an unstated precision goal?

Conversely hypot(x, y) is often more precise than norm = sqrt(x*x + y*y);. Is that OK?. OP has left no test harness and OP simply has not specified.

Range of data

State the range.

When double s[] and double c[] use a subset of the double range, additional optimizations are possible. @user555045 This should be part of the test harness.

Use restrict and const

Assuming s[] and c[] do not overlap, use restrict.

This allows the compiler to make additional optimizations.

// double f(double s[], double c[], double dt) { double f(const double * restrict s, const double * restrict c, double dt) { 

---

Avoid loss of precision when no speed is achieved

Consider 1. / sqrt(2). Why isn't it sqrt(0.5)?

Review the effect:

 printf("%a\n", 1. / sqrt(2)); printf("%a\n", sqrt(0.5)); 

Output is not the same. If a speed difference exist, sqrt(0.5) is likely faster.

0x1.6a09e667f3bccp-1 0x1.6a09e667f3bcdp-1 

Of course we could use

// double complex alpha = 1. / sqrt(2) + 0.0 * I; double complex alpha = 0.70710678118654752440084436210485 + 0.0 * I; 

---

The premise is that I have some continuous-time signal s(t) ...
c = [c(0), c(dt), c(2dt), ...]

Architecture re-write

Consider performing the calculation as the data comes in doing a little bit of the long calculation during each dt.

For real applications, I have found this to be the fastest: distributing the calculation rather than do it all after all the samples are acquired.

If dt is sufficient for worse case 1 iteration of the for (int j = 0; j < SEQUENCE_LENGTH; j++) {, then it should be easy to implement.

Else code might use 2 threads: one to acquire, one to compute.

---

Watch out for pitfalls

sin(norm * dt) / norm; is perhaps a big problem when norm == 0.0as it is 0.0/0.0.

Consider safer code calculations approaches and avoid code that can crash quickly.

Good luck.