加速大量与数组相关的计算，Visual Studio

Question

我想知道加速大量数组计算的最佳方法可能是什么。假设我有这种情况：

int template_t[] = {1, 2, 3, 4, 5, 6, ...., 125};
int image[3200][5600];
int template_image[3200][5600];

for(int i = 0; i < 3200; i++) {
    for(int j = 0; j < 5600; j++) {

        // iterate over template to find template value per pixel
        for(int h = 0; h < template_length; h++)
            template_image[i][j] += template_t[h] * image[i][j];

    }
}

当然，我的情况要复杂得多，但同样的想法也适用。我有一些代表图像中像素的大数组，我需要对每个像素应用一些模板数组来计算要放置在模板图像中的值。

我想了几种方法来加快速度：

• SIMD 指令？但是，我似乎无法在 Visual Studio 中找到任何用于编写 SIMD 特定代码的资源。
• 并行化 - 尽管我已经并行化了整个执行本身，所以程序会基于 X 核运行自身的 X 个实例。程序的输入是大量的图像文件，因此这些 X 实例都将处理单独的文件。

什么能让我物超所值？感谢您的任何建议！

Answer 1

首先，对类型使用_t 名称，而不是数组。让我们调用数组template_multipliers[]。

如果template_multipliers 是const，变量是unsigned，编译器可以在编译时对其求和并完全优化掉内循环。

对于 gcc，我们还可以通过将 template_t 的总和提升出循环来获得更好的代码。在这种情况下，即使使用int，而不是unsigned int，它也能在编译时求和。

有符号溢出是未定义的行为，这可能就是为什么 gcc 有时不知道它正在优化的目标机器上实际会发生什么。 （例如，在 x86 上，这不是您需要避免的事情。一系列加法的最终结果不取决于操作的顺序，即使某些命令会在临时结果中产生带符号溢出，而有些则不会。gcc 不会t 总是在有符号的情况下利用加法的关联性）。

这纯粹是 gcc 的限制。您的代码必须避免源级操作顺序中的符号溢出，但如果编译器知道您会通过执行其他更快的操作获得相同的结果，那么它可以而且应该这样做。

---

// aligning the arrays makes gcc's asm output *MUCH* shorter: no fully-unrolled prologue/epilogue for handling unaligned elements
#define DIM1 320
#define DIM2 1000
alignas(32) unsigned int image[DIM1][DIM2];
alignas(32) unsigned int template_image[DIM1][DIM2];

// with const, gcc can sum them at compile time.
const
static unsigned int template_multipliers[] = {1, 2, 3, 4, 5, 6, 7, 8, 8, 10, 11, 12, 13,   125};
const static int template_length = sizeof(template_multipliers) / sizeof(template_multipliers[0]);


void loop_hoisted(void) {
  for(int i = 0; i < DIM1; i++) {
    for(int j = 0; j < DIM2; j++) {
        // iterate over template to find template value per pixel
        unsigned int tmp = 0;
        for(int h = 0; h < template_length; h++)
            tmp += template_multipliers[h];
        template_image[i][j] += tmp * image[i][j];

    }
  }
}

gcc 5.3 与 -O3 -fverbose-asm -march=haswell auto-vectorizes this 与内循环：

# gcc inner loop: ymm1 = set1(215) = sum of template_multipliers
.L2:
    vpmulld ymm0, ymm1, YMMWORD PTR [rcx+rax] # vect__16.10, tmp115, MEM[base: vectp_image.8_4, index: ivtmp.18_90, offset: 0B]
    vpaddd  ymm0, ymm0, YMMWORD PTR [rdx+rax]   # vect__17.12, vect__16.10, MEM[base: vectp_template_image.5_84, index: ivtmp.18_90, offset: 0B]
    vmovdqa YMMWORD PTR [rdx+rax], ymm0       # MEM[base: vectp_template_image.5_84, index: ivtmp.18_90, offset: 0B], vect__17.12
    add     rax, 32   # ivtmp.18,
    cmp     rax, 4000 # ivtmp.18,
    jne     .L2       #,

这是英特尔 Haswell 内部循环中的 9 个融合域微指令，因为 pmulld 在 Haswell 及更高版本上是 2 个微指令（即使使用单寄存器寻址模式也不能微融合）。这意味着循环只能每 3 个时钟运行一次迭代。 gcc 可以通过对目标使用指针增量以及为 src 使用dst + src-dst 2 寄存器寻址模式（因为它不能微熔无论如何）。

请参阅godbolt Compiler Explorer link 获取 OP 代码的较少修改版本的完整源代码，该代码不会提升 template_multipliers 的总和。它使 asm 很奇怪：

    unsigned int tmp = template_image[i][j];
    for(int h = 0; h < template_length; h++)
        tmp += template_multipliers[h] * image[i][j];
    template_image[i][j] = tmp;

.L8:  # ymm4 is a vector of set1(198)
    vmovdqa ymm2, YMMWORD PTR [rcx+rax]       # vect__22.42, MEM[base: vectp_image.41_73, index: ivtmp.56_108, offset: 0B]
    vpaddd  ymm1, ymm2, YMMWORD PTR [rdx+rax]   # vect__1.47, vect__22.42, MEM[base: vectp_template_image.38_94, index: ivtmp.56_108, offset: 0B]
    vpmulld ymm0, ymm2, ymm4  # vect__114.43, vect__22.42, tmp110
    vpslld  ymm3, ymm2, 3       # vect__72.45, vect__22.42,
    vpaddd  ymm0, ymm1, ymm0    # vect__2.48, vect__1.47, vect__114.43
    vpaddd  ymm0, ymm0, ymm3    # vect__29.49, vect__2.48, vect__72.45
    vpaddd  ymm0, ymm0, ymm3    # vect_tmp_115.50, vect__29.49, vect__72.45
    vmovdqa YMMWORD PTR [rdx+rax], ymm0       # MEM[base: vectp_template_image.38_94, index: ivtmp.56_108, offset: 0B], vect_tmp_115.50
    add     rax, 32   # ivtmp.56,
    cmp     rax, 4000 # ivtmp.56,
    jne     .L8       #,

每次循环都会对template_multipliers 进行一些求和。循环中加法的数量取决于数组中的值（而不仅仅是值的数量）。

---

这些优化中的大多数应该适用于 MSVC，除非整个程序链接时优化允许它进行求和，即使 template_multipliers 是非常量。

Answer 2

一个简单的优化，不应该编译器已经为你做了，是：

    int p = template_image[i][j], p2= image[i][j];
    // iterate over template to find template value per pixel
    for(int h = 0; h < template_length; h++)
        p += template_t[h] * p2;

    template[i][j]= p;

进一步看这个，你的模板定义为 1, 2, 3,..125，然后p2 乘以 1*2*3*4..*125，它是常数（让我们调用它CT) 等等：

for (h..
    template_image[i][j] += template_t[h] * image[i][j];

等价于

template_image[i][j] += CT * image[i][j];

所以算法变成：

#define CT 1*2*3*4*5*6*7...*125 // must stil lbe completed
int image[3200][5600];
int template_image[3200][5600];

for(int i = 0; i < 3200; i++) {
    for(int j = 0; j < 5600; j++) {
        template_image[i][j] += CT * image[i][j];
    }
}

这可以通过j 并行化。