suma de prefijo SIMD en la CPU de Intel

Question

Necesito implementar un algoritmo de suma de prefijos y necesitaría 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 una manera de hacer esto utilizando la instrucción de la CPU SSE/MMX/SIMD?

Mi primera idea es sumar cada par en paralelo recursivamente hasta que toda la suma se haya calculado como 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 un poco más claro "z" no es que la salida definitiva

sino que se utiliza para calcular que 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 prefijo paralelo más rápido que conozco es de atropellar a la suma en dos pasadas en paralelo y utilizar SSE así en la segunda pasada.

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

Suponiendo n elementos de una matriz, m núcleos, y una anchura SIMD de w el coste de tiempo debe ser

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

Desde el pase puño no utiliza SIMD el costo de tiempo siempre será mayor que n/m

Por ejemplo, para cuatro núcleos con SIMD_width de 4 (cuatro flotadores de 32 bits con SSE) el costo sería 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, es posible usar SIMD en la primera pasada también.Entonces el costo de tiempo es simplemente

2*n/(m*w) 

he publicado el código para el caso general que utiliza OpenMP para el roscado y los intrínsecos para el código de SSE y discutir detalles sobre el caso especial en el siguiente enlace parallel-prefix-cumulative-sum-with-sse

Editar: Logré encontrar una versión SIMD para el primer pase, que es aproximadamente el doble de rápido que el código secuencial. Ahora recibo un impulso total de aproximadamente 7 en mi sistema de puente de hiedra de cuatro núcleos.

Editar: Para más grandes arrays un problema es que después de la primera pasada de 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 fragmento, pero ejecuta cada fragmento en serie. El chunk_size es un valor que debe ajustarse. Por ejemplo, lo configuré en 1MB = 256K flotantes. Ahora el segundo pase se realiza mientras los valores están todavía dentro de la memoria caché de nivel 2. Hacer esto ofrece una gran mejora para las matrices grandes.

Aquí está el código para SSE. El código AVX es 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áselo una matriz de a de flotadores y rellena la matriz s con la suma acumulativa.

__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

prefijo de suma puede ser calculada de forma paralela , 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 si hacerlo en paralelo realmente lo beneficiará mucho, pero eche un vistazo a este documento de nvidia sobre la implementación de prefix-sum en paralelo (solo mire los algoritmos e ignore) CUDA): http://http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html

Answer 3

Puede aprovechar un pequeño paralelismo para longitudes grandes de registro y pequeñas sumas. Por ejemplo, la suma de 16 valores de 1 byte (que casualmente encajan en un registro sse) requiere solo log 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 largos, las dependencias puede estar escondido tras explotar el paralelismo a nivel de instrucción y el aprovechamiento de reordenamiento de instrucciones.

Editar: algo así 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

Para un arsenal de 1000 enteros de 32 bits, yo era capaz de conseguir un pequeño aumento de velocidad de aproximadamente 1,4x método de un solo subproceso, de usar @ Hirschhornsalz en un bucle en Intel SandyBridge. Con un buffer de entrada 60kiB, la aceleración es de aproximadamente 1.37. Con 8MiB de ints, la aceleración sigue siendo 1.13. (I5-2500K en turbo 3,8 GHz, con DDR3-1600.)

más pequeño elementos (int16_t o uint8_t, o las versiones sin firmar) tomaría una etapa extra de cambio/añadir para cada duplicación del número de elementos por vector . El desbordamiento es malo, por lo tanto, 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; 
} 

---

Testing con n = 1.000, con la lista compilada en el binario. (Y sí, comprobé que en realidad se está ejecutando un bucle, sin tomar ningún atajo en tiempo de compilación que haga que el vector o la prueba no vectorial carezcan de sentido).

Tenga en cuenta que compilar con AVX para obtener instrucciones vectoriales no destructivas de 3 operandos ahorra muchas instrucciones movdqa, pero solo se ahorra una pequeña cantidad de ciclos. Esto se debe a que shuffle y vector-int-add solo se pueden ejecutar en los puertos 1 y 5, en SnB/IvB, por lo que port0 tiene muchos ciclos de repuesto para ejecutar las instrucciones mov. Los cuellos de botella de rendimiento de la caché de uops pueden ser la razón por la cual la versión que no es AVX es un poco más lenta. (Todas esas instrucciones extra mov nos empujan hasta 3.35 insn/ciclo). El frontend solo está inactivo al 4.54% de los ciclos, por lo que apenas se mantiene al día.

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 aproximadamente el mismo, pero tal vez un poco más lento por ciclo de reloj, porque aleatoria sólo se puede ejecutar en el puerto 5, no el puerto 1. (vector-int añadir es todavía P1/5 sobre Haswell.)

OTOH, IACA cree que Haswell será ligeramente más rápido que SnB por una iteración, si compila sin -funroll-loops (lo que sí ayuda en SnB). Haswell puede hacer ramas en port6, pero en SnB las ramas están en port5, 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 compila el bucle de utilizar un modo de direccionamiento de un registro aun cuando no 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 pueden micro fusibles, entonces cambio la fuente a esa condición de ciclo para el beneficio del clang.