インテルCPU上のSIMDプリフィックスの合計

Question

プリフィックスサムアルゴリズムを実装する必要があり、できるだけ高速にする必要があります。 例：

[3, 1, 7, 0, 4, 1, 6, 3] 
should give 
[3, 4, 11, 11, 15, 16, 22, 25] 

この使用SSE/MMX/SIMD CPUの命令を実行する方法はありますか？

私の最初のアイデアは、すべての合計が以下のように計算されるまで、各ペアを並行して再帰的に合計することです！

 //in parallel do 
     for (int i = 0; i<z.length; i++){ 
      z[i] = x[i<<1] + x[(i<<1)+1]; 
     } 

アルゴリズムもう少し明確な「Z」は、最終ouputを

ではなく、代わりにouputを

 int[] w = computePrefixSum(z); 
     for (int i = 1; i<ouput.length; i++){ 
      ouput[i] = (i%2==0) ? (x[i] + ouput[i-1]) : w[(i-1)>>1]; 
     } 

Answer 1

私の知っている最速の並列プレフィックス合計アルゴリズムを並列に2回のパスで合計上で実行され、第二のパスでも同様にSSEを使用することです。

最初のパスでは、部分和を並列に計算し、各部分合計の合計を格納します。 2回目のパスでは、前の部分合計から次の部分合計までの合計を加算します。 OpenMPを使用して複数のスレッドを使用して、両方のパスを並行して実行できます。 2番目のパスでは、SIMDを使用することもできます。これは、一定の値が各部分合計に追加されるためです。

拳パスがSIMDを使用しないので、時間コストは常にn/m

より大きくなる時間コストが

n/m + n/(m*w) = (n/m)*(1+1/w) 

なければならないアレイのn要素、mコア、及びwのSIMD幅を仮定

たとえば、SIMD_widthが4の4つのコア（SSEを搭載した4つの32ビット浮動小数点数）の場合、コストは5n/16になります。または、時間コストがnのシーケンシャルコードより約3.2倍高速です。ハイパースレッディングを使用すると、スピードアップはさらに大きくなります。

特殊なケースでは、最初のパスでもSIMDを使用できます。そして、時間コストは、単に

2*n/(m*w) 

である私は、SSEコードのためのスレッドと組み込み関数のためのOpenMPを使用し、以下のリンク parallel-prefix-cumulative-sum-with-sse

で特殊なケースの詳細を議論する一般的なケースのためのコードを掲載しました編集： 私は、シーケンシャルコードの約2倍の速さで最初のパスのSIMDバージョンを見つけることができました。今度は私の4つのコア・アイビー・ブリッジ・システムで約7倍のブーストが得られます。

編集：より大きなアレイの場合 一つの問題は、最初のパスの後にほとんどの値がキャッシュから追い出されていることです。私は、チャンク内で並列に実行されるが、各チャンクを連続して実行するソリューションを考え出した。 chunk_sizeは、調整する必要がある値です。たとえば、1MB = 256Kに設定します。値がまだレベル2のキャッシュ内にある間に、2番目のパスが実行されます。これを行うことで、大きな配列の大きな改善が得られます。

ここにSSEのコードを示します。 AVXコードはほぼ同じ速度ですので、私はここに投稿しませんでした。接頭辞の合計を行う関数はscan_omp_SSEp2_SSEp1_chunkです。 floatの配列aを渡し、配列sを累積合計で満たします。 1000個の32ビット整数の配列のために

__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_shuffle_ps(_mm_setzero_ps(), x, 0x40)); 
    return x; 
} 

float pass1_SSE(float *a, float *s, const int n) { 
    __m128 offset = _mm_setzero_ps(); 
    #pragma omp for schedule(static) nowait 
    for (int i = 0; i < n/4; i++) { 
     __m128 x = _mm_load_ps(&a[4 * i]); 
     __m128 out = scan_SSE(x); 
     out = _mm_add_ps(out, offset); 
     _mm_store_ps(&s[4 * i], out); 
     offset = _mm_shuffle_ps(out, out, _MM_SHUFFLE(3, 3, 3, 3)); 
    } 
    float tmp[4]; 
    _mm_store_ps(tmp, offset); 
    return tmp[3]; 
} 

void pass2_SSE(float *s, __m128 offset, const int n) { 
    #pragma omp for schedule(static) 
    for (int i = 0; i<n/4; i++) { 
     __m128 tmp1 = _mm_load_ps(&s[4 * i]); 
     tmp1 = _mm_add_ps(tmp1, offset); 
     _mm_store_ps(&s[4 * i], tmp1); 
    } 
} 

