为什么gcc自动向量化不能处理大于3x3的卷积矩阵？

Question

我已经实现了以下卷积矩阵程序。

#include <stdio.h>
#include <time.h>

#define NUM_LOOP 1000
#define N 128   //input or output dimention 1
#define M N     //input or output dimention 2
#define P 5 //convolution matrix dimention 1 if you want a 3x3 convolution matrix it must be 3
#define Q P     //convolution matrix dimention 2
#define Csize P*Q   
#define Cdiv  1     //div for filter 
#define Coffset 0   //offset 

//functions
void unusual(); //unusual implementation of convolution
void naive();
//data
unsigned short int input[N][M] __attribute__(( aligned(32))); // input data
unsigned short int output[N][M] __attribute__(( aligned(32))); // out put data
unsigned short int kernel[P][Q] __attribute__(( aligned(32)));//convolution coefficients

int main(){
    struct timespec tStart, tEnd;//used to record the processiing time
    double tTotal , tBest=10000;//minimum of toltal time will asign to the best time

    int w=0;
    do{// this loop repeat the body to record the best time
        clock_gettime(CLOCK_MONOTONIC,&tStart);

        //function to be executed here :

        unusual();

        clock_gettime(CLOCK_MONOTONIC,&tEnd);
        tTotal = (tEnd.tv_sec - tStart.tv_sec);
        tTotal += (tEnd.tv_nsec - tStart.tv_nsec) / 1000000000.0;

        if(tTotal<tBest)
            tBest=tTotal;
    } while(w++ < NUM_LOOP);

    printf(" The best time: %lf sec in %d repetition for %dX%d matrix\n",tBest,w, MAX1, MAX2);

    return 0;
}

//unusual sequential convolution
void unusual(){
    int i, j,k,temp;

    for (i=P/2; i< N-P/2; i++){
        for(j=Q/2; j< M-Q/2; j++){
            temp=0;
            for(k=0; k< Csize; k++){
                temp += (kernel[k/P][k%Q]) * (input[i - (P/2) + (k/Q)][j - (Q/2) + (k%Q)]);

            }
            output[i][j]=((temp/(Cdiv))+Coffset);
        }
    }
}
//The naive implementation
inline void naive(){
    int i, j,k,l,temp;
    for (i=P/2; i< N-P/2; i++){
        for(j=Q/2; j< M-Q/2; j++){
            temp=0;

            for(k = 0; k <  P; k++){ 
                for(l = 0; l <  Q; l++){
                    temp += (kernel[k][l]) * (input[i - (P/2)+k][j - (Q/2)+l]);
                }
            }
            output[i][j]=((temp/(Cdiv))+Coffset);
        }
    }
}

问题在于当我使用-O3进行自动向量化时，它仅适用于3x3卷积矩阵。我已经查看了汇编输出，并且自动向量化只对3x3内核进行了一些更改，并合理地提高了性能（20倍的速度，注意：异常函数的标量版本比朴素函数慢），但是对于5x5卷积矩阵没有任何改进。

- Amiri

---

你能把这个完成到可以编译的程度吗？ - harold

将那个矩阵放在单个寄存器中似乎并不是非常有用，因为它是二维的，所以需要进行一些疯狂的洗牌才能将数据与之对齐。但我会再次查看那个巨大的汇编代码混乱，并了解它的作用。 - harold

它对于conv的每个元素使用一个向量，因此向量大小不是问题。寄存器的数量仍然可能是问题，但我真的不确定。当然，在5x5情况下，它不能使用25个向量来保存条目，但可以即时广播一些。尽管如此，GCC仍然不合作，我花了很多时间试图说服它为5x5制作好代码，但它根本不行，即使使用内部函数也不行。嗯，也许如果我变得非常明确并且不留任何东西给优化器，但那样我必须手动展开并失去通用性。 - harold

此外，GCC正在使用dwords进行所有计算，这并不必要（如果Cdiv不再为1，除了某些特定情况和稍微舍入一下，它肯定必须是2的幂，因为没有整数向量除法）。 因此，short temp会产生更好的代码。但对于5x5来说仍然有点太困难了，显然。 GCC还非常小心地不覆盖填充，因此每行都以15个标量卷积开始和结束，使用Intrinsics很容易避免这种情况（顶部和底部仍然应该具有这种情况，否则读取将超出输入）。 - harold

看起来你不是第一个研究这个的人，自动向量化非常复杂，它的工作效果（或者说不工作）非常依赖于你编写代码的方式。重构你的代码结构可能会有很大帮助，我认为在stackoverflow上链接页面并不总是最好的方法，但我想给你链接以下页面：https://locklessinc.com/articles/vectorize/这个人在向量化方面做了相当多的工作，我希望它能帮助你解决问题。 - koldewb

