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

Question

这个问题真的很奇怪。

我正在将一个例程转换为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依赖于它以前的值的正确方法是什么？

反过来，类似于函数编程的折叠实现也同样有用，因为我试图理解如何更好地提升数据依赖性，而不是为了优化而进行优化。

最后，如果你能推荐一些文献，请做。不知道如何谷歌这个话题。

Answer 1

我唯一能想到的就是水平加法。假设您有内容为{pc，0，pi，pis}的\uu m128i向量。然后，第一个HADD将把它变成{pc，pi-pis}，第二个HADD将把它变成pc-pi-pis。

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

但是交错指令使得管道总是满的不会是一个微不足道的练习。链接到HADD：https://msdn.microsoft.com/en-us/library/bb531452（v=vs.120）. aspx

让我添加一个链接到非常有用的关于浮点数的讨论wrt HADD。很多代码和结论可以直接应用于整数HADD：在x86上进行水平浮点向量和的最快方法

Answer 2

因此，您只是在寻找一种生成4 phase\u电流值向量的方法，可以将其作为参数传递给任意函数。

TL:DR：设置增量和步长的初始向量，使每个向量元素在序列中移动4步，得到相位电流[i 0..i 3]的向量，仍然只有两个向量相加操作（垂直，而不是水平）。这种串行依赖关系可以通过代数/数学来解决。

这有点像前缀和（您可以使用log2（vector_width）SIMD对具有vector_width元素的向量进行随机添加操作。）您还可以使用两步计算并行化多个线程的前缀和，其中每个线程前缀和数组的一个区域，然后组合结果并让每个线程将其目标数组的区域偏移一个常量（该区域第一个元素的总和。有关多线程，请参阅链接问题。

但是你有一个巨大的简化，即相位增量步长（你想要的值的二阶导数）是常数。我假设有用的FUNC（phase\u current）通过值而不是非常量引用获取其参数，因此对phase\u current的唯一修改是通过循环中的=。而且，有用的函数不能以某种方式改变增量或增量步长。

实现这一点的一个选项是，只需在SIMD向量的4个独立元素中独立运行标量算法，每次偏移1次迭代。对于整数加法，尤其是在Intel CPU上，矢量整数加法延迟仅为1个周期，运行总计的4次迭代是很便宜的，我们可以在调用有用的FUNC之间这样做。这将是一种为有用的\u FUNC生成向量输入的方法，它所做的工作与标量代码一样多（假设SIMD integer add与标量integer add一样便宜，如果我们受数据依赖性的限制，每个时钟只能进行2次加法，这一点通常是正确的）。

上面的方法有点通用，对于存在真正的串行依赖关系的问题的变体可能很有用，我们无法用简单的数学廉价地消除这种依赖关系。

如果我们聪明的话，我们甚至可以做得比前缀和或暴力一步一步运行4个序列更好。理想情况下，我们可以推导出一种封闭形式的方法，按值的顺序逐步执行4步（或者，无论SIMD向量宽度是多少，乘以多个累加器所需的展开因子，以获得有用的函数）。

对一系列步骤求和，<代码>步骤*2，<代码>步骤*3。。。将给我们一个常数乘以高斯闭式公式，用于整数之和小于等于n：和（1..n）=n*（n 1）/2。这个序列是0，1，3，6，10，15，21，28。。。(https://oeis.org/A000217)（我已经计算出了初始阶段增量）。

在这个序列中，诀窍是4<代码>（n 4）*（n 5）/2-n*（n 1）/2简化为4*n 10。再对其求导，我们得到4。但要在第二个积分中进行4步，我们需要4*4=16。因此，我们可以保持一个向量，我们使用SIMD add将其增量，向量为16*phase\u increment\u 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

因此，我们可以使用Intel的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，尤其是《失眠症游戏》（GDC 2015）中SIMD的幻灯片，其中有一些关于如何看待SIMD的好东西，以及如何布局数据以便使用它。