與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。