有没有机会用 SIMD 加速循环代码？答案

【问题标题】：Any chance to accelerate recurrent code with SIMD?有没有机会用 SIMD 加速循环代码？
【发布时间】：2019-07-06 22:38:05
【问题描述】：

考虑以下代码，其中a 是float 的参数数组，s 是float 的初始未初始化结果数组：

s[n - 1] = mu * a[n - 1];
for (int j = n - 2; j >= 0; j--)
    s[j] = mu * (a[j] + s[j + 1]);
return s;

是否有机会使用 SIMD (AVX2) 提高此类循环代码的性能？

编辑：后来我发现这个公式/算法被称为“折扣和”，但在互联网上找不到它的并行版本。

【问题讨论】：

您确定您的编译器尚未在优化构建中使用 SIMD 指令吗？
@DanielLangr：嗯，什么？ mulss 是标量单一的。我没有看到任何 SIMD 指令，唯一的 ps 指令是寄存器副本（它当然使用 movaps，而不是 movss 与错误依赖项合并）。由于您的i <= n-2 在签名int i 和size_t n-2 之间，ICC 甚至无法优化从int-> 指针宽度的符号扩展。虽然签名溢出 UB 应该允许这样做。

标签： c++ performance vectorization simd avx2

【解决方案1】：

相关：Is it possible to use SIMD on a serial dependency in a calculation, like an exponential moving average filter? - 如果前面有 n 步的封闭式公式，您可以使用它来回避串行依赖项。但我认为情况并非如此。

这看起来像一个前缀和类型的串行依赖，在一个垂直添加的顶部a[j]。有一些方法可以加快速度，例如
O( SIMD_width / log(SIMD_width) )。

【讨论】：

不过，我不知道如何将前缀和算法扩展到折扣 mu。
@SergeRogatch：我还没想好。在某些时候，您可能需要乘以 set(1, mu, mu*mu, mu*mu*mu) 之类的向量。或者在每个洗牌步骤之后乘以set1(mu)，所以当它下降到向量的底部时，你最终得到了j+3元素的mu立方？我可能不会最终回到这个问题，但对于你或其他想要尝试回答这个问题的人来说，这可能是一个很好的起点。

【解决方案2】：

有时将表达式写成矩阵向量积会有所帮助。假设您已经知道 sₖ₊₈，您可以使用

从 aₖ 到 aₖ₊₇ 计算 sₖ 到 sₖ₊₇

 [ µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷ µ⁸]   [aₖ₊₀     ]
 [ 0  µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷]   [aₖ₊₁     ]
 [ 0  0  µ  µ² µ³ µ⁴ µ⁵ µ⁶]   [aₖ₊₂     ]
 [ 0  0  0  µ  µ² µ³ µ⁴ µ⁵]   [aₖ₊₃     ]
 [ 0  0  0  0  µ  µ² µ³ µ⁴] * [aₖ₊₄     ]
 [ 0  0  0  0  0  µ  µ² µ³]   [aₖ₊₅     ]
 [ 0  0  0  0  0  0  µ  µ²]   [aₖ₊₆     ]
 [ 0  0  0  0  0  0  0  µ ]   [aₖ₊₇+sₖ₊₈]

由于sₖ₊₈ 在计算时可能会有一些延迟，因此将其移出产品是有意义的。这可以通过一次广播和一次融合多加来计算：

 [ µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷ µ⁸]   [aₖ₊₀]   [ µ⁸]
 [ 0  µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷]   [aₖ₊₁]   [ µ⁷]
 [ 0  0  µ  µ² µ³ µ⁴ µ⁵ µ⁶]   [aₖ₊₂]   [ µ⁶]
 [ 0  0  0  µ  µ² µ³ µ⁴ µ⁵]   [aₖ₊₃]   [ µ⁵]
 [ 0  0  0  0  µ  µ² µ³ µ⁴] * [aₖ₊₄] + [ µ⁴] * sₖ₊₈
 [ 0  0  0  0  0  µ  µ² µ³]   [aₖ₊₅]   [ µ³]
 [ 0  0  0  0  0  0  µ  µ²]   [aₖ₊₆]   [ µ²]
 [ 0  0  0  0  0  0  0  µ ]   [aₖ₊₇]   [ µ ]

