从二阶导数计算的曲线的 SIMD 优化答案

【问题标题】：SIMD optimization of a curve computed from the second derivative从二阶导数计算的曲线的 SIMD 优化
【发布时间】：2018-06-07 14:17:54
【问题描述】：

这个问题实在是太好奇了。

我正在将例程转换为 SIMD 指令（而且我对 SIMD 编程很陌生），并且在使用以下代码时遇到了问题：

// args:
uint32_t phase_current;
uint32_t phase_increment;
uint32_t phase_increment_step;

for (int i = 0; i < blockSize; ++i)
{
    USEFUL_FUNC(phase_current);
    phase_increment += phase_increment_step;
    phase_current += phase_increment;
}

问题：假设USEFUL_FUNC 有一个SIMD 实现，我只是想计算一个正确的phase_current 向量进行处理，那么处理依赖的phase_current 的正确方法是什么它以前的值？

反过来，类似fold 的函数式编程实现也同样有用，因为我试图了解如何更多地提升数据依赖项，而我正在尝试优化以进行优化。

最后，如果你能推荐一些文学作品，请做。不知道如何通过 Google 搜索该主题。

【问题讨论】：

我不明白这个问题，
更简洁，如何处理 pi += pis;个人计算机 += 圆周率；使用 SIMD 向量指令？（具体来说，是否可以矢量化 phase_current 及其增量，因为 phase_increment 取决于其先前的状态，所有值都是正确的。）
我不确定这对于一般的 USEFUL_FUNC 是否可行，因为 phase_current 取决于其先前的值。但是，如果您实际上是在询问计算函数的第二个反导数，正如您的标题所暗示的那样，您的代码可能无法反映您的意思。也许您希望 USEFUL_FUNC 计算一些东西并将结果存储在一个数组中？
@TobiasRibizel 我需要弄清楚如何重写问题。这里的问题是矢量化phase_current，最好使用 SIMD 指令。我在问是否有办法避免phase_current 在计算同一事物时对其先前值的依赖。（很明显，phase_increment 斜率是phase_increment_step，所以我们可以向量化phase_increment 值。）
那么我的回答是，对于一般的USEFUL_FUNC，这是不可能的。如果将 pc0 称为 pc 的初始值，pc1 取决于 pc0 和 pi，pc2 取决于 pc1 和 pi，因此您需要在计算 pc2 之前完成 pc1 的计算，这使得 pc1 和 pc2 的任何并行计算（尤其是通过 SIMD ) 不可能。

标签： c++ optimization vectorization simd

【解决方案1】：

我唯一能想到的是水平添加。想象一下，您有 __m128i 向量，其内容为 {pc, 0, pi, pis}。然后首先 HADD 将它变成 {pc, pi + pis} 和第二个 HADD 将使其变为 pc + pi + pis。

HADD 同时在两个 __m128i 上运行，因此可以进行一些加速。

但是交错指令使得管道总是满的并不是一件容易的事。 HADD 链接：https://msdn.microsoft.com/en-us/library/bb531452(v=vs.120).aspx

让我添加指向非常有用的讨论 wrt HADD for floats 的链接。很多代码和结论可以直接应用于整数HADD：Fastest way to do horizontal float vector sum on x86

【讨论】：

谢谢，非常感谢，因为我没有考虑过。两个增量的问题对于向量化/并行化来说并不是一件小事。
用于浮动讨论的 HADD 非常有用。谢谢！
我对水平总和问题的回答最大的一点是，如果需要，您应该不使用hadd / phadd 指令 (_mm_hadd_epi32)快速运行（就像这里它在一个内部循环中，@ilzxc）。如果您需要移入零，您应该使用像 pshufd 或 psrldq (_mm_srli_si128) 这样的常规随机播放。我建议 not 链接到 MSVC 的 hadd 内在函数。这个问题看起来更像是一个前缀和。请参阅 stackoverflow.com/questions/10587598/… 了解每个 4x 32 位元素向量仅使用 2 次随机播放的版本。
哦，等等，它不像前缀和那么难，因为increment_step 是恒定的。使用恒定的二阶导数，数学允许仅使用垂直加法就可以跨过序列而无需额外工作.见my answer

【解决方案2】：

所以您只是在寻找一种方法来生成 4 个 phase_current 值的向量，您可以将其作为 arg 传递给任意函数。

TL:DR：为增量和步进设置初始向量，使每个向量元素在序列中跨步 4，为您提供 phase_current[i+0..i+3] 的向量，仍然只有两个向量加法操作（垂直，而不是水平的）。 这种串行依赖关系可以通过代数/数学来分解。

这有点像前缀和 (which you can SIMD with log2(vector_width) shuffle+add operations for vectors with vector_width elements。) 您还可以使用两步计算将前缀和与多个线程并行化，其中每个线程对数组的一个区域进行前缀求和，然后组合结果和让每个线程将其目标数组的区域偏移一个常数（该区域第一个元素的总和。也请参阅多线程的链接问题。

但是你有一个巨大的简化，phase_increment_step（你想要的值的二阶导数）是常数。我假设 USEFUL_FUNC(phase_current); 按值而不是通过非常量引用获取其 arg，因此对 phase_current 的唯一修改是循环中的 +=。而 useful_func 不能以某种方式改变 increment 或 increment_step。

实现这一点的一个选项是在 SIMD 向量的 4 个独立元素中独立运行标量算法，每次偏移 1 次迭代。对于整数加法，尤其是在向量整数加法延迟仅为 1 个周期的 Intel CPU 上，运行总运行次数的 4 次迭代很便宜，我们可以在调用 USEFUL_FUNC 之间执行此操作。这将是一种为USEFUL_FUNC生成向量输入的方法，其工作量与标量代码一样多（假设SIMD整数加法与标量整数加法一样便宜，如果我们受到数据依赖性限制为每个时钟2次加法，这基本上是正确的)。

