与SSE并行前缀（累计）总和

Question

我正在寻找一些关于如何与SSE做一个并行前缀总和的建议。我有兴趣通过一系列整数，浮点数或双精度来做这件事。

我已经想出了两个解决方案。特例和一般情况。在这两种情况下，解决方案都以两次与OpenMP并行的方式在阵列上运行。对于特殊情况，我在两次通行证上都使用SSE。对于一般情况，我只在第二阶段使用它。

我的主要问题是我如何在一般情况下的第一次使用SSE？以下链接simd-prefix-sum-on-intel-cpu显示字节的改进，但不显示32位数据类型。

特殊情况被称为特殊的原因是它需要数组是特殊的格式。例如，让我们假设浮点数只有16个元素的数组a。然后，如果阵列重新排列是这样的（结构的数组为结构阵列）：

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; 
} 

Answer 1

这是我第一次回答我的问题，但它似乎是适当的。基于hirschhornsalz 答案前缀总和在16个字节simd-prefix-sum-on-intel-cpu我已经想出了一个解决方案，使用SIMD在第一遍4,8和16 32位字。

一般理论如下。对于n单词的顺序扫描，需要n添加（n-1扫描n个单词，并且从先前扫描的单词集中执行另外一个添加）。然而，使用SIMD n个字可以在log （n）中进行扫描，并且相等数量的移位再加上另外一个加法和广播以从先前的SIMD扫描携带。所以对于一些值为n的SIMD方法会赢。

让我们看一下使用SSE，AVX和32位字AVX-512：

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不会为32位字的扫描，直到AVX-有用512。这也假定只能在一条指令中完成转换和广播。这对于上海证券交易所是如此，但not for AVX and maybe not even for 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阵列中的主循环中使用另一个加法和一个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); 
    } 
} 

三班倒需要7个内在。广播需要2种内在因素。所以，这4个增加了13个内在因素。对于AVX2，只需要5个内在函数就可以完成11个内部函数。顺序总和只需要8个添加。因此AVX和AVX2都不适用于第一遍。

编辑：

所以最后这个基准测试，结果是出乎意料的。上证所和AVX指令都约快两倍以下顺序代码：

void scan(float a[], float s[], int n) { 
    float sum = 0; 
    for (int i = 0; i<n; i++) { 
     sum += a[i]; 
     s[i] = sum; 
    } 
} 

我想这是由于指令级paralellism。

这样回答我自己的问题。在一般情况下，我成功地将SIMD用于pass1。当我在我的4核心常春藤桥系统上将它与OpenMP结合起来时，对于512k浮点数，总加速大约为7。