如何使用内部函数 C++ 将 3 个加法和 1 个乘法转换为矢量化 SIMD

Posted 2023-02-16

【中文标题】如何使用内部函数 C++ 将 3 个加法和 1 个乘法转换为矢量化 SIMD

【英文标题】：How to convert 3 addition and 1 multiply into vectorized SIMD using intrinsic functions C++

【发布时间】：2020-07-30 11:32:11

【问题描述】：

我正在使用 2D 前缀总和（也称为 Summed-Area Table S）来解决问题。对于一个二维数组I（灰度图/矩阵/等），其定义为：

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

计算具有两个角(top,left) 和(bot,right) 的子矩阵之和可以在 O(1) 中完成：

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

我的一个问题是计算所有可能的子矩阵和，其大小恒定(bot-top == right-left == R)，然后用于计算它们的均值/方差。我已将其矢量化为以下形式。

lineSize 是一次要处理的元素数。我选择lineSize = 16 是因为 Intel CPU AVX 指令可以同时处理 8 个双精度。可以是 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];

如何优化上述循环以获得尽可能快的速度？可读性并不重要。我听说可以使用 AVX2 和内在函数，但我不知道如何。

编辑：CPU 是 i7-7700HQ，kabylake = skylake 家族

编辑 2：忘了提到 lineSize, R, ... 已经是 const

【问题讨论】：

int64? AVX2 确实有 64 位整数减法，但您必须从 32 位 SIMD 操作中合成 64 位整数乘法。由于您可以做大量的 SIMD 添加工作，这最终可能仍然值得，但请参阅 Fastest way to multiply an array of int64_t?

糟糕，您在乘法之前将 int64_t 转换为 double，所以实际上您需要 How to efficiently perform double/int64 conversions with SSE/AVX?（在 AVX-512 之前没有直接的硬件支持；仿真对于AVX2 的 +-2^51 范围）。 int32_t 会更有效，如果您的数据可以适合它而不会溢出。另外float 而不是double 会节省内存带宽...

不幸的是，平方值意味着它会溢出 4K 及以上的图像。更不用说 float 在任何大于 2^23 的东西上都具有可怕的精度。 @PeterCordes

哦，我明白了，是的，您要减去两个乘积，因此灾难性取消意味着您需要保持很高的精度才能留下任何有效位。不清楚为什么 4k 分辨率意味着S[] 中的数字更大，但可能来自您的用例中的某些内容。但是，在安全的情况下，您可以通过使用更窄的整数来显着加快较小图像的速度，例如有 2 个版本的函数。 （每个 SIMD 向量的元素数量是原来的两倍，与 double 之间的转换要简单得多。）

Sqr[] 包含平方值，因此 int 仅支持 2^15 像素或更小的图像。你能写一个评论来说明如何只做加/减部分吗？我稍后会尝试乘法部分。如果可能的话，你能在一切都是 32 位 int/float 的情况下编写 SIMD 版本吗？

【参考方案1】：

您的编译器可以为您生成 AVX/AVX2/AVX-512 指令，但您需要：

编译时选择最新的可用架构。例如对于 GCC，如果您知道您的代码将在 Skylake 及更高版本上运行，但不需要支持较旧的 CPU，您可能会说 -march=skylake。没有它，就无法生成 AVX 指令。 将restrict 或__restrict 添加到您的指针输入中，以告知编译器它们不重叠。这适用于 S 和 Sqr，以及 mean 和 var（两个对具有相同的类型，因此编译器假定它们可能重叠，但您知道它们不会重叠）。 确保您的数据“过度对齐”。例如，如果您希望编译器使用 256 位 AVX2 指令，您应该将数组对齐到 256 位。有几种方法可以做到这一点，例如使用对齐方式创建 typedef，或使用 alignas() 或 std::assume_aligned()（在 C++20 之前作为 GCC 属性提供）。关键是您需要编译器知道 S、Sqr、mean 和 var 与您的目标架构上可用的最大 SIMD 向量大小对齐，这样它就不必生成尽可能多的修复代码。 尽可能使用constexpr，例如lineSize。

最重要的是，配置文件以在您进行更改时比较性能，并查看生成的代码（例如g++ -S），看看它是否符合您的要求。

【讨论】：

请注意，restrict 不是标准 C++ 中的东西。

在这里拆分循环听起来是个糟糕的主意。 （此答案的旧版本将其作为第 5 点）。您通常需要 更多 计算强度（每次加载数据时的 ALU 工作量），而不是更少。 mean[i+0..3] 已经在向量寄存器中，准备好与做独立减法工作并行相乘！裂变这些循环将花费更多的总指令和更多的缓存流量。 （或者更糟糕的是，如果数组很大并且不适合 L2 缓存，则内存或 L3 流量。）

@PeterCordes：感谢您的反馈，我删除了第 5 点关于拆分循环的内容。

如果您没有 AVX512，您可能需要手动对其进行矢量化处理，因为仅 AVX2 很难打包 int64_t double。最有效的模拟它的方法只适用于有限的范围。 （编译器无法假设，因此即使它确实进行了矢量化，它也可能比通过做出该假设可以做的要慢）How to efficiently perform double/int64 conversions with SSE/AVX?。如果幸运的话，您可能会获得接近全速的自动矢量化。

对于可能未对齐的indexTopLeft 等，对齐数组可能无关紧要。 AVX 不需要对齐向量，而 GCC8 及更高版本，如 clang，通常会使用未对齐的 SIMD 加载自动向量化，而不是尝试到达对齐边界。 （对齐的数据很好；它避免了加载和存储中的缓存行拆分，因此代码生成策略是“乐观的”：在检查对齐的情况下不会浪费任何开销，然后全速运行。但如果未对齐，让硬件相当有效地处理它，所有支持 AVX 的硬件都会这样做。）

【参考方案2】：

由于求和的依赖性，我认为您不能使用 SIMD 有效地执行这种类型的求和。

相反，您可以进行不同的计算，这可以通过 SIMD 进行简单优化：

仅计算行部分求和。您可以通过同时计算多行来将其与 SIMD 并行化。 现在对行求和，通过使用相同的 SIMD 优化计算仅对输出的 cols 部分求和，您可以获得所需的求和面积表。

您可以对求和和平方和执行相同的操作。

唯一的问题是您需要额外的内存，而这种类型的计算需要更多的内存访问。额外的内存可能是一件小事，但可以通过以缓存友好的方式存储临时数据（行的总和）来改善更多的内存访问。您可能需要对此进行试验。

【讨论】：

这实际上是我的第 7 次改进，甚至是在多线程之后。我发现内存访问是压倒性的瓶颈。另外，如何同时将 SIMD 用于多行？每个值都是垂直的，在内存中很远。

@HuyĐứcLê 内存问题取决于平台。在普通笔记本电脑上，单线程能够使内存饱和，但在工作站上，您需要相当多的线程才能做到这一点。

@HuyĐứcLê 在这种算法的情况下，您应该知道缓存行是 64 字节。因此，每当您加载附近 64 字节的内容时，都会在 chache 中高高上传。因此，我相信如果行是连续的，您可以同时加载多个行，然后有效地计算整个缓存。

@HuyĐứcLê 对于临时表，您可能应该以混合行-列数据的格式存储...考虑使用 simd 存储/加载列和行的有效方法（只要它是完成正确对齐）。

@ALX32z 我会想到的。另外，我已经使用了多线程，这肯定会使内存使用饱和。