是否有机会使用SIMD加速循环代码？

Question

考虑以下代码，其中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）提高这种循环代码的性能？

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

Answer 1

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

此外，未对齐使用可能会导致问题。也可能从末端循环。

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

Answer 2

将表达式编写为矩阵向量积有时会有所帮助。假设您已经知道ₖ₊₈ 您可以计算ₖ 至<代码>sₖ₊₇ 自<代码>aₖ 至<代码>aₖ₊₇ 使用

 [ µ  µ² µ³ µ⁴ µ⁵ µ⁶ µ⁷ µ⁸]   [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ₖ₊₇]   [ µ ] 

第一个矩阵可以分解为三个矩阵，每个矩阵可以使用一个随机和一个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

如果我分析正确，倍数的计算ₖ 可以交错，这将抵消延迟。唯一的热路径是最终的vbroadcastss vfmadd，用于计算ₖₖ₊₇ 自<代码>dₖₖ₊₇ 和<代码>sₖ₊₈ 。计算16个s的块可能是值得的ₖ 和fmaddⁱ * sₖ₊₁₆ 到该值。此外，用于计算<代码>d的FMAₖ 仅使用一半的元素。通过一些复杂的旋转，可以用相同数量的FMA计算两个块（我认为这不值得努力，但可以自由尝试）。

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

N、 B.如果您计算的不是公式，而是：

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

此外，在标量版本中，您将拥有融合多加法，而不是先加法再乘法。结果只会相差µ的一个因子。

Answer 3

相关：是否可以在计算中对串行依赖项使用SIMD，如指数移动平均滤波器？-如果前面有一个n步的封闭式公式，您可以使用它来回避序列依赖项。但我认为情况并非如此。

这看起来像是一种串行依赖的前缀和类型，位于垂直加法的顶部，带有a[j]。有一些方法可以加速，例如获得O（SIMD\u width/log（SIMD\u width））的加速。

• Intel cpu上的SIMD前缀和
• 与SSE的并行前缀（累积）和