void scan_omp_SSEp2_SSEp1_chunk(float a[], float s[], int n) { 
    float *suma; 
    const int chunk_size = 1<<18; 
    const int nchunks = n%chunk_size == 0 ? n/chunk_size : n/chunk_size + 1; 
    //printf("nchunks %d\n", nchunks); 
    #pragma omp parallel 
    { 
     const int ithread = omp_get_thread_num(); 
     const int nthreads = omp_get_num_threads(); 

     #pragma omp single 
     { 
      suma = new float[nthreads + 1]; 
      suma[0] = 0; 
     } 

     float offset2 = 0.0f; 
     for (int c = 0; c < nchunks; c++) { 
      const int start = c*chunk_size; 
      const int chunk = (c + 1)*chunk_size < n ? chunk_size : n - c*chunk_size; 
      suma[ithread + 1] = pass1_SSE(&a[start], &s[start], chunk); 
      #pragma omp barrier 
      #pragma omp single 
      { 
       float tmp = 0; 
       for (int i = 0; i < (nthreads + 1); i++) { 
        tmp += suma[i]; 
        suma[i] = tmp; 
       } 
      } 
      __m128 offset = _mm_set1_ps(suma[ithread]+offset2); 
      pass2_SSE(&s[start], offset, chunk); 
      #pragma omp barrier 
      offset2 = s[start + chunk-1]; 
     } 
    } 
    delete[] suma; 
} 

Answer 2

プレフィックス合計を計算するために使用さを並列に計算することができるようにするに実際には、GPUプログラミングの基本的なアルゴリズムの1つです。 Intelプロセッサ上でSIMD拡張を使用しているのであれば、並行して実行すると実際には多くの利点があるのか​​どうかはわかりませんが、nvidiaのこの記事で並列プレフィックス和を実装してみましょう（アルゴリズムを見て無視してください） CUDA）：http://http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html

Answer 3

大きなレジスタの長さと小さな合計に対して、多少の並列性を利用することができます。たとえば、1バイトの16個の値を加算すると（これは1つのsseレジスタに収まる）、ログの数はで、シフト数は同じです。
あまり多くはありませんが、15個以上の追加と追加のメモリアクセスが必要です。

__m128i x = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3); 
x = _mm_add_epi8(x, _mm_srli_si128(x, 1)); 
x = _mm_add_epi8(x, _mm_srli_si128(x, 2)); 
x = _mm_add_epi8(x, _mm_srli_si128(x, 4)); 
x = _mm_add_epi8(x, _mm_srli_si128(x, 8)); 

// x == 3, 4, 11, 11, 15, 16, 22, 25, 28, 29, 36, 36, 40, 41, 47, 50 

合計が長い場合は、命令レベルの並列性を活用し、命令の並べ替えを利用することで依存関係を隠すことができます。

編集：

__m128i x0 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3); 
__m128i x1 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3); 
__m128i x2 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3); 
__m128i x3 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3); 

__m128i mask = _mm_set_epi8(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0); 

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 1)); 
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 1)); 
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 1)); 
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 1)); 

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 2)); 
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 2)); 
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 2)); 
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 2)); 

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 4)); 
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 4)); 
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 4)); 
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 4)); 

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 8)); 
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 8)); 
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 8)); 
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 8)); 

x1 = _mm_add_epi8(_mm_shuffle_epi8(x0, mask), x1); 
x2 = _mm_add_epi8(_mm_shuffle_epi8(x1, mask), x2); 
x3 = _mm_add_epi8(_mm_shuffle_epi8(x2, mask), x3); 

のようなもの

Answer 4

、私はインテルのSandybridge上のループの@ hirschhornsalzの方法を使用して、およそ1.4倍シングルスレッドの小さなスピードアップを取得することができました。 intの60キロバイトのバッファでは、スピードアップは約1.37です。 8MiBの整数では、スピードアップは1.13です。 （DDR3-1600と3.8GHzのターボでi5-2500k、。）

小さい要素（int16_t又はuint8_t、または符号なしのバージョン）は、ベクトル当たりの要素の数のそれぞれ2倍にする追加/シフトの余分な段階を取ります。オーバーフローが悪いので、SSEに大きな利点があるとしても、すべての要素の合計を保持できないデータ型を使用しないでください。

#include <immintrin.h> 

