从二阶导数计算出的曲线的 SIMD 优化
SIMD optimization of a curve computed from the second derivative
这个问题真是好奇
我正在将例程转换为 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-like 实现同样有用，因为我试图了解如何提升数据依赖性，而我为了优化而尝试优化.
最后，如果你能推荐一些文学作品，请推荐。不确定如何 Google 这个主题。
我唯一能想到的就是水平添加。假设您有 __m128i 个内容为 {pc, 0, pi, pis} 的向量。然后首先HADD将它变成
{pc, pi + pis}, 第二个 HADD 将进入 pc + pi + pis.
HADD 同时对两个 __m128i 进行操作，因此可以加快速度。
但是交错指令使得流水线总是满的并不是一件容易的事。 Link 到 HADD：https://msdn.microsoft.com/en-us/library/bb531452(v=vs.120).aspx
让我将 link 添加到非常有用的关于花车 HADD 的讨论中。很多代码和结论可以直接应用于整数 HADD：Fastest way to do horizontal float vector sum on x86
所以您只是在寻找一种方法来生成 4 个 phase_current 值的向量，您可以将其作为 arg 传递给任意函数。
TL:DR：设置增量和步长的初始向量，以便每个向量元素在序列中跨步 4，为您提供 phase_current[i+0..i+3] 的向量只有两个向量 ADD 操作（垂直的，不是水平的）。 您可以使用代数/数学分解出这种序列相关性。
这有点像前缀和（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_step.
实现这个的一个选择是 运行 标量算法在 SIMD 向量的 4 个独立元素中独立，每次偏移 1 次迭代。对于整数加法，尤其是在向量整数加法延迟只有 1 个周期的 Intel CPU 上，运行ning 4 次迭代 运行ning-total 很便宜，我们可以在调用 USEFUL_FUNC。这将是一种生成矢量输入的方法 USEFUL_FUNC 做与标量代码一样多的工作（假设 SIMD 整数加法与标量整数加法一样便宜，如果我们受到数据依赖性的限制，这通常是正确的每个时钟添加 2 个）。
上述方法更通用一些，可用于此问题的变体，其中存在真正的串行依赖性，我们无法通过简单的数学运算以低廉的成本消除。
如果我们够聪明，我们甚至可以做得比前缀和或蛮力 运行一次处理 4 个序列更好。理想情况下 我们可以推导出一种封闭形式的方法，在值序列 中步进 4（或者无论 SIMD 向量宽度是多少，乘以 [=21 的多个累加器所需的任何展开因子=]).
对 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. This sequence goes 0, 1, 3, 6, 10, 15, 21, 28, ... (https://oeis.org/A000217)。 （我已经排除了最初的 phase_increment）。
技巧是按照这个顺序进行 4 次。 (n+4)*(n+5)/2 - n*(n+1)/2 simplifies 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://whosebug.com/tags/sse/info, especially slides from SIMD at Insomniac Games (GDC 2015)，其中有一些关于如何一般地考虑 SIMD 以及如何布置您的数据，以便您可以使用它。
这个问题真是好奇
我正在将例程转换为 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-like 实现同样有用，因为我试图了解如何提升数据依赖性，而我为了优化而尝试优化.
最后，如果你能推荐一些文学作品，请推荐。不确定如何 Google 这个主题。
我唯一能想到的就是水平添加。假设您有 __m128i 个内容为 {pc, 0, pi, pis} 的向量。然后首先HADD将它变成
{pc, pi + pis}, 第二个 HADD 将进入 pc + pi + pis.
HADD 同时对两个 __m128i 进行操作，因此可以加快速度。
但是交错指令使得流水线总是满的并不是一件容易的事。 Link 到 HADD：https://msdn.microsoft.com/en-us/library/bb531452(v=vs.120).aspx
让我将 link 添加到非常有用的关于花车 HADD 的讨论中。很多代码和结论可以直接应用于整数 HADD：Fastest way to do horizontal float vector sum on x86
所以您只是在寻找一种方法来生成 4 个 phase_current 值的向量，您可以将其作为 arg 传递给任意函数。
TL:DR：设置增量和步长的初始向量，以便每个向量元素在序列中跨步 4，为您提供 phase_current[i+0..i+3] 的向量只有两个向量 ADD 操作（垂直的，不是水平的）。 您可以使用代数/数学分解出这种序列相关性。
这有点像前缀和（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_step.
实现这个的一个选择是 运行 标量算法在 SIMD 向量的 4 个独立元素中独立，每次偏移 1 次迭代。对于整数加法，尤其是在向量整数加法延迟只有 1 个周期的 Intel CPU 上，运行ning 4 次迭代 运行ning-total 很便宜，我们可以在调用 USEFUL_FUNC。这将是一种生成矢量输入的方法 USEFUL_FUNC 做与标量代码一样多的工作（假设 SIMD 整数加法与标量整数加法一样便宜，如果我们受到数据依赖性的限制，这通常是正确的每个时钟添加 2 个）。
上述方法更通用一些，可用于此问题的变体，其中存在真正的串行依赖性，我们无法通过简单的数学运算以低廉的成本消除。
如果我们够聪明，我们甚至可以做得比前缀和或蛮力 运行一次处理 4 个序列更好。理想情况下 我们可以推导出一种封闭形式的方法，在值序列 中步进 4（或者无论 SIMD 向量宽度是多少，乘以 [=21 的多个累加器所需的任何展开因子=]).
对 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. This sequence goes 0, 1, 3, 6, 10, 15, 21, 28, ... (https://oeis.org/A000217)。 （我已经排除了最初的 phase_increment）。
技巧是按照这个顺序进行 4 次。 (n+4)*(n+5)/2 - n*(n+1)/2 simplifies 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://whosebug.com/tags/sse/info, especially slides from SIMD at Insomniac Games (GDC 2015)，其中有一些关于如何一般地考虑 SIMD 以及如何布置您的数据，以便您可以使用它。