用于乘以特定数组元素的等效 SIMD 指令答案

【问题标题】：Equivalent SIMD instruction for multiplying specific array elements用于乘以特定数组元素的等效 SIMD 指令
【发布时间】：2015-11-12 00:47:23
【问题描述】：

我刚刚了解了如何获得 2 个数组的点积（如下代码所示）：

int A[8] = {1,2,3,4,5,1,2,3};
int B[8] = {2,3,4,5,6,2,3,4};

float result = 0;

for (int i = 0; i < 8; i ++) {
    result += A[i] * B[i];
}

相当于（在 SIMD 中）：

int A[8] = {1,2,3,4,5,1,2,3};
int B[8] = {2,3,4,5,6,2,3,4};

float result = 0;

__m128 r1 = {0,0,0,0};
__m128 r2 = {0,0,0,0};
__m128 r3 = {0,0,0,0};

for (int i = 0; i < 8; i += 4) {
  float C[4] = {A[i], A[i+1], A[i+2], A[i+3]};
  float D[4] = {B[i], B[i+1], B[i+2], B[i+3]};
  __m128 a = _mm_loadu_ps(C);
  __m128 b = _mm_loadu_ps(D);

  r1 = _mm_mul_ps(a,b);
  r2 = _mm_hadd_ps(r1, r1);
  r3 = _mm_add_ss(_mm_hadd_ps(r2, r2), r3);
  _mm_store_ss(&result, r3);
}

如果我想将数组中不连续的元素相乘，我现在很好奇如何在 SIMD 中获取等效代码。例如，如果我想执行以下操作，SIMD 中的等价物是什么？

int A[8] = {1,2,3,4,5,1,2,3};
int B[8] = {2,3,4,5,6,2,3,4};

float result = 0;
for (int i = 0; i < 8; i++) {
    for (int j = 0; j < 8; j++) {
        result += A[foo(i)] * B[foo(j)]
    }
}

foo 只是一些函数，它返回一个 int 作为输入参数的某个函数。

【问题讨论】：

你自己破解过吗？尝试和失败会比得到答案更快地教会你。
int 和 float 和 __m128 之间的数据类型是否存在混淆？
另外，当您使用它时，请尝试使用 OpenCL。在 Intel 系统上，这让您无需为每个处理器编译代码，同时获得与内部处理器相同的速度。
您应该专注于获得第一个示例，然后再进行更复杂的操作。有几个错误，即使这些错误被修复，代码也远非最佳。由于 SIMD 旨在接近最大的理论最佳案例性能，因此您确实必须更加关注细节，并准备好不断迭代优化过程，直到代码接近完美。
好吧，这可能本身就是一个 Stack Overflow 问题，但您可能应该摆脱 int 到 float 的转换，在整数域中执行所有操作，并将水平操作移到循环，仅用于初学者。您还需要收紧要求，例如最多有多少个元素，输入数据的范围是多少，等等。

标签： c simd code-translation

【解决方案1】：

一般来说，SIMD不喜欢随机访问单个元素之类的事情。但是，仍然有一些技巧可以使用。

如果foo 提供的索引在编译时是已知的，您可能可以打乱两个向量以正确对齐它们的元素。只需查看 swizzle category in the Intrinsics Guide 中的内在函数即可。你肯定需要_mm_shuffle_ps 和_mm_unpackXX_ps 之类的东西。各种移位/对齐指令也可能有用。

使用 AVX2，您可以尝试使用收集指令。对于 float 输入 32 位模式，您可以使用 _mm_i32gather_ps 或 _mm256_i32gather_ps 内在函数。然而，@PaulR 写道 here 说它们并不比普通的标量加载快。

SSSE3 中的_mm_shuffle_epi8 instrinsic 可能会提供另一种解决方案。这是一个很棒的指令，允许以单个字节的粒度执行寄存器内收集操作。但是，创建 shuffle 掩码并不是一项简单的任务。 This paper（阅读第 3.1 和 4 节）展示了如何将此方法扩展到大于一个 XMM 寄存器的输入数组，但对于 64 个或更多元素，它似乎不再比标量代码更好。

【讨论】：

【解决方案2】：

如果我必须做这个任务，我会这样做：

int A[8] = {1,2,3,4,5,1,2,3};
int B[8] = {2,3,4,5,6,2,3,4};

float PA[8], PB[8];
for (int i = 0; i < 8; i++) 
{
    PA[i] = A[foo(i)]; 
    PB[i] = B[foo(i)]; 
}

__m128 sums = _mm_set1_ps(0);
for (int i = 0; i < 8; i++) 
{
    __m128 a = _mm_set1_ps(PA[i]);
    for (int j = 0; j < 8; j += 4) 
    {
        __m128 b = _mm_loadu_ps(PB + j);
        sums = _mm_add_ps(sums, _mm_mul_ps(a, b));
    }
}
float results[4];
_mm_storeu_ps(results, sums);
float result = results[0] + results[1] + results[2] + results[3];

【讨论】：