// In-place rewrite an array of values into an array of prefix sums. 
// This makes the code simpler, and minimizes cache effects. 
int prefix_sum_sse(int data[], int n) 
{ 
// const int elemsz = sizeof(data[0]); 
#define elemsz sizeof(data[0]) // clang-3.5 doesn't allow compile-time-const int as an imm8 arg to intrinsics 

    __m128i *datavec = (__m128i*)data; 
    const int vec_elems = sizeof(*datavec)/elemsz; 
    // to use this for int8/16_t, you still need to change the add_epi32, and the shuffle 

    const __m128i *endp = (__m128i*) (data + n - 2*vec_elems); // don't start an iteration beyond this 
    __m128i carry = _mm_setzero_si128(); 
    for(; datavec <= endp ; datavec += 2) { 
     IACA_START 
     __m128i x0 = _mm_load_si128(datavec + 0); 
     __m128i x1 = _mm_load_si128(datavec + 1); // unroll/pipeline by 1 
//  __m128i x2 = _mm_load_si128(datavec + 2); 
//  __m128i x3; 

     x0 = _mm_add_epi32(x0, _mm_slli_si128(x0, elemsz)); // for floats, use shufps not bytewise-shift 
     x1 = _mm_add_epi32(x1, _mm_slli_si128(x1, elemsz)); 

     x0 = _mm_add_epi32(x0, _mm_slli_si128(x0, 2*elemsz)); 
     x1 = _mm_add_epi32(x1, _mm_slli_si128(x1, 2*elemsz)); 

    // more shifting if vec_elems is larger 

     x0 = _mm_add_epi32(x0, carry); // this has to go after the byte-shifts, to avoid double-counting the carry. 
     _mm_store_si128(datavec +0, x0); // store first to allow destructive shuffle (non-avx pshufb if needed) 

     x1 = _mm_add_epi32(_mm_shuffle_epi32(x0, _MM_SHUFFLE(3,3,3,3)), x1); 
     _mm_store_si128(datavec +1, x1); 

     carry = _mm_shuffle_epi32(x1, _MM_SHUFFLE(3,3,3,3)); // broadcast the high element for next vector 
    } 
    // FIXME: scalar loop to handle the last few elements 
    IACA_END 
    return data[n-1]; 
    #undef elemsz 
} 

int prefix_sum_simple(int data[], int n) 
{ 
    int sum=0; 
    for (int i=0; i<n ; i++) { 
     IACA_START 
     sum += data[i]; 
     data[i] = sum; 
    } 
    IACA_END 
    return sum; 
} 

// perl -we '$n=1000; sub rnlist($$) { return map { int rand($_[1]) } (1..$_[0]);} @a=rnlist($n,127); $"=", "; print "$n\n@a\n";' 

int data[] = { 51, 83, 126, 11, 20, 63, 113, 102, 
     126,67, 83, 113, 86, 123, 30, 109, 
     97, 71, 109, 86, 67, 60, 47, 12, 
     /* ... */ }; 


int main(int argc, char**argv) 
{ 
    const int elemsz = sizeof(data[0]); 
    const int n = sizeof(data)/elemsz; 
    const long reps = 1000000 * 1000/n; 
    if (argc >= 2 && *argv[1] == 'n') { 
     for (int i=0; i < reps ; i++) 
      prefix_sum_simple(data, n); 
    }else { 
     for (int i=0; i < reps ; i++) 
      prefix_sum_sse(data, n); 
    } 
    return 0; 
} 

---

バイナリにコンパイルされたリストのn = 1000とテスト、。 （そして、実際にループしていることを確認しました。コンパイル時のショートカットを取らずに、ベクトルまたは非ベクトルテストを無意味にします）。 movdqa命令がたくさんありますが、わずかなサイクルしか保存しません。これは、shuffleとvector-int-addはどちらもSnB/IvBのポート1と5でのみ実行できるため、port0にはmov命令を実行するための余裕のあるサイクルがたくさんあるためです。 uop-cacheのスループットのボトルネックは、非AVXバージョンの方がやや遅い理由かもしれません。 （これらの余分なmov命令は私たちを3.35 insn/cycleまで押し上げる）。フロントエンドはサイクルの4.54％しかアイドル状態ではないので、ほとんど維持していません。

gcc -funroll-loops -DIACA_MARKS_OFF -g -std=c11 -Wall -march=native -O3 prefix-sum.c -mno-avx -o prefix-sum-noavx 

    # gcc 4.9.2 

