C - have a simple loop that does arithmetic calculation; profiler shows this is a bottleneck. How to speed it up?

Question

My first post here. A great site and resource.

I did search a bit and looked at the questions with similar titles, but couldn't find something specifically about this.

I'm trying to remove any redundancy and bloat from a C astronomical calculation library that my C++ program uses. I ran a simple profiler (VerySleepy).

Here is the code that the profiler showed as using the most time (aside from C library functions sprintf, etc.):

double swi_echeb(const double x, const double* const coef, const int ncf)
{
    int j = ncf - 1;
    double x2, br, brp2, brpp;
    x2 = x * 2.;
    br = 0.;
    brp2 = 0.;  /* dummy assign to silence gcc warning */
    brpp = 0.;

    for (; j >= 0; --j) {                 // <-- 0.39s
        brp2 = brpp;                      // <-- 0.01s
        brpp = br;                        // <-- 0.32s
        br = x2 * brpp - brp2 + coef[j];  // <-- 3.49s ***
    }                                     // <-- 0.14s

    return (br - brp2) * .5;              // <-- 0.06s
}                                         // <-- 0.05s

This particular function is deeply embedded within others and the main "kick-off" function that my program calls is called thousands of times.

You can see the standout statement with 3.49s as much higher than all the other statement times. I know there are ways to speed C arithmetic with using multiplication over division when possible. But I don't know much more than that.

Like:

1. Would it be better to split this statement up into smaller pieces?:

   br = x2 * brpp;

   br -= brp2;

   br += coef[j];

Any other ideas or critiques. I did not write this code, though I did add the const to the function parameters as I love const-correctness.

I've never tried using registers or other fancy tricks to speed things up before. Anyone think something like that can work here?

I know people will say, "Try it!" So I will, and will update what I get if it helps anyone with similar arithmetic questions.

EDIT: Posting Results I've Tested From Suggestions

In order from fastest to slowest, here are what I've found so far. Profiler is VerySleepy. Compiler is Visual Studio 2008 Pro Ed. Compile options for both library and my application are:

Debug, C7 format, /O2 /Ob2 /Oi /Ot /Oy /GT /GL /GF /FD /MTd /GS- /Gy /fp:fast /FAs

The following is Andrew's suggestion about doing "4 iterations per loop". It was the fastest so far.

TOTAL TIME spent in function (times from the other statements in the function are not shown here) = 2.08 seconds

for (; index >= 3; index -= 4) {                    // 0.02s
    brp2    = brpp;
    brpp    = br;                                   // 0.02s
    br      = x2 * brpp - brp2 + coef[index];       // 0.25s
    brp2    = brpp;
    brpp    = br;                                   // 0.13s
    br      = x2 * brpp - brp2 + coef[index - 1];   // 0.33s
    brp2    = brpp;
    brpp    = br;                                   // 0.13s
    br      = x2 * brpp - brp2 + coef[index - 2];   // 0.34s
    brp2    = brpp;
    brpp    = br;                                   // 0.14s
    br      = x2 * brpp - brp2 + coef[index - 3];   // 0.42s
}

for (; index >= 0; --index) {                 // 0.03s
    brp2    = brpp;                           // 0.03s
    brpp    = br;
    br      = x2 * brpp - brp2 + coef[index]; // 0.11s
}

The next fastest was the original unaltered code, with a total time of 2.39 seconds inside the function, again including the statements outside the loop. Note that this is less than my original post. My original post was unoptimized code, but since everyone suggested it, all of my tests were subsequently as optimized as I could get in VS08:

for (j = ncf - 1; j >= 0; j--) {      // 0.02s
    brp2 = brpp;                      // 0.03s
    brpp = br;                        // 0.07s
    br = x2 * brpp - brp2 + coef[j];  // 2.14s
}

After this original code, the next fastest was Drew's idea of setting the pointer in advance and using that. Total time spent inside function was 2.49 seconds, including times from statements outside loop:

for (; index >= coef; --index) {         // 0.01s
    brp2    = brpp;
    brpp    = br;                        // 0.06s
    br      = x2 * brpp - brp2 + *index; // 2.24s
}

I also tried a mix of both Andrew's loop unrolling and Drew's pointer usage, but that took 2.39 seconds, the same as the unaltered code.

Based on the results, the loop-unrolling is the way to go so far for my usage.

Answer 1

First thing I would try would be to iterate in steps of 4, ie: j+=4 (or in your case, j -=4) and semi-unroll the loop. The reason for this is it will help the compiler to make SSE optimisations and to batch memory access from main memory to cache. Just be aware that you will have to cater for the last few elements in case the loop count is not divisible by 4. For example:

// Disclaimer: I have not tested this code!
for (; j >= 3; j -= 4) {              
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef[j]; 
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef[j-1]; 
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef[j-2]; 
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef[j-3]; 
}                          
// if (j % 4) != 0 before the loop operation, 
// handle 1, 2 or 3 remaining elements here

Second thing I would try would be to preload coeff[j] into a register immediate prior to the calculation. The reason for this is floating point calculations are pipelined, meaning that a memory access in the wrong place can have adverse effects on performance. The calculation itself can be very fast but might take 14 instructions just to queue up the data from cache into the FPU. Add to that an access from main memory it can get even worse. For instance, try this (could also be tried with and without the -=4 unrolling)

// Disclaimer: I have not tested this code!
register double coef1, coef2, coef3, ceof4;
for (; j >= 3; j -= 4) {           
    coef1 = coef[j];    // Preloads the 4 sequential coeffs from 
    coef2 = coef[j-1];  // main memory to cache (if available)
    coef3 = coef[j-2];  
    coef4 = coef[j-3];  
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef1; 
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef2; 
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef3; 
    brp2 = brpp;                      
    brpp = br;                        
    br = x2 * brpp - brp2 + coef4; 
} 

