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.
sqrt()
are in SSE/AVX, but they take much longer than addition and multilication. Also, they are poorly vectorized on Sandy Bridge AVX, taking twice as long as the SSE instruction (with half the width). For example, double precision AVX (4 doubles wide) can do a packed multiply and packed add every cycle (assuming no dependencies or stalls on memory) which is 8 flops per cycle. The divide takes between 20 and 44 cycles to do those "4 flops". $\endgroup$ – Jed Brown Jan 24 '12 at 20:57