Microbenchmarking of algorithm

Question

I want to compare the speed of functions in C. This is my approach.

I measure the time of each implementation by making the calculation on the input array a certain number of times and adding up the time of each iteration using the monotonic clock and clock_gettime.

    totalTime = 0.0;
    for (int i = 0; i < ITERATIONS; ++i) {
        startTime = (float)clock()/CLOCKS_PER_SEC;

        doFun(input, output);
        endTime = (float)clock()/CLOCKS_PER_SEC;

        timeElapsed = endTime - startTime;

        totalTime += timeElapsed;
    }
    printf("Time:       %f \n", totalTime);

-O0
Native:       26.296245 
Quick:   34.653388 
SSE:        21.701504 

-O1
Native:       11.036241 
Quick:   7.885706 
SSE:        2.043576 

-O2
Native:       0.000776 
Quick:   0.000795 
SSE:        0.000761 

I am really unsure if this is even the right approach to compare the algorithm and alternatives.

Especially for -O0 I don't understand how the algorithm is slower than the native function although it should be up to 4 times faster and the SSE variant again ideally 4 times faster than that.

At what optimization level should I compare the algorithms to have meaningful results and why do I not get the expected times at -O0?

Answer 1

Benchmarking like this isn't really all that helpful. Your benchmark results for -O2 are very close together - the reason is that you're benchmarking your memory bandwidth, and not your CPU performance. Benchmarking -O0 and -O1 doesn't make a great deal of sense tbh. Also ensure you have -ffast-math enabled.

Anyhow, using boring old std::sqrt:

void test_sqrt(float* __restrict a, const float* __restrict b) {
    for(int i = 0; i < 256; ++i) {
        a[i] = 1.0f / sqrtf(b[i]);
    }
}

It unrolls the loop so that each iteration 32 floats are processed (I've only included the asm for 8 of those floats).

        // load
        vmovups ymm2, ymmword ptr [rsi + 4*rax]

        // 1/sqrt
        vrsqrtps        ymm3, ymm2 // latency 4, TP = 1

        // refine
        vmulps  ymm2, ymm2, ymm3 // latency 4, TP = 2
        vfmadd213ps     ymm2, ymm3, ymm0 // latency 4, TP = 2
        vmulps  ymm3, ymm3, ymm1 // latency 4, TP = 2
        vmulps  ymm2, ymm3, ymm2 // latency 4, TP = 2

        // output
        vmovups ymmword ptr [rdi + 4*rax], ymm2

Things to note:

1. The compiler has removed the division, and is using rsqrt + a couple of newton-raphson iterations to refine the result.
2. The compiler has vectorised the code to be 8 wide (so an improvement over your SSE version).

The quake3 rsqrt (original version):

        vmovdqu ymm3, ymmword ptr [rsi + 4*rcx]
    
        vpsrad  ymm7, ymm3, 1 // latency 4, TP = 1
        vpsubd  ymm7, ymm0, ymm7  // latency 1, TP = 3
        vmulps  ymm3, ymm3, ymm1  // latency 4, TP = 2
        vmulps  ymm11, ymm7, ymm7  // latency 4, TP = 2
        vfnmadd213ps    ymm11, ymm3, ymm2   // latency 4, TP = 2
        vmulps  ymm3, ymm11, ymm7  // latency 4, TP = 2

        vmovups ymmword ptr [rdx + 4*rcx], ymm3

Again, the compiler has vectorised this code to use 8 wide AVX ops.

The two versions from a performance POV are almost identical (the Quake3 version has 1 additional instruction).

The big difference though is wrt to the accuracy. The quake 3 version is not brilliant on that front.

If you want absolute performance, use vrsqrtps. It's one instruction, and it has more accuracy than the Quake3 routine.

void test_rsqrt(float* __restrict a, const float* __restrict b) {
    for(int i = 0; i < 256; i += 8) {
        __m256 A = _mm256_load_ps(b + i);
        __m256 B = _mm256_rsqrt_ps(A);
        _mm256_store_ps(a + i, B);
    }
}