Somme du préfixe SIMD sur processeur Intel

Question

J'ai besoin d'implémenter un algorithme de somme de préfixe et je voudrais qu'il soit le plus rapide possible. Ex:

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

Est-il possible de le faire en utilisant l'instruction cpu SSE/MMX/SIMD?

Ma première idée est de sommer chaque paire en parallèle de façon récursive jusqu'à ce que toute la somme ait été calculée comme ci-dessous!

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

Pour l'algorithme un peu plus clair « z » est pas la ouput finale

mais utilisé pour calculer le 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

Le plus rapide algorithme de somme préfixe parallèle que je connaisse est de courir sur la somme en deux passes en parallèle et utiliser SSE aussi bien dans le deuxième passage.

Lors de la première passe, vous calculez des sommes partielles en parallèle et stockez la somme totale pour chaque somme partielle. Dans la deuxième passe, vous ajoutez la somme totale de la somme partielle précédente à la somme partielle suivante. Vous pouvez exécuter les deux passes en parallèle en utilisant plusieurs threads (par exemple avec OpenMP). La seconde passe, vous pouvez également utiliser SIMD car une valeur constante est ajoutée à chaque somme partielle.

En supposant n éléments d'un tableau, m noyaux, et une largeur de SIMD de w le coût du temps devrait être

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

Depuis le col de poing n'utilise pas SIMD le coût du temps sera toujours supérieure à n/m

Par exemple pour quatre cœurs avec une largeur SIMD_4 (quatre flotteurs 32 bits avec SSE) le coût serait 5n/16. Ou environ 3,2 fois plus rapide que le code séquentiel qui a un coût en temps de n. En utilisant l'hyper-threading, l'accélération sera encore plus grande.

Dans certains cas, il est également possible d'utiliser SIMD lors du premier passage.Ensuite, le coût du temps est tout simplement

2*n/(m*w) 

J'ai posté le code pour le cas général qui utilise OpenMP pour le filetage et intrinsics pour le code SSE et discuter des détails sur le cas particulier sur le lien suivant parallel-prefix-cumulative-sum-with-sse

Editer: J'ai réussi à trouver une version SIMD pour le premier passage qui est environ deux fois plus rapide que le code séquentiel. Maintenant, je reçois un total d'environ 7 sur mon système de quatre noyaux de pont de lierre.

Edit: Pour les tableaux plus grand problème est que, après la première passe la plupart des valeurs ont été expulsées du cache. Je suis venu avec une solution qui fonctionne en parallèle à l'intérieur d'un morceau, mais exécute chaque morceau en série. Le chunk_size est une valeur qui doit être ajustée. Par exemple je l'ai mis à 1MB = 256K flotteurs. Maintenant, la deuxième passe est effectuée alors que les valeurs sont toujours dans le cache de niveau 2. Faire cela donne une grande amélioration pour les grandes baies.

Voici le code pour SSE. Le code AVX est à peu près la même vitesse donc je ne l'ai pas posté ici. La fonction qui fait la somme du préfixe est scan_omp_SSEp2_SSEp1_chunk. Passez-lui un tableau a de flottants et il remplit le tableau s avec la somme cumulée.

__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

préfixe somme

peut être calculée en parallèle , c'est en fait l'un des algorithmes fondamentaux de la programmation GPU. Si vous utilisez des extensions SIMD sur un processeur Intel, je ne suis pas sûr que le faire en parallèle vous sera très utile, mais jetez un oeil à ce document de nvidia sur l'implémentation de parallel-sum-sum (regardez les algorithmes et ignorez le CUDA): http://http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html

Answer 3

Vous pouvez exploiter un parallélisme mineur pour les grandes longueurs de registre et les petites sommes. Par exemple, l'addition de 16 valeurs de 1 octet (qui se rangent dans un registre sse) ne nécessite que l'ajout de 16 additions et d'un nombre égal de décalages.
Pas beaucoup, mais plus rapide que 15 ajouts dépendants et les accès mémoire supplémentaires.

__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 

Si vous avez des sommes plus, les dépendances pourraient être cachés en exploitant le parallélisme au niveau des instructions et en tirant parti de réordonnancement d'instruction.

Edit: quelque chose comme

__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

Pour un tableau de 1000 entiers de 32 bits, j'ai pu obtenir un petit gain de vitesse d'environ mono-thread 1.4x, méthode de l'aide @ hirschhornsalz dans une boucle sur Intel Sandybridge. Avec un tampon de 60 Ko d'ints, l'accélération est d'environ 1,37. Avec 8MiB d'ints, l'accélération est toujours de 1.13. (I5-2500K à turbo 3.8GHz, avec DDR3-1600.)

éléments plus petits (ou int16_tuint8_t, ou les versions non signés) prendrait une étape supplémentaire de décalage/ajouter pour chaque doublement du nombre d'éléments par vecteur . Le débordement est mauvais, donc n'essayez pas d'utiliser un type de données qui ne peut pas contenir la somme de tous les éléments, même si cela donne un plus grand avantage à 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; 
} 

---

Test avec n = 1000, avec la liste établie dans le binaire. (Et oui, j'ai vérifié qu'il est en boucle, ne prenant aucun raccourci de compilation qui rend le test vectoriel ou non-vectoriel inutile.)

Notez que la compilation avec AVX pour obtenir des instructions vectorielles non-destructives à 3 opérandes beaucoup d'instructions movdqa, mais enregistre seulement une petite quantité de cycles. En effet, shuffle et vector-int-add ne peuvent s'exécuter que sur les ports 1 et 5, sur SnB/IvB, donc port0 dispose de beaucoup de cycles de rechange pour exécuter les instructions mov. Les goulets d'étranglement du débit de uop-cache pourraient être la raison pour laquelle la version non-AVX est légèrement plus lente. (Toutes ces instructions supplémentaires nous poussent jusqu'à 3.35 insn/cycle). Le frontend n'est que 4.54% de cycles inutilisés, donc il se maintient à peine.

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 

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

Haswell doit être sur le même, mais peut-être un peu plus lent par horloge, parce que lecture aléatoire ne peut fonctionner que sur le port 5, pas le port 1. (ajouter un vecteur int est encore p1/5 sur Haswell.)

OTOH, IACA pense que Haswell sera légèrement plus rapide que SnB pour une itération, si vous compilez sans -funroll-loops (ce qui aide sur SnB). Haswell peut faire des branches sur port6, mais sur les branches SnB sont sur port5, que nous saturons déjà.

# 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 

BTW, gcc compilé la boucle d'utiliser un mode d'adressage d'un registre, même quand j'avais un compteur de boucle et faisais load(datavec + i + 1). C'est le meilleur code, esp. sur la famille SnB où les modes d'adressage à 2 registres ne peuvent pas fusionner, donc je change la source à cette condition de boucle pour le bénéfice du clang.