somma cumulativa di prefisso (cumulativo) con SSE

Question

Sto cercando alcuni consigli su come eseguire una somma di prefisso parallelo con SSE. Sono interessato a farlo su una schiera di Ints, Floating o Double.

Ho trovato due soluzioni. Un caso speciale e un caso generale. In entrambi i casi la soluzione viene eseguita sull'array in due passaggi in parallelo con OpenMP. Per il caso speciale uso SSE su entrambe le passate. Per il caso generale lo uso solo al secondo passaggio.

La mia domanda principale è come utilizzare SSE al primo passaggio nel caso generale? Il seguente collegamento simd-prefix-sum-on-intel-cpu mostra un miglioramento per i byte ma non per i tipi di dati a 32 bit.

Il motivo per cui il caso speciale è chiamato speciale è che richiede che l'array sia in un formato speciale. Ad esempio, supponiamo che esistano solo 16 elementi di un array a di float. Poi, se la matrice fu rinviato simili (array di struct a struct di array):

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 somme verticali possono essere usati in entrambi i passaggi. Tuttavia, ciò sarebbe efficace solo se gli array fossero già nel formato speciale e l'output potesse essere utilizzato nel formato speciale. In caso contrario, sarebbe necessario effettuare un riorganizzazione costosa su input e output, il che renderebbe molto più lento rispetto al caso generale.

Forse dovrei considerare un algoritmo diverso per la somma del prefisso (ad esempio un albero binario)?

Codice per il caso generale:

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

Questa è la prima volta che sto rispondendo alla mia domanda, ma sembra opportuno. Basato su hirschhornsalz risposta per somma prefisso su 16 byte simd-prefix-sum-on-intel-cpu Ho trovato una soluzione per l'utilizzo di SIMD al primo passaggio per 4, 8 e 16 parole a 32 bit.

La teoria generale è la seguente. Per una scansione sequenziale delle parole n, sono necessarie le aggiunte n (n-1 per eseguire la scansione delle n parole e un'altra aggiunta effettuata dal gruppo di parole precedente scansionato). Tuttavia utilizzando SIMD n parole possono essere scansionate nel registro (n) aggiunte e un numero uguale di turni più un'aggiunta aggiuntiva e una trasmissione da trasportare dalla precedente ricerca SIMD. Quindi per un valore di n vincerà il metodo SIMD. sguardo

Let a parole a 32 bit con SSE, AVX, e 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 

sulla base di tale sembra SIMD non sarà utile per una scansione di parole a 32 bit fino AVX- 512. Ciò presuppone anche che gli spostamenti e le trasmissioni possano essere eseguiti solo in 1 istruzione. Questo è vero per SSE ma not for AVX and maybe not even for AVX2.

In ogni caso ho messo insieme un codice funzionante e testato che fa una somma di prefisso usando 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)); 
} 

Si noti che la funzione scan_SSE ha due aggiunte (_mm_add_ps) e due turni (_mm_slli_si128). I cast sono utilizzati solo per rendere felice il compilatore e non vengono convertiti in istruzioni. Quindi all'interno del ciclo principale sopra l'array in prefix_sum_SSE viene utilizzata un'altra aggiunta e una casuale. Si tratta di 6 operazioni totali rispetto alle sole 4 aggiunte con la somma sequenziale.

ecco una soluzione che lavora per 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); 
    } 
} 

I tre turni hanno bisogno di 7 intrinseche. La trasmissione ha bisogno di 2 elementi intrinseci. Quindi con le 4 aggiunte sono 13 gli intrinisti. Per AVX2 sono necessari solo 5 intrinsechi per i turni in modo da avere 11 intrinseci totali. La somma sequenziale ha bisogno solo di 8 aggiunte. Quindi probabilmente né AVX né AVX2 saranno utili per il primo passaggio.

Edit:

Così ho finalmente benchmark questo ei risultati sono inaspettati. Il codice di SSE e AVX sono entrambi circa due volte più veloce come il seguente codice sequenziale:

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

Credo che questo è dovuto al livello di istruzione paralellism.

In modo che risponda alla mia stessa domanda. Sono riuscito a utilizzare SIMD per pass1 nel caso generale. Quando combino questo con OpenMP sul mio sistema di bridge ivy a 4 core la velocità totale è di circa sette per 512k float.