Answer 1

看起来没有人对此问题感兴趣。因此，我将分享我的发现，并在未来更新我的答案。

第一个更新：根据我的经验，gcc -fopt-info-vec 报告对于 Csize <= 16 的向量化。这是因为向量化因子为 16，这也是 gcc 不会向量化其他内核大小的不寻常实现的原因之一。向量化因子是指可以放入矢量中的元素数量。在本例中，short integer 等于 16位 元素。

来自维基百科的信息：

在第一步中，编译器查找可能阻止向量化的障碍。阻碍向量化的一个主要障碍是小于矢量长度的真正数据依赖关系。其他障碍包括函数调用和短迭代计数。

- Amiri

Answer 2

自动矢量化器的主要障碍是非常数循环变量。如果您的实现中使用int Csize = P * Q;，则无法进行矢量化。因此，在协助自动向量化器时，请考虑这一点。这不是问题，因为您像#define Csize那样声明了Csize。但请在您的工作中注意这一点。 然后，您的非常规实现是编译器中的优化方法的朴素实现的循环变换。看起来您破坏了朴素实现。您的发现表明它受到16的限制，因此我展开了您的非常规函数，自动向量化器表示已对其进行了矢量化。

for(k=0; k< P*Q; k+=2){
                temp += (kernel[k/Q][k%Q]) * (input[i - (P/2) + (k/Q)][j - (Q/2) + (k%Q)]);
                temp += (kernel[k/Q][k%Q]) * (input[i - (P/2) + ((k+1)/Q)][j - (Q/2) + ((k+1)%Q)]);
}

它也适用于7x7内核：

for(k=0; k< P*Q; k+=4){//IACA_START
                temp += (kernel[k/Q][k%Q]) * (input[i - (P/2) + (k/Q)][j - (Q/2) + (k%Q)]);
                temp += (kernel[k/Q][k%Q]) * (input[i - (P/2) + ((k+1)/Q)][j - (Q/2) + ((k+1)%Q)]);
                temp += (kernel[k/Q][k%Q]) * (input[i - (P/2) + ((k+2)/Q)][j - (Q/2) + ((k+2)%Q)]);
                temp += (kernel[k/Q][k%Q]) * (input[i - (P/2) + ((k+3)/Q)][j - (Q/2) + ((k+3)%Q)]);
}

你不需要自己展开它，你可以通过使用#pragma属性来强制编译器展开或更改循环结构。这是因为编译器使用SLP概念进行自动向量化，有趣的是SLP基于展开！

- ADMS

Answer 3

我猜测它无法进行优化，可能是由于内存对齐问题导致的。您已经指定卷积为2字节短整型。大多数SSE函数喜欢使用128位向量，而AVX则喜欢512位向量。

在我的机器上，我像这样声明了conv：

uint16_t conv[Cdim1][8] = {0}; //You need to pad extra fields with zeroes

然后将内部循环替换为以下内容：

for(ki = 0; ki < Cdim; ++ki) 
    for(kj = 0; kj < 8; ++kj)
        temp += (conv[ki][kj]) * (input[i - (Cdim1/2) + ki][j - (Cdim2/2) + kj]);

使用以下命令编译：gcc so.c -Wall -Wextra -Ofast -mtune=native，可以让我获得向量优化！

注意事项：

• 不要使用 8，尝试找到最小所需填充并制作宏，使其适用于维度≥8的卷积矩阵。
• 用一些零来填充输入，以消除末尾的未定义行为。
• 请注意，这实际上并没有提高性能。事实上它会变慢！
• 请注意，如果您按照以下顺序执行循环 for(ki) for(i) for(j) for(kj)，则可以挤出几个周期。这可能是因为每个卷积行可以存储更长时间，从而减轻寄存器压力。这也可能是我的 CPU 上的故障。
• 在声明变量时考虑使用__attribute__ ((aligned (8)))。在这种情况下，它没改变任何事情，但在优化时也要考虑这个。当然，这只适用于 GCC，您需要其他技巧来应对 MSVC。

- vguberinic

---

谢谢，但你改变了内部循环。正如我在评论中提到的，我会更新问题。gcc可以向量化四个循环实现。我的gcc 6.2将3x3和5x5内核向量化，速度提高了约6倍。即使是7x7和9x9，gcc也有所改进。顺便说一下，问题是关于不寻常的实现方式，称为“seqconv”，而不是朴素的实现。 - Amiri