sum paralela del prefijo (acumulativa) con SSE

Estoy buscando algunos consejos sobre cómo hacer una sum de prefijo paralela con SSE. Estoy interesado en hacer esto en una variedad de ints, floats o dobles.

He encontrado dos soluciones. Un caso especial y un caso general. En ambos casos, la solución se ejecuta en la matriz en dos pasos en paralelo con OpenMP. Para el caso especial, uso SSE en ambos pases. Para el caso general, lo uso solo en el segundo pase.

Mi pregunta principal es cómo puedo usar SSE en el primer pase en el caso general. El siguiente enlace simd-prefix-sum-on-intel-cpu muestra una mejora para los bytes pero no para los tipos de datos de 32 bits.

La razón por la cual el caso especial se llama especial es que requiere que la matriz esté en un formato especial. Por ejemplo, supongamos que solo hay 16 elementos de una matriz de flotantes. Entonces, si la matriz se reorganizó de esta manera (array de estructuras a struct de matrices):

 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] 

Sumas verticales de SSE podrían usarse en ambas pasadas. Sin embargo, esto solo sería eficiente si las matrices ya estuvieran en el formato especial y la salida se pudiera usar en el formato especial. De lo contrario, la reordenación costosa tendría que hacerse tanto en la entrada como en la salida, lo que lo haría mucho más lento que el caso general.

Tal vez debería considerar un algoritmo diferente para la sum del prefijo (por ejemplo, un árbol binario)?

Código para el caso general:

 void prefix_sum_omp_sse(double a[], double s[], int n) { double *sum; #pragma omp parallel { const int ithread = omp_get_thread_num(); const int nthreads = omp_get_num_threads(); #pragma omp single { sum = new double[nthreads + 1]; sum[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; } 

Esta es la primera vez que respondo mi propia pregunta, pero parece apropiada. Basado en la respuesta de hirschhornsalz para prefix sum en 16 bytes simd-prefix-sum-on-intel-cpu He encontrado una solución para usar SIMD en la primera pasada para 4, 8 y 16 palabras de 32 bits.

La teoría general es la siguiente. Para un escaneo secuencial de n palabras, se necesitan n adiciones (n-1 para escanear las n palabras y una adición más del conjunto anterior de palabras escaneadas). Sin embargo, el uso de SIMD n palabras se puede escanear en log 2 (n) adiciones y un número igual de turnos más una adición más y difusión para llevar desde el escaneo SIMD anterior. Entonces, para algún valor de n el método SIMD ganará.

Miremos las palabras de 32 bits con SSE, AVX y 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 

Basado en eso, parece que SIMD no será útil para un escaneo de palabras de 32 bits hasta el AVX-512. Esto también supone que los cambios y la transmisión se pueden realizar en solo 1 instrucción. Esto es cierto para SSE pero no para AVX y tal vez ni siquiera para AVX2 .

En cualquier caso, armé un código de trabajo y probado que hace una sum de prefijos 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)); } 

Observe que la función scan_SSE tiene dos adiciones (_mm_add_ps) y dos cambios (_mm_slli_si128). Los moldes solo se usan para hacer feliz al comstackdor y no se convierten en instrucciones. Luego dentro del bucle principal sobre la matriz en prefix_sum_SSE otra adición y una mezcla. Eso es 6 operaciones en total en comparación con solo 4 adiciones con la sum secuencial.

Aquí hay una solución de trabajo para 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); } } 

Los tres turnos necesitan 7 intrínsecos. La transmisión necesita 2 intrínsecos. Entonces con las 4 adiciones eso es 13 intrínsecas. Para AVX2, solo se necesitan 5 intrínsecos para los turnos, por lo que se necesitan 11 intrínsecos en total. La sum secuencial solo necesita 8 adiciones. Por lo tanto, es probable que ni AVX ni AVX2 sean útiles para el primer pase.

Editar:

Así que finalmente hice una evaluación comparativa de esto y los resultados son inesperados. El código SSE y AVX son dos veces más rápidos que el siguiente código secuencial:

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

Supongo que esto se debe al paralelismo de nivel de instrucción.

Entonces eso responde mi propia pregunta. Logré usar SIMD para pass1 en el caso general. Cuando combino esto con OpenMP en mi sistema de puente de hiedra de 4 núcleos, la velocidad total es de aproximadamente siete para flotadores de 512k.