我正在编写一个页面排名程序。我正在编写一种更新排名的方法。我已经成功地让它与嵌套的循环和线程版本一起工作。但是，我想改用 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的数组/向量和一个大小为npages * npages的矩阵(按行主要顺序)的数据转换为可用于 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 操作 (add，mul) 来优化它，因为它会被调用很多。
考虑使用OpenMP4.x#pragma omp simd reduce 进行最内层循环。请记住，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+) 和英特尔编译器支持。一些LLVM和PGI编译器也可能支持它。
附言自动矢量化("auto"是指编译器的矢量化，没有任何编译指示，即没有开发人员的明确指导)有时可能适用于某些编译器变体(尽管由于数组元素作为归约变量，这不太可能)。但是，严格来说，自动矢量化此代码是不正确的。您必须使用显式 SIMD 编译指示来"解决"缩减依赖性，并(作为一个好的副作用)消除指针的歧义(如果数组是通过指针访问的)。
首先，EOF是正确的，你应该看到gcc/clang/icc在自动矢量化标量代码方面做得如何。 我无法为您检查，因为您只发布了代码片段，而不是我可以扔到 http://gcc.godbolt.org/上的任何内容。
你绝对不需要恶意处理任何东西。 请注意，您的内部函数版本仅在res[]时使用 32B，并且总是覆盖之前的任何内容。 因此，您不妨使用单个 32B 阵列。 或者更好的是，使用更好的方法来获得向量的水平和。
(有关矩阵的不同数据排列的建议，请参见底部)
计算每个temp[i]会使用每个P[j]，因此实际上从更智能地矢量化中可以获得一些好处。对于来自P[j]的每个载荷，使用该向量与该j的matrix_cap[]有4个不同的载荷，但有4个不同的i值。您将累积 4 个不同的向量，并且必须在最后将它们中的每一个相加到一个temp[i]值。
所以你的内部循环将有 5 个读取流(P[] 和 4 个不同的matrix_cap行)。 它将做 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 标记维基)。 即使你在内存上遇到瓶颈，而不是执行端口吞吐量，你也需要多个累加器，但它们都可以用于相同的temp[i](所以你可以将它们垂直相加为一个向量，然后求和)。
但是，一次累积多个temp[i]的结果具有在加载后多次重用P[j]的巨大优势。 您还可以在末尾保存垂直添加。 多个读取流实际上可能有助于隐藏任何一个流中缓存未命中的延迟。 (英特尔 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]并循环你的向量累加器，但你必须检查编译器是否完全展开它，并且仍然将所有内容保存在寄存器中。
如何格式化数据，以便对数据使用 AVX/SIMD 操作 (add，mul) 来优化它，因为它会被调用很多。
我之前错过了问题的这一部分。 首先，显然float会给你 2 倍的元素数量 每个向量(和每单位内存带宽)。 如果缓存命中率增加，内存/缓存占用量减少 2 的系数可能会比这提供更多的加速。
理想情况下，矩阵将被"条带化"以匹配矢量宽度。 来自矩阵的每个负载将获得 4 个相邻i值的matrix_cap[IDX(i,j)]向量，但下一个 32B 将是相同 4 个i值的下一个j值。 这意味着每个向量累加器正在累积每个元素中不同i的总和，因此最后不需要水平和。
P[j]保持线性，但您广播加载它的每个元素，用于 8 个向量，每个向量 4 个i值(或 8 个 vec 的 8 个i秒用于float)。 因此，您将P[j]载荷的重用系数增加矢量宽度的系数。 广播负载在 Haswell 和以后几乎是免费的(仍然只采用加载端口 uop)，并且在 SnB/IvB 上非常便宜，它们也采用随机端口 uop。