################# SSE (non-AVX) vector version ############ 
$ ocperf.py stat -e task-clock,cycles,instructions,uops_issued.any,uops_dispatched.thread,uops_retired.all,uops_retired.retire_slots,stalled-cycles-frontend,stalled-cycles-backend ./prefix-sum-noavx 
perf stat -e task-clock,cycles,instructions,cpu/event=0xe,umask=0x1,name=uops_issued_any/,cpu/event=0xb1,umask=0x1,name=uops_dispatched_thread/,cpu/event=0xc2,umask=0x1,name=uops_retired_all/,cpu/event=0xc2,umask=0x2,name=uops_retired_retire_slots/,stalled-cycles-frontend,stalled-cycles-backend ./prefix-sum-noavx 

Performance counter stats for './prefix-sum-noavx': 

     206.986720  task-clock (msec)   # 0.999 CPUs utilized   
     777,473,726  cycles     # 3.756 GHz      
    2,604,757,487  instructions    # 3.35 insns per cycle   
                # 0.01 stalled cycles per insn 
    2,579,310,493  uops_issued_any   # 12461.237 M/sec 
    2,828,479,147  uops_dispatched_thread # 13665.027 M/sec 
    2,829,198,313  uops_retired_all   # 13668.502 M/sec (unfused domain) 
    2,579,016,838  uops_retired_retire_slots # 12459.818 M/sec (fused domain) 
     35,298,807  stalled-cycles-frontend # 4.54% frontend cycles idle 
     1,224,399  stalled-cycles-backend # 0.16% backend cycles idle 

     0.207234316 seconds time elapsed 
------------------------------------------------------------ 


######### AVX (same source, but built with -mavx). not AVX2 ######### 
$ ocperf.py stat -e task-clock,cycles,instructions,uops_issued.any,uops_dispatched.thread,uops_retired.all,uops_retired.retire_slots,stalled-cycles-frontend,stalled-cycles-backend ./prefix-sum-avx 

Performance counter stats for './prefix-sum-avx': 

     203.429021  task-clock (msec)   # 0.999 CPUs utilized   
     764,859,441  cycles     # 3.760 GHz      
    2,079,716,097  instructions    # 2.72 insns per cycle   
                # 0.12 stalled cycles per insn 
    2,054,334,040  uops_issued_any   # 10098.530 M/sec     
    2,303,378,797  uops_dispatched_thread # 11322.764 M/sec     
    2,304,140,578  uops_retired_all   # 11326.509 M/sec     
    2,053,968,862  uops_retired_retire_slots # 10096.735 M/sec     
     240,883,566  stalled-cycles-frontend # 31.49% frontend cycles idle 
     1,224,637  stalled-cycles-backend # 0.16% backend cycles idle 

     0.203732797 seconds time elapsed 
------------------------------------------------------------ 


################## scalar version (cmdline arg) #############  
$ ocperf.py stat -e task-clock,cycles,instructions,uops_issued.any,uops_dispatched.thread,uops_retired.all,uops_retired.retire_slots,stalled-cycles-frontend,stalled-cycles-backend ./prefix-sum-avx n 

Performance counter stats for './prefix-sum-avx n': 

     287.567070  task-clock (msec)   # 0.999 CPUs utilized   
    1,082,611,453  cycles     # 3.765 GHz      
    2,381,840,355  instructions    # 2.20 insns per cycle   
                # 0.20 stalled cycles per insn 
    2,272,652,370  uops_issued_any   # 7903.034 M/sec     
    4,262,838,836  uops_dispatched_thread # 14823.807 M/sec     
    4,256,351,856  uops_retired_all   # 14801.249 M/sec     
    2,256,150,510  uops_retired_retire_slots # 7845.650 M/sec     
     465,018,146  stalled-cycles-frontend # 42.95% frontend cycles idle 
     6,321,098  stalled-cycles-backend # 0.58% backend cycles idle 

     0.287901811 seconds time elapsed 

------------------------------------------------------------  

shuffleが唯一のポート5、ないポート1上で実行することができますのでハスウェルは、ほぼ同じで、多分わずかに遅いごとのクロックなければならない（ベクトルのint addはまだハスウェル上/ 5 p1からさ。）

OTOH、IACAは、-funroll-loops（SnBを助ける）なしでコンパイルすると、1回の繰り返しでSnBより少し速くなると考えています。 Haswellはポート6でブランチを行うことができますが、SnBブランチではすでに飽和しているポート5にあります。

