与SSE并行的前缀(累计)总和
我正在寻找有关如何与SSE进行并行前缀和的一些建议。 我有兴趣在一系列的整数,浮点数或双精度数上做这个。
我想出了两个解决方案。 一个特例和一般情况。 在这两种情况下,解决方案在与OpenMP并行的两次传递中在arrays上运行。 对于特殊情况,我在两次传球时使用SSE。 对于一般情况,我只在第二遍使用它。
我的主要问题是如何在一般情况下的第一遍中使用SSE? 以下链接simd-prefix-sum-on-intel-cpu显示字节的改进,但不是32位数据类型。
特殊情况称为特殊情况的原因是它要求数组采用特殊格式。 例如,让我们假设浮点数组中只有16个元素。 然后,如果数组像这样重新排列(结构数组结构):
a[0] a[1] ...a[15] -> a[0] a[4] a[8] a[12] a[1] a[5] a[9] a[13]...a[3] a[7] a[11] a[15]
SSE垂直总和可用于两个通道。 但是,只有当数组已经采用特殊格式并且输出可以以特殊格式使用时,这才有效。 否则,必须在输入和输出上进行昂贵的重新排列,这将使其比一般情况慢得多。
也许我应该考虑一个不同的前缀和算法(例如二叉树)?
一般情况的代码:
void prefix_sum_omp_sse(double a[], double s[], int n) { double *suma; #pragma omp parallel { const int ithread = omp_get_thread_num(); const int nthreads = omp_get_num_threads(); #pragma omp single { suma = new double[nthreads + 1]; suma[0] = 0; } double sum = 0; #pragma omp for schedule(static) nowait //first parallel pass for (int i = 0; i<n; i++) { sum += a[i]; s[i] = sum; } suma[ithread + 1] = sum; #pragma omp barrier #pragma omp single { double tmp = 0; for (int i = 0; i<(nthreads + 1); i++) { tmp += suma[i]; suma[i] = tmp; } } __m128d offset = _mm_set1_pd(suma[ithread]); #pragma omp for schedule(static) //second parallel pass with SSE as well for (int i = 0; i<n/4; i++) { __m128d tmp1 = _mm_load_pd(&s[4*i]); tmp1 = _mm_add_pd(tmp1, offset); __m128d tmp2 = _mm_load_pd(&s[4*i+2]); tmp2 = _mm_add_pd(tmp2, offset); _mm_store_pd(&s[4*i], tmp1); _mm_store_pd(&s[4*i+2], tmp2); } } delete[] suma; }
这是我第一次回答我自己的问题,但这似乎是合适的。 基于hirschhornsalz回答16字节simd-prefix-sum-on-intel-cpu上的前缀和我得到了一个解决方案,在第一次使用4,8和16个32位字时使用SIMD。
一般理论如下。 对于n
字的顺序扫描,需要n
加法(n-1扫描n个字,并且从扫描的前一组字中携带另外一个加法)。 然而,使用SIMD n字可以在log 2 (n)加法中扫描,并且相同数量的移位加上一次加法和广播以从先前的SIMD扫描进行携带。 因此,对于n
的某个值,SIMD方法将获胜。
让我们看看SSE,AVX和AVX-512的32位字:
4 32-bit words (SSE): 2 shifts, 3 adds, 1 broadcast sequential: 4 adds 8 32-bit words (AVX): 3 shifts, 4 adds, 1 broadcast sequential: 8 adds 16 32 bit-words (AVX-512): 4 shifts, 5 adds, 1 broadcast sequential: 16 adds
基于此,看起来SIMD在AVX-512之前对32位字的扫描没有用。 这也假设移位和广播只能在1条指令中完成。 这适用于SSE但不适用于AVX,甚至可能不适用于AVX2 。
无论如何,我把一些工作和测试的代码放在一起,使用SSE进行前缀和。
inline __m128 scan_SSE(__m128 x) { x = _mm_add_ps(x, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(x), 4))); x = _mm_add_ps(x, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(x), 8))); return x; } void prefix_sum_SSE(float *a, float *s, const int n) { __m128 offset = _mm_setzero_ps(); for (int i = 0; i < n; i+=4) { __m128 x = _mm_load_ps(&a[i]); __m128 out = scan_SSE(x); out = _mm_add_ps(out, offset); _mm_store_ps(&s[i], out); offset = _mm_shuffle_ps(out, out, _MM_SHUFFLE(3, 3, 3, 3)); }
请注意, scan_SSE
函数有两个加法(_mm_add_ps)和两个移位(_mm_slli_si128)。 强制转换仅用于使编译器满意并且不会转换为指令。 然后在prefix_sum_SSE
的数组上的主循环内, prefix_sum_SSE
另一个加法和一个shuffle。 总共有6次操作,而顺序总和只增加了4次。
这是AVX的工作解决方案:
inline __m256 scan_AVX(__m256 x) { __m256 t0, t1; //shift1_AVX + add t0 = _mm256_permute_ps(x, _MM_SHUFFLE(2, 1, 0, 3)); t1 = _mm256_permute2f128_ps(t0, t0, 41); x = _mm256_add_ps(x, _mm256_blend_ps(t0, t1, 0x11)); //shift2_AVX + add t0 = _mm256_permute_ps(x, _MM_SHUFFLE(1, 0, 3, 2)); t1 = _mm256_permute2f128_ps(t0, t0, 41); x = _mm256_add_ps(x, _mm256_blend_ps(t0, t1, 0x33)); //shift3_AVX + add x = _mm256_add_ps(x,_mm256_permute2f128_ps(x, x, 41)); return x; } void prefix_sum_AVX(float *a, float *s, const int n) { __m256 offset = _mm256_setzero_ps(); for (int i = 0; i < n; i += 8) { __m256 x = _mm256_loadu_ps(&a[i]); __m256 out = scan_AVX(x); out = _mm256_add_ps(out, offset); _mm256_storeu_ps(&s[i], out); //broadcast last element __m256 t0 = _mm256_permute2f128_ps(out, out, 0x11); offset = _mm256_permute_ps(t0, 0xff); } }
这三个class次需要7个内在函数。 广播需要2个内在函数。 因此,有13个内在因素的4个增加。 对于AVX2,转换只需要5个内在函数,因此总共有11个内在函数。 顺序总和只需要8次加法。 因此,AVX和AVX2可能都不适用于第一次传递。
编辑:
所以我最终对此进行了基准测试,结果出乎意料。 SSE和AVX代码的速度都是以下顺序代码的两倍:
void scan(float a[], float s[], int n) { float sum = 0; for (int i = 0; i
我想这是由于指令级别的并行性。
所以这回答了我自己的问题。 在一般情况下,我成功地将SIMD用于pass1。 当我在我的4核常春藤桥系统上将它与OpenMP结合使用时,512k浮点数的总速度约为7。