第一个矩阵可以分解为三个矩阵，每个矩阵可以使用一个 shuffle 和一个 FMA 来计算：

 [ 1  0  0  0  µ⁴ 0  0  0 ]   [ 1  0  µ² 0  0  0  0  0 ]   [ µ  µ² 0  0  0  0  0  0 ]   [aₖ₊₀]   [ µ⁸]
 [ 0  1  0  0  µ³ 0  0  0 ]   [ 0  1  µ  0  0  0  0  0 ]   [ 0  µ  0  0  0  0  0  0 ]   [aₖ₊₁]   [ µ⁷]
 [ 0  0  1  0  µ² 0  0  0 ]   [ 0  0  1  0  0  0  0  0 ]   [ 0  0  µ  µ² 0  0  0  0 ]   [aₖ₊₂]   [ µ⁶]
 [ 0  0  0  1  µ  0  0  0 ]   [ 0  0  0  1  0  0  0  0 ]   [ 0  0  0  µ  0  0  0  0 ]   [aₖ₊₃]   [ µ⁵]
 [ 0  0  0  0  1  0  0  0 ] * [ 0  0  0  0  1  0  µ² 0 ] * [ 0  0  0  0  µ  µ² 0  0 ] * [aₖ₊₄] + [ µ⁴] * sₖ₊₈
 [ 0  0  0  0  0  1  0  0 ]   [ 0  0  0  0  0  1  µ  0 ]   [ 0  0  0  0  0  µ  0  0 ]   [aₖ₊₅]   [ µ³]
 [ 0  0  0  0  0  0  1  0 ]   [ 0  0  0  0  0  0  1  0 ]   [ 0  0  0  0  0  0  µ  µ²]   [aₖ₊₆]   [ µ²]
 [ 0  0  0  0  0  0  0  1 ]   [ 0  0  0  0  0  0  0  1 ]   [ 0  0  0  0  0  0  0  µ ]   [aₖ₊₇]   [ µ ]

最右边的矩阵向量乘积实际上是一个乘法。

总体而言，对于 8 个元素，您需要 4 个 FMA、一个乘法和 4 个随机播放/广播（$$$$ 表示任何东西（有限）都可以在这里——或者，如果这些保证为 0，µ 向量可以是部分共享。所有向量都标记为最低有效优先，所有乘法都是元素方式）：

bₖₖ₊₇  = [aₖ₊₀, aₖ₊₁, aₖ₊₂, aₖ₊₃, aₖ₊₄, aₖ₊₅, aₖ₊₆, aₖ₊₇] * [µ  µ  µ  µ  µ  µ  µ  µ ]     vmulps
bₖₖ₊₇ += [aₖ₊₁, $$$$, aₖ₊₃, $$$$, aₖ₊₅, $$$$, aₖ₊₆, $$$$] * [µ² 0  µ² 0  µ² 0  µ² 0 ]     vshufps (or vpsrlq) + vfmadd 

cₖₖ₊₇  =  bₖₖ₊₇ 
cₖₖ₊₇ += [bₖ₊₂, bₖ₊₂, $$$$, $$$$, bₖ₊₆, bₖ₊₆, $$$$, $$$$] * [µ² µ  0  0  µ² µ  0  0 ]     vshufps + vfmadd

dₖₖ₊₇  =  cₖₖ₊₇ 
dₖₖ₊₇ += [cₖ₊₄, cₖ₊₄, cₖ₊₄, cₖ₊₄, $$$$, $$$$, $$$$, $$$$] * [µ⁴ µ³ µ² µ  0  0  0  0 ]     vpermps + vfmadd

sₖₖ₊₇ = dₖₖ₊₇ 
       + [sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈, sₖ₊₈] * [µ⁸ µ⁷ µ⁶ µ⁵ µ⁴ µ³ µ² µ ]     vbroadcastss + vfmadd

如果我分析正确，多个dₖ 的计算可以交错，这会抵消延迟。唯一的热路径是最终的vbroadcastss + vfmadd，从dₖₖ₊₇ 和sₖ₊₈ 计算sₖₖ₊₇。可能值得计算 16 个 sₖ 和 fmadd µⁱ * sₖ₊₁₆ 的块。此外，计算dₖ 的 FMA 仅使用一半的元素。通过一些复杂的调配，可以计算出具有相同数量 FMA 的两个块（我认为这不值得付出努力——但请随意尝试）。

为了比较：纯标量实现需要对 8 个元素进行 8 次加法和 8 次乘法运算，每个操作都取决于之前的结果。

注意如果您计算的不是您的公式，您可以保存一个乘法：

sₖ = aₖ₊₁ + µ*sₖ₊₁

此外，在标量版本中，您将使用 Fused-Multiple-Adds，而不是先加法再乘法。结果只会相差µ。

【讨论】：

【解决方案3】：

如果您的意思是矢量化，那么不，不是直接。由于计算使用上一次迭代的值来计算下一次，因此它不是简单的矢量化。

还没有对齐使用s 可能会导致问题。也可能从最后循环。

如果您只是指 SIMD 集中的指令，那么也许可以使用它们，但不一定会带来巨大的好处，而且编译器通常知道如何做到这一点。

【讨论】：