# compile without -DIACA_MARKS_OFF 
$ iaca -64 -mark 1 -arch HSW prefix-sum-avx  
Intel(R) Architecture Code Analyzer Version - 2.1 
Analyzed File - prefix-sum-avx 
Binary Format - 64Bit 
Architecture - HSW 
Analysis Type - Throughput 

******************************************************************* 
Intel(R) Architecture Code Analyzer Mark Number 1 
******************************************************************* 

Throughput Analysis Report 
-------------------------- 
Block Throughput: 6.20 Cycles  Throughput Bottleneck: Port5 

Port Binding In Cycles Per Iteration: 
--------------------------------------------------------------------------------------- 
| Port | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | 
--------------------------------------------------------------------------------------- 
| Cycles | 1.0 0.0 | 5.8 | 1.4 1.0 | 1.4 1.0 | 2.0 | 6.2 | 1.0 | 1.3 | 
--------------------------------------------------------------------------------------- 

N - port number or number of cycles resource conflict caused delay, DV - Divider pipe (on port 0) 
D - Data fetch pipe (on ports 2 and 3), CP - on a critical path 
F - Macro Fusion with the previous instruction occurred 
* - instruction micro-ops not bound to a port 
^ - Micro Fusion happened 
# - ESP Tracking sync uop was issued 
@ - SSE instruction followed an AVX256 instruction, dozens of cycles penalty is expected 
! - instruction not supported, was not accounted in Analysis 

| Num Of |     Ports pressure in cycles      | | 
| Uops | 0 - DV | 1 | 2 - D | 3 - D | 4 | 5 | 6 | 7 | | 
--------------------------------------------------------------------------------- 
| 1 |   |  | 1.0 1.0 |   |  |  |  |  | | vmovdqa xmm2, xmmword ptr [rax] 
| 1 | 1.0  |  |   |   |  |  |  |  | | add rax, 0x20 
| 1 |   |  |   | 1.0 1.0 |  |  |  |  | | vmovdqa xmm3, xmmword ptr [rax-0x10] 
| 1 |   |  |   |   |  | 1.0 |  |  | CP | vpslldq xmm1, xmm2, 0x4 
| 1 |   | 1.0 |   |   |  |  |  |  | | vpaddd xmm2, xmm2, xmm1 
| 1 |   |  |   |   |  | 1.0 |  |  | CP | vpslldq xmm1, xmm3, 0x4 
| 1 |   | 1.0 |   |   |  |  |  |  | | vpaddd xmm3, xmm3, xmm1 
| 1 |   |  |   |   |  | 1.0 |  |  | CP | vpslldq xmm1, xmm2, 0x8 
| 1 |   | 1.0 |   |   |  |  |  |  | | vpaddd xmm2, xmm2, xmm1 
| 1 |   |  |   |   |  | 1.0 |  |  | CP | vpslldq xmm1, xmm3, 0x8 
| 1 |   | 1.0 |   |   |  |  |  |  | | vpaddd xmm3, xmm3, xmm1 
| 1 |   | 0.9 |   |   |  | 0.2 |  |  | CP | vpaddd xmm1, xmm2, xmm0 
| 2^ |   |  |   |   | 1.0 |  |  | 1.0 | | vmovaps xmmword ptr [rax-0x20], xmm1 
| 1 |   |  |   |   |  | 1.0 |  |  | CP | vpshufd xmm1, xmm1, 0xff 
| 1 |   | 0.9 |   |   |  | 0.1 |  |  | CP | vpaddd xmm0, xmm1, xmm3 
| 2^ |   |  | 0.3  | 0.3  | 1.0 |  |  | 0.3 | | vmovaps xmmword ptr [rax-0x10], xmm0 
| 1 |   |  |   |   |  | 1.0 |  |  | CP | vpshufd xmm0, xmm0, 0xff 
| 1 |   |  |   |   |  |  | 1.0 |  | | cmp rax, 0x602020 
| 0F |   |  |   |   |  |  |  |  | | jnz 0xffffffffffffffa3 
Total Num Of Uops: 20 

ところで、gccが、私はループカウンタを持っていたとload(datavec + i + 1)をしていたとしても1レジスタのアドレッシングモードを使用するようにループをコンパイル。それが最高のコードです。 2レジスタアドレッシングモードがマイクロヒューズできないSnBファミリでは、clangのためにソースをそのループ条件に変更します。