在 C++ 中使用 SSE2 SIMD 对两个数组求和的正确方法
The correct way to sum two arrays with SSE2 SIMD in C++
让我们从以下内容开始:
#include <vector>
#include <random>
using namespace std;
现在,假设一个人有以下三个std:vector<float>
:
N = 1048576;
vector<float> a(N);
vector<float> b(N);
vector<float> c(N);
default_random_engine randomGenerator(time(0));
uniform_real_distribution<float> diceroll(0.0f, 1.0f);
for(int i-0; i<N; i++)
{
a[i] = diceroll(randomGenerator);
b[i] = diceroll(randomGenerator);
}
现在,假设需要按元素对 a
和 b
求和并将结果存储在 c
中,其标量形式如下所示:
for(int i=0; i<N; i++)
{
c[i] = a[i] + b[i];
}
上面代码的 SSE2 向量化版本是什么,请记住输入是上面定义的 a
和 b
(即作为 float
的集合) ehe 输出是 c
(也是 float
的集合)?
经过大量研究,我得出以下结论:
for(int i=0; i<N; i+=4)
{
float a_toload[4] = { a[i], a[i + 1], a[i + 2], a[i + 3] };
float b_toload[4] = { b[i], b[i + 1], b[i + 2], b[i + 3] };
__m128 loaded_a = _mm_loadu_ps(a_toload);
__m128 loaded_b = _mm_loadu_ps(b_toload);
float result[4] = { 0, 0, 0, 0 };
_mm_storeu_ps(result, _mm_add_ps(loaded_a , loaded_b));
c[i] = result[0];
c[i + 1] = result[1];
c[i + 2] = result[2];
c[i + 3] = result[3];
}
然而,这似乎真的很麻烦,而且效率肯定很低:上面的 SIMD 版本实际上比初始标量版本慢三倍(当然,在 Microsoft VS15 的发布模式下进行了优化,并经过 100 万次迭代,而不仅仅是 12 次)。
您不需要中间数组来加载到 SSE 寄存器。只需直接从您的阵列加载。
auto loaded_a = _mm_loadu_ps(&a[i]);
auto loaded_b = _mm_loadu_ps(&b[i]);
_mm_storeu_ps(&c[i], _mm_add_ps(loaded_a, loaded_b));
您也可以省略两个 loaded
变量并将它们合并到 add 中,尽管编译器应该为您完成。
你需要注意这一点,因为如果向量大小不是 4 的倍数,它将无法正常工作(你将访问数组的末尾,导致未定义的行为,并且写入超过 c
的结尾可能会造成破坏)。
您的 for 循环可以简化为
const int aligendN = N - N % 4;
for (int i = 0; i < alignedN; i+=4) {
_mm_storeu_ps(&c[i],
_mm_add_ps(_mm_loadu_ps(&a[i]),
_mm_loadu_ps(&b[i])));
}
for (int i = alignedN; i < N; ++i) {
c[i] = a[i] + b[i];
}
补充说明:
- 处理最后几个浮点数的小循环很常见,当
N%4 != 0
或 N 在编译时未知时,这是强制性的。
- 我注意到您选择了未对齐的版本 load/store,与对齐版本相比,惩罚很小。我在 Whosebug 上找到了这个 link:Is the SSE unaligned load intrinsic any slower than the aligned load intrinsic on x64_64 Intel CPUs?
让我们从以下内容开始:
#include <vector>
#include <random>
using namespace std;
现在,假设一个人有以下三个std:vector<float>
:
N = 1048576;
vector<float> a(N);
vector<float> b(N);
vector<float> c(N);
default_random_engine randomGenerator(time(0));
uniform_real_distribution<float> diceroll(0.0f, 1.0f);
for(int i-0; i<N; i++)
{
a[i] = diceroll(randomGenerator);
b[i] = diceroll(randomGenerator);
}
现在,假设需要按元素对 a
和 b
求和并将结果存储在 c
中,其标量形式如下所示:
for(int i=0; i<N; i++)
{
c[i] = a[i] + b[i];
}
上面代码的 SSE2 向量化版本是什么,请记住输入是上面定义的 a
和 b
(即作为 float
的集合) ehe 输出是 c
(也是 float
的集合)?
经过大量研究,我得出以下结论:
for(int i=0; i<N; i+=4)
{
float a_toload[4] = { a[i], a[i + 1], a[i + 2], a[i + 3] };
float b_toload[4] = { b[i], b[i + 1], b[i + 2], b[i + 3] };
__m128 loaded_a = _mm_loadu_ps(a_toload);
__m128 loaded_b = _mm_loadu_ps(b_toload);
float result[4] = { 0, 0, 0, 0 };
_mm_storeu_ps(result, _mm_add_ps(loaded_a , loaded_b));
c[i] = result[0];
c[i + 1] = result[1];
c[i + 2] = result[2];
c[i + 3] = result[3];
}
然而,这似乎真的很麻烦,而且效率肯定很低:上面的 SIMD 版本实际上比初始标量版本慢三倍(当然,在 Microsoft VS15 的发布模式下进行了优化,并经过 100 万次迭代,而不仅仅是 12 次)。
您不需要中间数组来加载到 SSE 寄存器。只需直接从您的阵列加载。
auto loaded_a = _mm_loadu_ps(&a[i]);
auto loaded_b = _mm_loadu_ps(&b[i]);
_mm_storeu_ps(&c[i], _mm_add_ps(loaded_a, loaded_b));
您也可以省略两个 loaded
变量并将它们合并到 add 中,尽管编译器应该为您完成。
你需要注意这一点,因为如果向量大小不是 4 的倍数,它将无法正常工作(你将访问数组的末尾,导致未定义的行为,并且写入超过 c
的结尾可能会造成破坏)。
您的 for 循环可以简化为
const int aligendN = N - N % 4;
for (int i = 0; i < alignedN; i+=4) {
_mm_storeu_ps(&c[i],
_mm_add_ps(_mm_loadu_ps(&a[i]),
_mm_loadu_ps(&b[i])));
}
for (int i = alignedN; i < N; ++i) {
c[i] = a[i] + b[i];
}
补充说明:
- 处理最后几个浮点数的小循环很常见,当
N%4 != 0
或 N 在编译时未知时,这是强制性的。 - 我注意到您选择了未对齐的版本 load/store,与对齐版本相比,惩罚很小。我在 Whosebug 上找到了这个 link:Is the SSE unaligned load intrinsic any slower than the aligned load intrinsic on x64_64 Intel CPUs?