如何将 AVX/SIMD 与嵌套循环和 += 格式一起使用？
How to use AVX/SIMD with nested loops and += format?
我正在编写网页排名程序。我正在编写一种更新排名的方法。我已经成功地使用嵌套 for 循环和线程版本。但是我想改用 SIMD/AVX.
这是我想更改为 SIMD/AVX 实现的代码。
#define IDX(a, b) ((a * npages) + b) // 2D matrix indexing
for (size_t i = 0; i < npages; i++) {
temp[i] = 0.0;
for (size_t j = 0; j < npages; j++) {
temp[i] += P[j] * matrix_cap[IDX(i,j)];
}
}
对于此代码，P[] 的大小为 npages，matrix_cap[] 的大小为 npages * npages。 P[] 是页面的排名，temp[] 用于存储下一次迭代的页面排名，以便能够检查收敛。
我不知道如何用 AVX 解释 += 以及如何获取我的数据，其中涉及大小为 npages 的两个 arrays/vectors 和一个大小为 [=16] 的矩阵=]（按行主要顺序）转换为可与 SIMD/AVX 操作一起使用的格式。
就 AVX 而言，这是我目前所拥有的，尽管它非常非常不正确，只是我大致想要做的事情的尝试。
ssize_t g_mod = npages - (npages % 4);
double* res = malloc(sizeof(double) * npages);
double sum = 0.0;
for (size_t i = 0; i < npages; i++) {
for (size_t j = 0; j < mod; j += 4) {
__m256d p = _mm256_loadu_pd(P + j);
__m256d m = _mm256_loadu_pd(matrix_hat + i + j);
__m256d pm = _mm256_mul_pd(p, m);
_mm256_storeu_pd(&res + j, pm);
for (size_t k = 0; k < 4; k++) {
sum += res[j + k];
}
}
for (size_t i = mod; i < npages; i++) {
for (size_t j = 0; j < npages; j++) {
sum += P[j] * matrix_cap[IDX(i,j)];
}
}
temp[i] = sum;
sum = 0.0;
}
如何格式化我的数据，以便我可以对其使用 AVX/SIMD 操作（加、乘）来优化它，因为它会被调用很多。
首先，EOF 是正确的，您应该看看 gcc/clang/icc 在自动向量化您的标量代码方面做得如何。我无法为您检查，因为您只发布了代码片段，而不是我可以在 http://gcc.godbolt.org/.
上发布的任何内容
您绝对不需要 malloc 任何东西。请注意，您的内在函数版本仅在 res[] 时使用 32B，并且始终覆盖之前的任何内容。因此，您还不如使用单个 32B 数组。或者更好，use a better method to get a horizontal sum of your vector.
（有关矩阵不同数据排列的建议，请参阅底部）
计算每个 temp[i] 使用每个 P[j]，因此更智能地进行矢量化实际上可以获得一些好处。 对于来自 P[j] 的每个负载，使用该矢量与来自 matrix_cap[] 的 4 个不同负载用于那个 j，但有 4 个不同的 i 值。 您将累积 4 个不同的向量，最后必须将它们中的每一个相加到一个 temp[i] 值。
因此您的内部循环将有 5 个读取流（P[] 和 matrix_cap 的 4 个不同行）。它将进行 4 次水平求和，最后进行 4 次标量存储，最终结果为 4 个连续的 i 值。 （或者可能做两次洗牌和两次 16B 存储）。 （或者 ，这实际上是昂贵的 _mm256_hadd_pd (vhaddpd) 指令的洗牌能力的一个很好的用例，但要小心它的通道内操作）
并行累积 8 到 12 个 temp[i] 值可能会更好，因此来自 P[j] 的每个负载都会重复使用 8 到 12 次。 （尽管如此，请检查编译器输出以确保您没有 运行 超出向量规则并将 __m256d 向量溢出到内存中。）这将为清理循环留下更多工作。
FMA 吞吐量和延迟使得您需要 10 个矢量累加器来保持 10 个 FMA 处于飞行状态，以使 Haswell 上的 FMA 单元饱和。 Skylake 将延迟降低到 4c，因此您只需要 8 个向量累加器即可在 SKL 上使其饱和。 （参见 x86 标签 wiki）。即使您遇到内存瓶颈，而不是执行端口吞吐量瓶颈，您也会需要多个累加器，但它们可能都是相同的 temp[i]（因此您将它们垂直相加为一个向量，然后 hsum 那).
但是，一次累积多个 temp[i] 的结果具有加载后多次重复使用 P[j] 的巨大优势。您还可以在最后保存垂直添加。多个读取流实际上可能有助于隐藏任何一个流中缓存未命中的延迟。 （Intel CPU 中的硬件预取器可以跟踪每 4k 页一个正向流和一个反向流，IIRC）。您可能会取得平衡，并为 4 个 temp[i] 结果中的每一个并行使用两个或三个向量累加器，如果您发现多个读取流是一个问题，但这意味着您必须加载相同的 P[j] 总共更多次。
所以你应该做类似的事情
#define IDX(a, b) ((a * npages) + b) // 2D matrix indexing
for (size_t i = 0; i < (npages & (~7ULL)); i+=8) {
__m256d s0 = _mm256_setzero_pd(),
s1 = _mm256_setzero_pd(),
s2 = _mm256_setzero_pd(),
...
s7 = _mm256_setzero_pd(); // 8 accumulators for 8 i values
for (size_t j = 0; j < (npages & ~(3ULL)); j+=4) {
__m256d Pj = _mm256_loadu_pd(P+j); // reused 8 times after loading
//temp[i] += P[j] * matrix_cap[IDX(i,j)];
s0 = _mm256_fmadd_pd(Pj, _mm256_loadu_pd(&matrix_cap[IDX(i+0,j)]), s0);
s1 = _mm256_fmadd_pd(Pj, _mm256_loadu_pd(&matrix_cap[IDX(i+1,j)]), s1);
// ...
s7 = _mm256_fmadd_pd(Pj, _mm256_loadu_pd(&matrix_cap[IDX(i+7,j)]), s7);
}
// or do this block with a hsum+transpose and do vector stores.
// taking advantage of the power of vhaddpd to be doing 4 useful hsums with each instructions.
temp[i+0] = hsum_pd256(s0); // See the horizontal-sum link earlier for how to write this function
temp[i+1] = hsum_pd256(s1);
//...
temp[i+7] = hsum_pd256(s7);
// if npages isn't a multiple of 4, add the last couple scalar elements to the results of the hsum_pd256()s.
}
// TODO: cleanup for the last up-to-7 odd elements.
您可能会编写 __m256d sums[8] 并遍历您的向量累加器，但您必须检查编译器是否完全展开它并且实际上仍将所有内容保存在寄存器中。
How to can I format my data so I can use AVX/SIMD operations (add,mul) on it to optimise it as it will be called a lot.
我之前错过了这部分问题。首先，显然 float 会给你每个向量（和每单位内存带宽）的元素数量的 2 倍。与缓存命中率增加相比，内存/缓存占用空间减少 2 倍可能会带来更多加速。
理想情况下，矩阵为 "striped" 以匹配向量宽度 。对于 4 个相邻的 i 值，来自矩阵的每个负载都会得到一个 matrix_cap[IDX(i,j)] 的向量，但是对于相同的 4 个 i 值，下一个 32B 将是下一个 j 值。这意味着每个向量累加器都在为每个元素中的不同 i 累加总和，因此最后不需要水平总和。
P[j] 保持线性，但你广播加载它的每个元素，用于每个 4 个 i 值的 8 个向量（或 8 个 is 的 8 个向量float)。因此，您将 P[j] 负载的重用因子增加了向量宽度的一个因子。广播负载在 Haswell 和更高版本上几乎是免费的（仍然只采用加载端口 uop），并且在 SnB/IvB 上便宜很多，他们也采用 shuffle-port uop。
考虑对最内层循环使用OpenMP4.x#pragma omp simd reduction .请记住，omp 缩减不适用于 C++ 数组，因此您必须使用如下所示的临时缩减变量。
#define IDX(a, b) ((a * npages) + b) // 2D matrix indexing
for (size_t i = 0; i < npages; i++) {
my_type tmp_reduction = 0.0; // was: // temp[i] = 0.0;
#pragma omp simd reduction (+:tmp_reduction)
for (size_t j = 0; j < npages; j++) {
tmp_reduction += P[j] * matrix_cap[IDX(i,j)];
}
temp[i] = tmp_reduction;
}
对于 x86 平台，OpenMP4.x 目前由最新的 GCC (4.9+) 和 Intel 编译器支持。一些 LLVM 和 PGI 编译器可能也支持它。
P.S。 Auto-vectorization（"auto" 意味着编译器在没有任何编译指示的情况下进行矢量化，即没有来自开发人员的明确指导）可能 有时 对某些编译器有效变体（尽管由于数组元素作为缩减变量而不太可能）。但是，严格来说自动向量化此代码是不正确的。您必须使用显式 SIMD pragma 来 "resolve" 减少依赖性和（作为一个好的副作用）消除指针歧义（以防通过指针访问数组）。
我正在编写网页排名程序。我正在编写一种更新排名的方法。我已经成功地使用嵌套 for 循环和线程版本。但是我想改用 SIMD/AVX.
这是我想更改为 SIMD/AVX 实现的代码。
#define IDX(a, b) ((a * npages) + b) // 2D matrix indexing
for (size_t i = 0; i < npages; i++) {
temp[i] = 0.0;
for (size_t j = 0; j < npages; j++) {
temp[i] += P[j] * matrix_cap[IDX(i,j)];
}
}
对于此代码，P[] 的大小为 npages，matrix_cap[] 的大小为 npages * npages。 P[] 是页面的排名，temp[] 用于存储下一次迭代的页面排名，以便能够检查收敛。
我不知道如何用 AVX 解释 += 以及如何获取我的数据，其中涉及大小为 npages 的两个 arrays/vectors 和一个大小为 [=16] 的矩阵=]（按行主要顺序）转换为可与 SIMD/AVX 操作一起使用的格式。
就 AVX 而言，这是我目前所拥有的，尽管它非常非常不正确，只是我大致想要做的事情的尝试。
ssize_t g_mod = npages - (npages % 4);
double* res = malloc(sizeof(double) * npages);
double sum = 0.0;
for (size_t i = 0; i < npages; i++) {
for (size_t j = 0; j < mod; j += 4) {
__m256d p = _mm256_loadu_pd(P + j);
__m256d m = _mm256_loadu_pd(matrix_hat + i + j);
__m256d pm = _mm256_mul_pd(p, m);
_mm256_storeu_pd(&res + j, pm);
for (size_t k = 0; k < 4; k++) {
sum += res[j + k];
}
}
for (size_t i = mod; i < npages; i++) {
for (size_t j = 0; j < npages; j++) {
sum += P[j] * matrix_cap[IDX(i,j)];
}
}
temp[i] = sum;
sum = 0.0;
}
如何格式化我的数据，以便我可以对其使用 AVX/SIMD 操作（加、乘）来优化它，因为它会被调用很多。
首先，EOF 是正确的，您应该看看 gcc/clang/icc 在自动向量化您的标量代码方面做得如何。我无法为您检查，因为您只发布了代码片段，而不是我可以在 http://gcc.godbolt.org/.
上发布的任何内容您绝对不需要 malloc 任何东西。请注意，您的内在函数版本仅在 res[] 时使用 32B，并且始终覆盖之前的任何内容。因此，您还不如使用单个 32B 数组。或者更好，use a better method to get a horizontal sum of your vector.
（有关矩阵不同数据排列的建议，请参阅底部）
计算每个 temp[i] 使用每个 P[j]，因此更智能地进行矢量化实际上可以获得一些好处。 对于来自 P[j] 的每个负载，使用该矢量与来自 matrix_cap[] 的 4 个不同负载用于那个 j，但有 4 个不同的 i 值。 您将累积 4 个不同的向量，最后必须将它们中的每一个相加到一个 temp[i] 值。
因此您的内部循环将有 5 个读取流（P[] 和 matrix_cap 的 4 个不同行）。它将进行 4 次水平求和，最后进行 4 次标量存储，最终结果为 4 个连续的 i 值。 （或者可能做两次洗牌和两次 16B 存储）。 （或者 _mm256_hadd_pd (vhaddpd) 指令的洗牌能力的一个很好的用例，但要小心它的通道内操作）
并行累积 8 到 12 个 temp[i] 值可能会更好，因此来自 P[j] 的每个负载都会重复使用 8 到 12 次。 （尽管如此，请检查编译器输出以确保您没有 运行 超出向量规则并将 __m256d 向量溢出到内存中。）这将为清理循环留下更多工作。
FMA 吞吐量和延迟使得您需要 10 个矢量累加器来保持 10 个 FMA 处于飞行状态，以使 Haswell 上的 FMA 单元饱和。 Skylake 将延迟降低到 4c，因此您只需要 8 个向量累加器即可在 SKL 上使其饱和。 （参见 x86 标签 wiki）。即使您遇到内存瓶颈，而不是执行端口吞吐量瓶颈，您也会需要多个累加器，但它们可能都是相同的 temp[i]（因此您将它们垂直相加为一个向量，然后 hsum 那).
但是，一次累积多个 temp[i] 的结果具有加载后多次重复使用 P[j] 的巨大优势。您还可以在最后保存垂直添加。多个读取流实际上可能有助于隐藏任何一个流中缓存未命中的延迟。 （Intel CPU 中的硬件预取器可以跟踪每 4k 页一个正向流和一个反向流，IIRC）。您可能会取得平衡，并为 4 个 temp[i] 结果中的每一个并行使用两个或三个向量累加器，如果您发现多个读取流是一个问题，但这意味着您必须加载相同的 P[j] 总共更多次。
所以你应该做类似的事情
#define IDX(a, b) ((a * npages) + b) // 2D matrix indexing
for (size_t i = 0; i < (npages & (~7ULL)); i+=8) {
__m256d s0 = _mm256_setzero_pd(),
s1 = _mm256_setzero_pd(),
s2 = _mm256_setzero_pd(),
...
s7 = _mm256_setzero_pd(); // 8 accumulators for 8 i values
for (size_t j = 0; j < (npages & ~(3ULL)); j+=4) {
__m256d Pj = _mm256_loadu_pd(P+j); // reused 8 times after loading
//temp[i] += P[j] * matrix_cap[IDX(i,j)];
s0 = _mm256_fmadd_pd(Pj, _mm256_loadu_pd(&matrix_cap[IDX(i+0,j)]), s0);
s1 = _mm256_fmadd_pd(Pj, _mm256_loadu_pd(&matrix_cap[IDX(i+1,j)]), s1);
// ...
s7 = _mm256_fmadd_pd(Pj, _mm256_loadu_pd(&matrix_cap[IDX(i+7,j)]), s7);
}
// or do this block with a hsum+transpose and do vector stores.
// taking advantage of the power of vhaddpd to be doing 4 useful hsums with each instructions.
temp[i+0] = hsum_pd256(s0); // See the horizontal-sum link earlier for how to write this function
temp[i+1] = hsum_pd256(s1);
//...
temp[i+7] = hsum_pd256(s7);
// if npages isn't a multiple of 4, add the last couple scalar elements to the results of the hsum_pd256()s.
}
// TODO: cleanup for the last up-to-7 odd elements.
您可能会编写 __m256d sums[8] 并遍历您的向量累加器，但您必须检查编译器是否完全展开它并且实际上仍将所有内容保存在寄存器中。
How to can I format my data so I can use AVX/SIMD operations (add,mul) on it to optimise it as it will be called a lot.
我之前错过了这部分问题。首先，显然 float 会给你每个向量（和每单位内存带宽）的元素数量的 2 倍。与缓存命中率增加相比，内存/缓存占用空间减少 2 倍可能会带来更多加速。
理想情况下，矩阵为 "striped" 以匹配向量宽度 。对于 4 个相邻的 i 值，来自矩阵的每个负载都会得到一个 matrix_cap[IDX(i,j)] 的向量，但是对于相同的 4 个 i 值，下一个 32B 将是下一个 j 值。这意味着每个向量累加器都在为每个元素中的不同 i 累加总和，因此最后不需要水平总和。
P[j] 保持线性，但你广播加载它的每个元素，用于每个 4 个 i 值的 8 个向量（或 8 个 is 的 8 个向量float)。因此，您将 P[j] 负载的重用因子增加了向量宽度的一个因子。广播负载在 Haswell 和更高版本上几乎是免费的（仍然只采用加载端口 uop），并且在 SnB/IvB 上便宜很多，他们也采用 shuffle-port uop。
考虑对最内层循环使用OpenMP4.x#pragma omp simd reduction .请记住，omp 缩减不适用于 C++ 数组，因此您必须使用如下所示的临时缩减变量。
#define IDX(a, b) ((a * npages) + b) // 2D matrix indexing
for (size_t i = 0; i < npages; i++) {
my_type tmp_reduction = 0.0; // was: // temp[i] = 0.0;
#pragma omp simd reduction (+:tmp_reduction)
for (size_t j = 0; j < npages; j++) {
tmp_reduction += P[j] * matrix_cap[IDX(i,j)];
}
temp[i] = tmp_reduction;
}
对于 x86 平台，OpenMP4.x 目前由最新的 GCC (4.9+) 和 Intel 编译器支持。一些 LLVM 和 PGI 编译器可能也支持它。
P.S。 Auto-vectorization（"auto" 意味着编译器在没有任何编译指示的情况下进行矢量化，即没有来自开发人员的明确指导）可能 有时 对某些编译器有效变体（尽管由于数组元素作为缩减变量而不太可能）。但是，严格来说自动向量化此代码是不正确的。您必须使用显式 SIMD pragma 来 "resolve" 减少依赖性和（作为一个好的副作用）消除指针歧义（以防通过指针访问数组）。