FLOP counting for library functions

Question

When evaluating the number of FLOPs in a simple function, one can often just go down the expression tallying basic arithmetic operators. However, in the case of mathematical statements involving even division, one cannot do this and expect to be able to compare with FLOP counts from functions with only additions and multiplications. The situation is even worse when the operation is implemented in a library. Therefore, it is imperative to have some reasonable notion of the performance of the special functions.

By special functions, we mean things like:

• exp()
• sqrt()
• sin/cos/tan()

which are usually provided by system libraries.

Determining the complexity of these is confounded even further by the fact that many of them are adaptive and have input-dependent complexity. For example, numerically stable implementations of exp() often adaptively rescale and use look-ups. My initial impression here is that the best one may do in this case is ascertain the average behavior of the functions.

This entire discussion is, of course, highly dependent on the architecture. For this discussion we can restrict ourselves to traditional general purpose architectures and exclude those with special function units (GPUs, etc.)

One can find fairly simple attempts to standardize these for particular architectures for the sake of system vs. system comparison, but this is not acceptable if one cares about method vs. method performance. Which methodologies for determining the FLOP complexity of these functions are considered to be acceptable? Are there any major pitfalls?

Answer 1

It sounds like you want a way to evaluate how FPU-bound your code is, or how effectively you are using the FPU, rather than to count the number of flops according to same anachronistic definition of a "flop". In other words, you want a metric that reaches the same peak if every floating point unit is running at full capacity every cycle. Let's look at an Intel Sandy Bridge to see how this might shake out.

Hardware-supported floating point operations

This chip supports AVX instructions, so registers are 32 bytes long (holding 4 doubles). The superscalar architecture allows instructions to overlap, with most arithmetic instructions taking a few cycles to complete, even though a new instruction might be able to start on the next cycle. These semantics are usually abbreviated by writing latency/inverse throughput, a value of 5/2 would mean that the instruction takes 5 cycles to complete, but you can start a new instruction every other cycle (assuming that the operands are available, so no data dependence and not waiting for memory).

There are three floating point arithmetic units per core, but the third is not relevant to our discussion, we'll call the relevant two the A and M units because their primary functions are addition and multiplication. Example instructions (see Agner Fog's tables)

• vaddpd: packed addition, occupying unit A for 1 cycle, latency/inverse throughput is 3/1
• vmulpd: packed multiplication, unit M, 5/1
• vmaxpd: packed select pairwise maximum, unit A, 3/1
• vdivpd: packed divide, unit M (and some A), 21/20 to 45/44 depending on input
• vsqrtpd: packed square root, some A and M, 21/21 to 43/43 depending on input
• vrsqrtps: packed low-accuracy reciprocal square root for single precision input (8 floats)

The precise semantics for what can overlap with vdivpd and vsqrtpd are apparently subtle and AFAIK, not documented anywhere. In most uses, I think there is little possibility for overlap, though the wording in the manual suggests that multiple threads may offer more possibility for overlap in this instruction. We can hit peak flops if we start a vaddpd and vmulpd on every cycle, for a total of 8 flops per cycle. Dense matrix-matrix multiply (dgemm) can get reasonably close to this peak.

When counting flops for special instructions, I would look at how much of the FPU is occupied. Suppose for argument that in your range of input, vdivpd took an average of 24 cycles to complete, fully occupying unit M, but addition could (if it was available) be executed concurrently for half the cycles. The FPU is capable of performing 24 packed multiplies and 24 packed additions during those cycles (perfectly interleaved vaddpd and vmulpd), but with a vdivpd, the best we can do is 12 additional packed adds. If we suppose that the best possible way to do division is to use the hardware (reasonable), we might count the vdivpd as 36 packed "flops", indicating that we should count each scalar divide as 36 "flops".

With reciprocal square root, it is sometimes possible to beat the hardware, especially if full accuracy is not needed, or if the range of input is narrow. As mentioned above, the vrsqrtps instruction is very inexpensive, so (if in single precision) you could do one vrsqrtps followed by one or two Newton iterations to clean up. These Newton iterations are just

y *= (3 - x*y*y)*0.5;

If many of these operations need to be performed, this can be significantly faster than naive evaluation of y = 1/sqrt(x). Prior to the availability of hardware approximate reciprocal square root, some performance-sensitive code used infamous integer operations to find an initial guess for the Newton iteration.

Library-provided math functions

We can apply a similar heuristic to library-provided math functions. You can profile to determine the number of SSE instructions, but as we have discussed, that is not the whole story and a program that spends all its time evaluating special functions might not appear to get close to peak, which might be true, but isn't useful for telling you that all the time is spent out of your control on the FPU.

I suggest using a good vector math library as the baseline (e.g. Intel's VML, part of MKL). Measure the number of cycles for each call and multiply by peak achievable flops over that number of cycles. So if a packed exponential takes 50 cycles to evaluate, count it as 100 flops times the register width. Unfortunately, vector math libraries are sometimes hard to call and don't have all the special functions, so you might end up doing scalar math, in which case you would count our hypothetical scalar exponential as 100 flops (even though it probably still takes 50 cycles, so you would only be getting 25% of "peak" if all the time is spent evaluating these exponentials).

As others have mentioned, you can count cycles and hardware event counters using PAPI or various interfaces. For simple cycle counting, you can read the cycle counter directly using the rdtsc instruction with a snippet of inline assembly.

Answer 2

You could count them on real systems using PAPI, which grants access to hardware counters, and simple test programs. My favorite PAPI interface/wrapper is IPM (Integrated Performance Monitor) but other solutions exist (TAU, for example). This should give a fairly stable method-to-method comparison.

Answer 3

I am going to answer this question as if you asked:

"How do I analytically compare or predict the performance of algorithms that heavily rely on special functions, instead of the traditional multiply-add-carry FLOP counts that come from numerical linear algebra"

I agree with your first premise, that the performance of many special functions is architecture-dependent, and that although you can usually treat each of these functions as having constant cost, the size of the constant will vary, even between two processors from the same company but with different architectures (see Agner Fog's instruction timing table for reference).

I disagree, though, that the focus of the comparison should be on the costs of the individual floating point operations. I think that counting FLOPs is to some extent still useful, but that there are several much more important considerations that may make the cost of special functions less relevant when comparing two potential algorithms, and these should be explicitly examined first before going to a comparison of floating-point operations:

1. Scalability - Algorithms featuring tasks that can be implemented efficiently on parallel architectures will dominate the scientific computing arena for the foreseeable future. An algorithm with a better "scalability", be it through lower communication, less need for synchronization, or better natural load balance, may employ more slower special functions and therefore be slower for small numbers of processes, but will eventually catch up as the number of processors is increased.

2. Temporal Locality of Reference - Does the algorithm reuse data between tasks, allowing the processor to avoid unnecessary memory traffic? Each level of the memory hierarchy that an algorithm traverses adds another order of magnitude cost (roughly) to each memory access. As a result, an algorithm with high density of special operations will likely be significantly faster than an algorithm with the equivalent number of simple function operations over a larger region of memory.

3. Memory Footprint - This is strongly related to the previous points, but as computers grow larger and larger, the amount of memory per core is actually trending downward. There are two benefits to a small memory footprint. The first is that a small amount of program data will likely be able to fit completely within the processor cache. The second is that, for very large problems, an algorithm with a smaller memory footprint may be able to fit into processor memory, allowing problems to be solved that would otherwise exceed the capability of the computer.

Answer 4

Why bother counting flops? Just count cycles for every operation and you'll have something that is universal.