Tips to Finetuning to increase the GFLOPS of a SIMD kernel

Question

Given an Image I with dimensions as {rows, columns}, for every pixel in the image, the following needs to be computed -

1. calc1 = image[i][j - 1] - image[i][j];
2. calc2 = image[i-1][j] - image[i][j];
3. calc3 = image[i][j-1] - 2 * image[i][j] + image[i][j + 1]
   
4. calc4 = image[i-1][j] - 2 * image[i][j] + image[i+1][j];
   
5. calc5 = image[i-1][j-1] - image[i-1][j] - image[i][j - 1] + image[i][j];

My current implementation is as follows -

inline void _myFunc(float* img_ptr,
        float* dx_ptr, float* dy_ptr, float* dxx_ptr, float* dyy_ptr, float* dxy_ptr, 
                        size_t rows, size_t cols){

    __m256 negTwos = _mm256_set1_ps(-2);

    for (size_t i = 1; i <rows - 1; i++)
    {
#pragma unroll
        for (size_t j = 1; j <cols - 1; j += 8)
        {
            auto idx = i * cols + j;
            auto idx_before = (i - 1) * cols + j;
            auto idx_after = (i + 1) * cols + j;

            _mm_prefetch(img_ptr + idx - 1, _MM_HINT_T0);
            _mm_prefetch(img_ptr + idx_before - 1, _MM_HINT_T0);
            _mm_prefetch(img_ptr + idx_after - 1, _MM_HINT_T0);

            __m256 __ij = _mm256_loadu_ps(img_ptr + idx);
            __m256 __ijn1 = _mm256_loadu_ps(img_ptr + idx - 1);
            __m256 __ijp1 = _mm256_loadu_ps(img_ptr + idx + 1);

            __m256 __in1jn1 = _mm256_loadu_ps(img_ptr + idx_before - 1);
            __m256 __in1j = _mm256_loadu_ps(img_ptr + idx_before);
            __m256 __ip1j = _mm256_loadu_ps(img_ptr + idx_after);

            // Read takes a total of 9 cycles 

            __m256 _dx = _mm256_sub_ps(__ijn1, __ij); 
            __m256 _dy = _mm256_sub_ps(__in1j, __ij); 
            //__m256 negTwo_times_ij = _mm256_mul_ps(negTwos, __ij);
            
            __m256 _dxx = _mm256_add_ps(__ijn1, _mm256_fmadd_ps(negTwos, __ij, __ijp1)); // Will FMA be Better ? 
            __m256 _dyy = _mm256_add_ps(__in1j, _mm256_fmadd_ps(negTwos, __ij, __ip1j)); 
            /*
            __m256 _dxx = _mm256_add_ps(__ijn1, _mm256_add_ps(negTwo_times_ij, __ijp1));
            __m256 _dyy = _mm256_add_ps(__in1j, _mm256_add_ps(negTwo_times_ij, __ip1j)); 
           */
            __m256 _dxy = _mm256_sub_ps(_mm256_sub_ps(__in1jn1, __in1j), _dx);

            _mm256_storeu_ps(dx_ptr  + idx, _dx);
            _mm256_storeu_ps(dy_ptr  + idx, _dy);
            _mm256_storeu_ps(dxx_ptr + idx, _dxx);
            _mm256_storeu_ps(dyy_ptr + idx, _dyy);
            _mm256_storeu_ps(dxy_ptr + idx, _dxy);

        }
        
    }
    
}

Currently, for a 9024 x 12032 image size, I obtain an average runtime of 220 ms and GFlops of 4.1, on a single thread (compiled with -O3 -mavx2 -march=native).

However, I am not satisfied with the current performance, and I suppose it can be better, specially by better computation and memory access.

Using a profiling tool,, I got to know the parts which take the maximum time(Among other things) -

__m256 _dxy = _mm256_sub_ps(_mm256_sub_ps(__in1jn1, __in1j), _dx);  :: 68.025ms
__m256 __ij = _mm256_loadu_ps(img_ptr + idx); :: 51.98 ms

All the Store Functions :: A cumulative of 42.065ms

__m256 _dxx = _mm256_add_ps(__ijn1, _mm256_add_ps(negTwo_times_ij, __ijp1)); :: 12.003ms
__m256 _dx = _mm256_sub_ps(__ijn1, __ij);  :: 7ms

A couple of observation (or rather questions)..

