c++ - Suma de prefijo SIMD en la CPU de Intel

sse mmx

Question

Necesito implementar un algoritmo de suma de prefijo y lo necesitaré para que sea lo más rápido posible. Ex:

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

¿Hay alguna forma de hacerlo utilizando la instrucción cpu SSE / mmx / SIMD?

¡Mi primera idea es sumar cada par en paralelo de forma recursiva hasta que toda la suma se haya calculado como se muestra a continuación!

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

Para que el algoritmo sea un poco más claro, "z" no es la salida final

pero en su lugar se utiliza para calcular la salida

 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

El algoritmo de suma de prefijos paralelos más rápido que conozco es sobrepasar la suma en dos pasadas en paralelo y usar SSE también en la segunda pasada.

En la primera pasada, usted calcula sumas parciales en paralelo y almacena la suma total para cada suma parcial. En el segundo paso, agregue la suma total de la suma parcial anterior a la siguiente suma parcial. Puede ejecutar ambas pasadas en paralelo utilizando varios subprocesos (por ejemplo, con OpenMP). La segunda pasada también puede usar SIMD ya que se agrega un valor constante a cada suma parcial.

Suponiendo que n elementos de una matriz, m núcleos y un ancho SIMD de w el costo de tiempo debe ser

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

Como la primera pasada no usa SIMD, el costo del tiempo siempre será mayor que n/m

Por ejemplo, para cuatro núcleos con un SIMD_width de 4 (cuatro flotantes de 32 bits con SSE) el costo sería de 5 n 5n/16 . O aproximadamente 3.2 veces más rápido que el código secuencial que tiene un costo de tiempo de n . Usando Hyper Threading la velocidad será aún mayor.

En casos especiales, también es posible usar SIMD en la primera pasada. Entonces el costo del tiempo es simplemente

2*n/(m*w)

Publiqué el código para el caso general que usa OpenMP para los subprocesos e intrínsecos para el código SSE y discuto los detalles sobre el caso especial en el siguiente enlace parallel-prefix-cumulative-sum-with-sse

Edición: logré encontrar una versión SIMD para la primera pasada, que es aproximadamente el doble de rápida que el código secuencial. Ahora recibo un aumento total de aproximadamente 7 en mi sistema de cuatro núcleos de puente ivy.

Edición: para matrices más grandes, un problema es que, después del primer paso, la mayoría de los valores se han desalojado de la memoria caché. Se me ocurrió una solución que se ejecuta en paralelo dentro de un trozo, pero ejecuta cada trozo en serie. El chunk_size es un valor que debe ajustarse. Por ejemplo lo puse a 1MB = 256K flotantes. Ahora el segundo paso se realiza mientras los valores aún están dentro del caché de nivel 2. Hacer esto da una gran mejora para los arreglos grandes.

Aquí está el código para SSE. El código AVX tiene aproximadamente la misma velocidad, así que no lo publiqué aquí. La función que hace la suma del prefijo es scan_omp_SSEp2_SSEp1_chunk . Pásale una matriz de flotadores y llena la matriz con la suma acumulada.

__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

Para una matriz de 1000 enteros de 32 bits, pude obtener una pequeña aceleración de aproximadamente 1.4x de un solo hilo, usando el método de @ hirschhornsalz en un bucle en Intel Sandybridge. Con un búfer de 60 kB de ints, la aceleración es de aproximadamente 1.37. Con 8MiB de ints, la aceleración sigue siendo 1.13. (i5-2500k a 3.8GHz turbo, con DDR3-1600.)

Los elementos más pequeños ( int16_t o uint8_t , o las versiones sin firma) llevarán una etapa adicional de cambio / adición para cada duplicación del número de elementos por vector. El desbordamiento es malo, así que no intente utilizar un tipo de datos que no pueda contener la suma de todos los elementos, a pesar de que le da a SSE una mayor ventaja.

#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;
}

Prueba con n = 1000, con la lista compilada en el binario. (Y sí, verifiqué que en realidad es un bucle, sin tomar ningún atajo en tiempo de compilación que haga que la prueba vectorial o no vectorial no tenga sentido).

Tenga en cuenta que la compilación con AVX para obtener instrucciones vectoriales no destructivas de 3 operandos guarda una gran cantidad de instrucciones movdqa , pero solo guarda una pequeña cantidad de ciclos. Esto se debe a que shuffle y vector-int-add solo pueden ejecutarse en los puertos 1 y 5, en SnB / IvB, por lo que port0 tiene muchos ciclos de repuesto para ejecutar las instrucciones de movimiento. Los cuellos de botella en el rendimiento de uop-cache pueden ser la razón por la que la versión no AVX es ligeramente más lenta. (Todas esas instrucciones de movimiento adicional nos empujan hasta 3.35 insn / ciclo). El frontend solo está inactivo al 4,54% de los ciclos, por lo que apenas se mantiene.

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 debería ser más o menos igual, pero quizás un poco más lento por hora, ya que la reproducción aleatoria solo puede ejecutarse en el puerto 5, no en el puerto 1. (vector-int add aún es p1 / 5 en Haswell).

OTOH, IACA cree que Haswell será un poco más rápido que SnB para una iteración, si compilas sin -funroll-loops (lo que ayuda en SnB). Haswell puede hacer sucursales en el puerto 6, pero en SnB las sucursales están en el puerto 5, que ya saturamos.

 # 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

Por cierto, gcc compiló el bucle para utilizar un modo de direccionamiento de registro único incluso cuando tenía un contador de bucle y estaba haciendo load(datavec + i + 1) . Ese es el mejor código, esp. en la familia SnB, donde los modos de direccionamiento de 2 registros no se pueden fusionar, así que cambio la fuente a esa condición de bucle para el beneficio del ruido.

Answer 3

Puede explotar un pequeño paralelismo para grandes longitudes de registro y pequeñas sumas. Por ejemplo, sumar 16 valores de 1 byte (que se ajustan a un registro sse) solo requiere el registro de ₂ 16 adiciones y un número igual de turnos.
No mucho, pero más rápido que 15 adiciones dependientes y los accesos de memoria adicionales.

__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 tiene sumas más largas, las dependencias podrían ocultarse explotando el paralelismo de nivel de instrucción y aprovechando el reordenamiento de la instrucción.

Edición: algo como

__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

prefix-sum se puede calcular en paralelo, en realidad es uno de los algoritmos fundamentales en la programación de GPU. Si está utilizando extensiones SIMD en un procesador Intel, no estoy seguro de que hacerlo en paralelo realmente le beneficie mucho, pero eche un vistazo a este documento de nvidia sobre la implementación de la suma de prefijos paralelos (solo mire los algoritmos e ignórelos). la CUDA): http://http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html