上述方法更通用一些，对于这个问题的变体可能很有用，因为我们无法通过简单的数学廉价地消除真正的串行依赖。

如果我们够聪明，我们甚至可以做得比前缀求和或蛮力一次运行 4 个序列更好。理想情况下，我们可以推导出一种封闭形式的方法，在值序列中逐步增加 4（或者无论 SIMD 向量宽度是多少，乘以您想要的多个累加器的展开因子，用于 USEFUL_FUNC）。

对step、step*2、step*3、...的序列求和将给我们一个恒定的时间Gauss's closed-form formula for the sum of integers up to n:sum(1..n) = n*(n+1)/2。此序列为 0、1、3、6、10、15、21、28，... (https://oeis.org/A000217)。（我已经排除了最初的phase_increment）。

诀窍是在这个序列中增加 4。 (n+4)*(n+5)/2 - n*(n+1)/2simplifies down to 4*n + 10。再次对其求导，我们得到 4。但是要在第二个积分中走 4 步，我们有4*4 = 16。所以我们可以维护一个向量 phase_increment，我们用 SIMD 添加一个向量 16*phase_increment_step。

我不完全确定我的计步推理是否正确（4 到 16 的额外因素有点令人惊讶）。 计算出正确的公式并在向量序列中取一阶和二阶差分可以非常清楚地知道这是如何计算的：

 // design notes, working through the first couple vectors
 // to prove this works correctly.

S = increment_step (constant)
inc0 = increment initial value
p0 = phase_current initial value

// first 8 step-increases:
[ 0*S,  1*S,   2*S,  3*S ]
[ 4*S,  5*S,   6*S,  7*S ]

// first vector of 4 values:
[ p0,  p0+(inc0+S),  p0+(inc0+S)+(inc0+2*S),  p0+(inc0+S)+(inc0+2*S)+(inc0+3*S) ]
[ p0,  p0+inc0+S,  p0+2*inc0+3*S,  p0+3*inc0+6*S ]  // simplified

// next 4 values:
[ p0+4*inc0+10*S,  p0+5*inc0+15*S,  p0+6*inc0+21*S,  p0+7*inc0+28*S ]

使用这个和之前的4*n + 10 公式：

// first 4 vectors of of phase_current
[ p0,              p0+1*inc0+ 1*S,  p0+2*inc0+3*S,   p0+ 3*inc0+ 6*S ]
[ p0+4*inc0+10*S,  p0+5*inc0+15*S,  p0+6*inc0+21*S,  p0+ 7*inc0+28*S ]
[ p0+8*inc0+36*S,  p0+9*inc0+45*S,  p0+10*inc0+55*S, p0+11*inc0+66*S ]
[ p0+12*inc0+78*S,  p0+13*inc0+91*S,  p0+14*inc0+105*S, p0+15*inc0+120*S ]

 first 3 vectors of phase_increment (subtract consecutive phase_current vectors):
[ 4*inc0+10*S,     4*inc0 + 14*S,   4*inc0 + 18*S,   4*inc0 + 22*S  ]
[ 4*inc0+26*S,     4*inc0 + 30*S,   4*inc0 + 34*S,   4*inc0 + 38*S  ]
[ 4*inc0+42*S,     4*inc0 + 46*S,   4*inc0 + 50*S,   4*inc0 + 54*S  ]

 first 2 vectors of phase_increment_step:
[        16*S,              16*S,            16*S,            16*S  ]
[        16*S,              16*S,            16*S,            16*S  ]
Yes, as expected, a constant vector works for phase_increment_step

所以我们可以使用英特尔的 SSE/AVX 内部函数编写这样的代码：

#include <stdint.h>
#include <immintrin.h>

void USEFUL_FUNC(__m128i);

// TODO: more efficient generation of initial vector values
void double_integral(uint32_t phase_start, uint32_t phase_increment_start, uint32_t phase_increment_step, unsigned blockSize)
{

    __m128i pstep1 = _mm_set1_epi32(phase_increment_step);

    // each vector element steps by 4
    uint32_t inc0=phase_increment_start, S=phase_increment_step;
    __m128i pincr  = _mm_setr_epi32(4*inc0 + 10*S,  4*inc0 + 14*S,   4*inc0 + 18*S,   4*inc0 + 22*S);

    __m128i phase = _mm_setr_epi32(phase_start,  phase_start+1*inc0+ 1*S,  phase_start+2*inc0+3*S,   phase_start + 3*inc0+ 6*S );
     //_mm_set1_epi32(phase_start); and add.
     // shuffle to do a prefix-sum initializer for the first vector?  Or SSE4.1 pmullo by a vector constant?

    __m128i pstep_stride = _mm_slli_epi32(pstep1, 4);  // stride by pstep * 16
    for (unsigned i = 0; i < blockSize; ++i)  {
        USEFUL_FUNC(phase);
        pincr = _mm_add_epi32(pincr, pstep_stride);
        phase = _mm_add_epi32(phase, pincr);
    }
}

进一步阅读：更多关于 SIMD 的一般信息，但主要是 x86 SSE/AVX，请参阅 https://stackoverflow.com/tags/sse/info，尤其是来自 SIMD at Insomniac Games (GDC 2015) 的幻灯片，其中有一些关于如何总体上考虑 SIMD 以及如何布局的好东西您的数据，以便您可以使用它。

【讨论】：