• Consider the lines __m256 _dxy = _mm256_sub_ps(_mm256_sub_ps(__in1jn1, __in1j), _dx); and __m256 _dxx = _mm256_add_ps(__ijn1, _mm256_add_ps(negTwo_times_ij, __ijp1));. Even though both essentially are the same in terms of cumulative latency, vaddps and vsubps both have a latency 4 each, but there is a stark difference in the time spent in these lines, that too almost 5 times more than the other. Why so ?

• All the store operations take ~ equal amount of time, the initial load operation __m256 __ij = _mm256_loadu_ps(img_ptr + idx); takes a majority of the time and other are relatively way faster. Is this due to the fact that the data would already be in the cache. If so, should not this be circumvented by _mm_prefetch and initial load should be fast as well. Is there any way load / store operations could be made faster.

I realise this kernel is both computationally as well as bandwidth bound(more computationally I suppose).

How could I make this faster.

I am running on Intel i7 9750h and the average frequency of CPU during the run is 4.1Ghz. I have used Intel Advisor tool to profile the code and obtain the runtimes.

TIA

Answer 1

Fine-grained profiling results normally need to be taken with a grain of salt. A very common situation is that (eg) a load takes a while, but the time is attributed to a later instruction instead. I don't think we should conclude that computing _dxy is the problem here. A way to verify that is taking that whole _dxy calculation out, and then seeing where the time goes. If it goes away, then yes, it really took that time. If the total time barely changes and just "moves" to something else, then it was probably misattributed (and likely still will be).

Prefetch

Prefetching has a significant latency, unless the thing being prefetched is already in L1 (and then why prefetch it). In the code here, data is prefetched and then immediately loaded, when the prefetch has no chance to be done yet. The address to prefetch should therefore be offset somewhat from the address that is being loaded. Given that the access pattern here is very linear, prefetching may not help much anyway, but with an offset it at least stands a chance.

Reordering to reduce cache misses

The very low flops seems to indicate cache misses (you could find out for real with a profiler). A row of the image (~36KB) does not fit in L1D cache (32KB), so I expect cache thrashing: each time the data is loaded into L1, it is used only once, then it gets booted from the cache, and it has to be loaded 2 more separate times for the other two rows that it is used in. Prefetching (with sufficient but not overly large offset) may alleviate the misses, but there is another solution: reorder the computation so that more of it happens when the data is in L1.

For example, instead of computing entire rows in one go, chop them into smaller pieces, maybe 1024 pixels across or so (size should be set by experiment). That way, for each time the data is loaded into L1, it can hopefully be used all 3 times that it will be used, because this way the data would be needed again before it gets kicked out due to reading too much data in the meantime.

Reordering can also be attained by unrolling.

Unrolling

Loop unrolling could help. Not to reduce loop overhead, but to group together computations that re-use the same data, so that we don't have to rely on it still being in L1 cache. The unrolling would mainly happen in the i direction, down the image, less so "across" the image.

Each extra set of 5 results that is computed by unrolling "downwards", only costs 3 more loads (if I counted them correctly), the other data was already loaded previously. That sounds good but it shouldn't be done too much: accessing too many rows in quick succession without iterating horizontally could lead to TLB thrashing. Unrolling "across" the image could help to reduce the amount of waste due to touching part of a cache line but not using the whole thing at once, but that mainly becomes relevant when the "down" unroll factor is very large which I don't think it should be. I expect that the "across" unroll factor should be quite small, maybe 1 (no unrolling) or 2, maybe up to 4 but sounds extreme to me, and that the "down" unroll factor should be medium, maybe up to 16, but you shouldn't trust these numbers too much and set them by experiment.

Answer 2

Get rid of ineffective optimizations

The #pragma unroll is not doing anything on any compiler I tried. Forcing unrolling is often not the right thing to do, compilers have pretty good heuristics and will usually do a good job of it themselves.

I would also try removing the _mm_prefetch() statements. The access pattern is very regular, and the prefetcher in contemporary CPUs will pick up that pattern and automatically fetch ahead, even if it is from multiple different lines of the image simultaneously. When that happens, the prefetch instructions are just dead weight. It also doesn't help to prefetch something right before you do the actual load. If anything, you would want to use this pattern:

1. Load the values you need in this iteration.
2. Prefetch the values you need in the next iteration (or as harold suggested, a few iterations ahead if memory access latency is high compared to the work you need to do per iteration).
3. Do the calculations on the values you loaded in step 1.

Don't use leading underscores for your own variables

Some uses of leading underscores are reserved. In particular, anything with two leading underscores is reserved. In general I would recommend not declaring anything with a leading underscore yourself.