代码之家 › 专栏 › 技术社区 › joshuatvernon imotov

如何在嵌套循环和+=格式中使用AVX/SIMD?

pagerank avx simd multithreading c

1

joshuatvernon imotov · 技术社区 · 9 年前

我正在写一个页面排名程序。我正在写一个更新排名的方法。我已经成功地将其用于嵌套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以及如何获取包含两个大小数组/向量的数据 n页 和一个尺寸矩阵 n页*n页 (按行主要顺序)转换为可用于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)对其进行优化,因为它将被称为很多。

2 回复 | 直到 9 年前

1

3

zam 9 年前

考虑使用 OpenMP4.x版 #pragma omp-simd约简 用于最里面的循环。请记住,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”是指编译器在没有任何杂注的情况下进行矢量化,即没有开发人员的明确指示)可以有时适用于某些编译器变体(尽管由于数组元素是reduction变量,所以不太可能)。然而,严格来说,自动矢量化此代码是不正确的。您必须使用显式SIMD pragma来“解决”缩减依赖关系,并消除指针的歧义(作为一个好的副作用)(以防通过指针访问数组)。

2

Community CDub 8 年前

首先,EOF是正确的,您应该看看gcc/clang/icc在自动矢量化标量代码方面做得有多好。我不能帮你检查,因为你只发布了代码片段,没有我能提供的任何东西 http://gcc.godbolt.org/ .

你绝对不需要分配任何东西。请注意,您的intrinsic版本在 res[] ,并始终覆盖以前存在的内容。因此,您不妨使用单个32B阵列。或者更好, use a better method to get a horizontal sum of your vector .

(有关矩阵不同数据安排的建议,请参见底部)

计算每个 temp[i] 使用频率 P[j] ,所以更聪明地进行矢量化实际上可以获得一些好处。 对于来自的每个负载 P【j】 ,将该矢量用于4种不同的载荷 matrix_cap[] 为了那个 j ,但4个不同 i 价值观 您将累积4个不同的矢量,并且必须将每个矢量加和为 温度[i] 值。

因此,您的内部循环将有5个读取流(P[]和4个不同的行 matrix_cap ). 它将进行4次水平求和,并在末尾进行4次标量存储,最终结果为连续4次 我 价值观(或者做两次洗牌和两个16B商店)。(或者可能 transpose-and-sum together ,这实际上是昂贵的 _mm256_hadd_pd ( vhaddpd )指示,但小心其车道内操作)

积累8到12可能更好 温度[i] 值并行,因此来自 P【j】 重复使用8至12次。(检查编译器输出,确保没有用完矢量reg并溢出 __m256d 但是,内存向量。)这将为清理循环留下更多工作。

FMA吞吐量和延迟使得您需要10个矢量累加器来保持10个FMA处于飞行状态,以使哈斯韦尔上的FMA单元饱和。Skylake将延迟降低到4c,因此您只需要8个矢量累加器就可以在SKL上饱和它。(请参见 x86 标签维基)。即使您的内存瓶颈,而不是执行端口吞吐量瓶颈,您也需要多个累加器,但它们可能都是相同的 温度[i] (所以你可以把它们垂直加起来,得到一个向量,然后再加上这个)。

然而,累积多个 温度[i] 同时具有重用的巨大优势 P【j】 加载后多次。还可以在末尾保存垂直添加。多个读取流实际上可能有助于隐藏任何一个流中缓存未命中的延迟。(Intel CPU中的硬件预取器可以跟踪每4k页一个正向流和一个反向流,IIRC)。您可以平衡一下,对4个向量累加器中的每一个使用两个或三个向量累加器 温度[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倍的因素可能会带来更大的加速。

理想情况下,矩阵将被“分条”以匹配矢量宽度 。矩阵中的每个加载都会得到一个向量 matrix_cap[IDX(i,j)] 4个相邻 我 值,但下一个32B将是下一个 j 相同4的值 我 价值观这意味着每个矢量累加器正在为不同的 我 在每个元素中,因此不需要在末尾进行水平求和。

P【j】 保持线性,但广播加载它的每个元素,用于4的8个矢量 我 每个值(或8个vec,共8个 我 s用于 浮动 ). 因此,您增加了 P【j】 按矢量宽度系数加载。在Haswell和以后的版本中,广播加载几乎是免费的(仍然只使用加载端口uop),而在SnB/IvB上,这一点非常便宜,因为它们也使用洗牌端口uop。