How to convert 3 addition and 1 multiply into vectorized SIMD using intrinsic functions C++

user_12892937 0
Viewed 388

I'm working with a problem using 2D prefix sum, also called Summed-Area Table S. For an 2D array I (grayscale image/matrix/etc), its definition is:

S[x][y] = S[x-1][y] + S[x][y-1] - S[x-1][y-1] + I[x][y]
Sqr[x][y] = Sqr[x-1][y] + Sqr[x][y-1] - Sqr[x-1][y-1] + I[x][y]^2

Calculating the sum of a sub-matrix with two corners (top,left) and (bot,right) can be done in O(1):

sum = S[bot][right] - S[bot][left-1] - S[top-1][right] + S[top-1][left-1]

One of my problem is to calculate all possible sub-matrix sum with a constant size (bot-top == right-left == R), which are then used to calculate their mean/variance. And I've vectorized it to the form below.

lineSize is the number of elements to be processed at once. I choose lineSize = 16 because Intel CPU AVX instructions can work on 8 doubles at the same time. It can be 8/16/32/...

#define cell(i, j, w) ((i)*(w) + (j))
const int lineSize = 16; 
const int R = 3; // any integer
const int submatArea = (R+1)*(R+1);
const double submatAreaInv = double(1) / submatArea;
void subMatrixVarMulti(int64* S, int64* Sqr, int top, int left, int bot, int right, int w, int h, int diff, double submatAreaInv, double mean[lineSize], double var[lineSize])
{
  const int indexCache = cell(top, left, w),
        indexTopLeft = cell(top - 1, left - 1, w),
        indexTopRight = cell(top - 1, right, w),
        indexBotLeft = cell(bot, left - 1, w),
        indexBotRight = cell(bot, right, w);
  
  for (int i = 0; i < lineSize; i++) {
    mean[i] = (S[indexBotRight+i] - S[indexBotLeft+i] - S[indexTopRight+i] + S[indexTopLeft+i]) * submatAreaInv;
    var[i] = (Sqr[indexBotRight + i] - Sqr[indexBotLeft + i] - Sqr[indexTopRight + i] + Sqr[indexTopLeft + i]) * submatAreaInv
         - mean[i] * mean[i];
}

How can I optimize the above loop to have the highest possible speed? Readability doesn't matter. I heard it can be done using AVX2 and intrinsic functions, but I don't know how.

Edit: the CPU is i7-7700HQ, kabylake = skylake family

Edit 2: forgot to mention that lineSize, R, ... are already const

avx avx2 intrinsics simd c++
2 Answers

Your compiler can generate AVX/AVX2/AVX-512 instructions for you, but you need to:

  1. Select the latest available architecture when compiling. For example with GCC you might say -march=skylake if you know your code will run on Skylake and later, but does not need to support older CPUs. Without this, AVX instructions cannot be generated.
  2. Add restrict or __restrict to your pointer inputs to tell the compiler they do not overlap. This applies to S and Sqr, as well as mean and var (both pairs have the same type, so the compiler assumes they might overlap, but you know they do not).
  3. Make sure your data is "over-aligned." For example if you want the compiler to use 256-bit AVX2 instructions, you should align your arrays to 256 bits. There are a few ways to do this, such as making a typedef with the alignment, or using alignas() or std::assume_aligned() (available as a GCC attribute prior to C++20). The point is you need the compiler to know that S, Sqr, mean and var are aligned to the largest SIMD vector size available on your target architecture, so that it does not have to generate as much fixup code.
  4. Use constexpr where possible, such as lineSize.

Most importantly, profile to compare performance as you make changes, and look at the generated code (e.g. g++ -S) to see if it looks the way you want it to.

John_Zwinck_43230

I don't think you can perform efficiently this type of sum using SIMD due to the dependencies of the summation.

Instead you can do the computation differently which can be trivially optimized with SIMD:

  1. Compute row-only partial summation. You parallelize it with SIMD by computing simultaneously for multiple rows.
  2. Now with rows summed up, by computing cols-only partial summation to the output using the same SIMD optimization you obtain your desired Summed-Area Table.

You can do the same for both summation and summation of squares.

The only issue is you need extra memory and this type of computation requires more memory accesses. The extra memory is probably a minor thing but more memory access perhaps can be improved by storing the temporary data (the sums of rows) in a cache friendly manner. You'll probably need to experiment with this.

user_108036710
Related