In this case the variables double x2, br, brp2, brpp, coef1, coef2, coef3, coef4 should be registers if at all possible.

Finally, using the above, can you apply SSE/SSE2 optimisation to it? Make sure this is enabled in the GCC compiler (I'm used to VS so the equivalent would be Release mode on, debug symbols off, optimization on, SSE2 on) and benchmark your code without the debugger attached. This alone can have a dramatic affect on performance.

Let us know the results. Performance tuning is trial and error!

Best of luck,

Answer 2

This appears to be a cache issue, not an arithmetic one.

for (; j >= 0; --j) {
    ...
    ... coef[j];
}

You're accessing an array here, and you are decrementing an index to do so. This action can really disrupt the cache-friendly locality inherent in a simple loop.

Is it possible to count forward? Ie,

for (int i = 0; i <= j; i++) {
    ...
    ... coef[i];
}

Would your calculation be valid?

Answer 3

I've never tried using registers or other fancy tricks to speed things up before. Anyone think something like that can work here?

There's a very easy register trick that anyone can do. Build the project for a recent CPU. Does this code need to run on a computer from 1995? 2000? 2005? If the program can count on a newer CPU, it can count on having more registers at its disposal.

Also, The integer indexing is unnecessary. You could instead make j a pointer directly to the double of interest. This may make a difference if your optimizing compiler isn't already doing it.

double swi_echeb(const double x, const double* const coef, const int ncf)
{
    const double *j = &coef[ncf - 1];
    // (stuff...)

    while (true) { 
        // (stuff...)

        br = x2 * brpp - brp2 + *j;
        if ( j == coef )
            break;
        --j;
    }  
}

Answer 4

The main 'problem' with that code is that you have a critical path along br. You can not begin to calculate the next iteration before you have completely finished the previous one. This also prohibits the of vector instructions: There is nothing to vectorize.

I have the impression that the number coefficients is always rather (single digit?) and the runtime stems from the amount of calls to that function.

One way to mitigate that is to calculate evaluate multiple polynomials at once. Of course this depends on a special layout of your data structures: The coefficients of a certain degree have to be in a linear array, so they can be loaded by a single vector instruction.

Answer 5

Well, unless there are some special issues -- like is your array of coefficients big enough that you could be swapping? -- you're probably pretty close.

• Andrew's notion of loop unrolling should be tried.
• Definitely make sure you have the optimization turned up with -O3

After that, you're going to need to look at the assembly language, or parallelism.

Answer 6

This operation is a slight variation of a prefix sum/scan. (It's a 3-op, 2-history scan). The key limiter to performance here is more than likely the serialization (of your math ops in the instruction pipe) caused by the cross loop dependencies, so serial loop unrolling is unlikely to help much here.

There are standard ways to parallelize prefix sums (see wikipedia), that could be used to accelerate this. Even with one thread, you would be able to greatly improve your efficiency by subdividing the coefficient array into 4 sub arrays, and computing the prefix some for each of them per loop iteration - the four streams of computation are independent, and will be properly pipelined by your hardware. Furthermore, since they are independent, you (or your compiler if you are luckly, but I doubt it) can utilize SSE or AVX on an x86 to process the array in parallel.

Once you have your four accumulated results (the results will likely be pairs since you have a 2-history prefix sum), you can combine them in the mathematically appropriate way for your sequence.

Answer 7

What are typical values for ncf? Main reason I ask is that you are iterating coef backwards. Non-sequential access is not make good use of the cache.

Answer 8

This is a case where the profiler finds what you would expect. Look at the code: you've got some loop setup ("oh, that runs once, it won't take up anything"), and a Loop. Inside the loop, you've got two assignments ("nope, those are cheap") and precisely one line of code that does a multiply, two additions, and an array reference.

You're not going to get this function to run a lot faster with micro-optimizations. The processor is actually spending its time doing the work that you want the function to do-- yeah, I know, shocking.

Your best bet is to go up a level or two. How can you reduce the number times this function is called? Is it getting called with the same parameters multiple times, so that you can cache the results? Are there places where you can use fewer coefficents, reducing the number of times the loop runs?

Answer 9

If you're going with Drew's suggestion of using a pointer rather than an array indexing, it might be necessary to restrict the pointer to see any benefit:

...
double *restrict index = coef + ncf - 1;
...
for (; index >= coef; --index) {
    brp2    = brpp;
    brpp    = br;
    br      = x2 * brpp - brp2 + *index;
}

This might help cache optimization, because the compiler can be certain that nobody will change the value that index is pointing to.

---

Also, I posted a similar problem last year, which got a number of great answers. Be sure to have a look.

Answer 10

Do you need the full precision of a double? Moving to float instead could save some time.

Although the optimizer should figure it out, adding an explicit register hint to your variable declarations couldn't hurt:

register double x2, br, brp2, brpp;

You could also try moving coef into a register:

register double* rc;
rc = coef;
. . .
br = x2 * brpp - brp2 + rc[j];

I don't know about this case in particular, but I have seen a surprising number of cases where the compiler doesn't correctly optimize compound expressions. You may be able to increase the odds of it doing the right thing by breaking it down to simple two-component expressions:

brp2 = brpp;     
brpp = br;
br = x2 * brpp;
br -= brp2;
br += rc[j];

You might have a look at the generated code, to see if there are any obvious inefficiencies.