Parallel algorithm to make three-layered independent operations with nontrivial indexing vectorizable using icc

user_13241857 0
Viewed 70

The current example shows an execution of N=10 independent operations inside a three-layered loop, but unfortunately the intel compiler autovectorization computes the cost at loop level, and when successfull at the innermost level it refuses to consider the vectorization of the two 'external' loops.

A solution would be to express the same with only one loop, where 'if-else' conditions are allowed but perhaps not encouraged.

The indexing covers a subset of the total possible combinations, if someone knows the formal name of it it would be nice.

A working minimal example for L=5 (i.e. 10 possible combinations) follows. The production-code makes heavy use of the case with L=6, which has 20 operations that are performed ~100k times (i.e. for different arrays of length 6).

Note: auto-vectorization of the innermost loop has been tested without the prints using std::cout at each iteration :-)

Note: for curiosity, the (current) strategy of precomputing the elements (with only the inner loop vectorized, intel report showing "estimated potential speedup: 1.8") as opposed to computing them on the fly only when needed is making the section 12% slower.

#include <iostream>

int vectorizable_function(int &A, int &B, int &C) { return A + B - C; }


void populate_combinations(int *container, int L, int *data){

    // Intel vectorization report provides feedback
    // 
    //

    int counter=0;
    for (int i=0; i<L-2; i++){
        for (int j=i+1; j<L-1; j++){
            for (int k=j+1; k<L; k++){
                // std::cout << "i,j,k are " << i << " " << j << " " << k << '\n';
                container[counter] = vectorizable_function( data[i], data[j], data[k] );
                counter++;
            }
        }
    }
    std::cout << "counter was " << counter << std::endl;
}


// Function we want to populate with the parallelizable algorithm
//
void vectorizable_populate_combinationsL5(int * container, int *data){
    #pragma vector always
    for (int i=0; i<10; i++){
        // Smart algorithm goes here?
    }
}

int main(){

    const int L = 5;

    const int TOTALNCOMBINATIONS = 10;
    int combinations[10];

    int data[5] = {5,6,8,12,3};

    populate_combinations(combinations, L, data);

    std::cout << "Results are:\n";
    for (int i=0; i<TOTALNCOMBINATIONS; i++) std::cout << combinations[i] << std::endl;

}

Responding to a comment:

Is your real function also separable the way this one is, where you can compute data[i] + data[j] and reuse it for all data[k]? Also commutative between i and j since they're both getting added, only the k element is getting subtracted. Although that probably doesn't matter, since j runs from [i+1 .. L).

The real function is shown in the following fragment:

double dxAll[18]; // declared and populated long before it's needed

int counter = 0;
    for (int i=0; i<4; i++){
      for (int j=i+1; j<5; j++){
        for (int k=j+1; k<6; k++){
              denomAll[counter] = dxAll[i]*dxAll[6 + j]*dxAll[12 + k] +
                dxAll[6 + i]*dxAll[12 + j]*dxAll[k] +
                dxAll[12 + i]*dxAll[j]*dxAll[6 + k] -
                dxAll[i]*dxAll[12 + j]*dxAll[6 + k] -
                dxAll[6 + i]*dxAll[j]*dxAll[12 + k] -
                dxAll[12 + i]*dxAll[6 + j]*dxAll[k];
      counter++;
        }
      }
    }
auto-vectorization avx512 icc intel c++
1 Answers

Write all your combinations on paper:

0 1 2
0 1 3
0 1 4
0 2 3
0 2 4
0 3 4
1 2 3
1 2 4
1 3 4
2 3 4

Note that there are 10 rows with 3 numbers, 30 numbers total. You need 30 bytes to encode that.

void vectorizable_populate_combinationsL5(int *__restrict container, int *__restrict ){
    static const uint8_t index_i[] = {0, 0, 0, 0, 0, 0, 1, 1, 1, 2};
    static const uint8_t index_j[] = {1, 1, 1, 2, 2, 3, 2, 2, 3, 3};
    static const uint8_t index_k[] = {2, 3, 4, 3, 4, 4, 3, 4, 4, 4};

    #pragma vector always
    for (int counter = 0; counter < 10; counter++)
    {
        int i = index_i[counter];
        int j = index_j[counter];
        int k = index_k[counter];
        container[i] = vectorizable_function(data[i], data[j], data[k]);
    }
}

Using __restrict pointers will let the compiler know that there's no aliasing between input and output, letting it reuse temporaries in registers and otherwise do CSE (Common subexpression elimination).

This does fully unroll to scalar add instructions with clang/LLVM (including the ICX OneAPI compiler which is based on LLVM).

But with the classic ICC, it not only unrolls, it also actually vectorizes. (https://godbolt.org/z/1f7rj1jKc). But unfortunately, to twice as many instructions than LLVM's scalar version, including 3x gather and 1x scatter (with -march=skylake-avx512)! However, it might not need many more instructions for a more complex vectorizable_function with the same input data pattern.

(Without __restrict, ICC doesn't vectorize, and does 2 loads + a store for every output element, not reusing loads.)

user_5098680
Related