¿Es posible vectorizar myNum + = a [b [i]] * c [i]; en x86_64?

Question

¿Qué intrínsecos usaría para vectorizar lo siguiente (si es que es posible vectorizar) en el x86_64?

double myNum = 0; 
for(int i=0;i<n;i++){ 
    myNum += a[b[i]] * c[i]; //b[i] = int, a[b[i]] = double, c[i] = double 
} 

Answer 1

Aquí es mi ir en él, completamente optimizado y probado:

#include <emmintrin.h> 

__m128d sum = _mm_setzero_pd(); 
for(int i=0; i<n; i+=2) { 
    sum = _mm_add_pd(sum, _mm_mul_pd(
     _mm_loadu_pd(c + i), 
     _mm_setr_pd(a[b[i]], a[b[i+1]]) 
    )); 
} 

if(n & 1) { 
    sum = _mm_add_pd(sum, _mm_set_sd(a[b[n-1]] * c[n-1])); 
} 

double finalSum = _mm_cvtsd_f64(_mm_add_pd(
    sum, _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1)) 
)); 

Esto produce el código de montaje muy bonito usando gcc -O2 -msse2 (4.4.1).

Como se puede decir, que tiene una aún n hará que este circuito sea más rápido, así como un c alineados. Si puede alinear c, cambiar _mm_loadu_pd-_mm_load_pd durante los tiempos de ejecución aún más rápidas.

Answer 2

respuesta corta no. Larga respuesta, sí, pero no de manera eficiente. Usted incurrirá en la multa por realizar cargas no alineadas que anularán cualquier tipo de beneficio. A menos que pueda garantizar que b [i] índices sucesivos estén alineados, lo más probable es que tenga un peor rendimiento después de la vectorización

si sabe de antemano qué índices son, lo mejor es desenrollar y especificar índices explícitos. Hice algo similar usando la especialización de plantilla y la generación de código. si te interesa, puedo compartir

para responder a su comentario, que básicamente tienen que concentrarse en una matriz. Lo más fácil para tratar de inmediato es el de bloquear el bucle por un factor de dos, la carga baja y alta una por separado, y luego usar mm * _pd como normalmente. Pseudocódigo:

__m128d a, result; 
for(i = 0; i < n; i +=2) { 
    ((double*)(&a))[0] = A[B[i]]; 
    ((double*)(&a))[1] = A[B[i+1]]; 
    // you may also load B using packed integer instruction 
    result = _mm_add_pd(result, _mm_mul_pd(a, (__m128d)(C[i]))); 
} 

No recuerdo exactamente los nombres de las funciones, es posible que desee verificarlos dos veces. Además, utilice la palabra clave restrict con los punteros si sabe que no puede haber problemas de aliasing. Esto permitirá que el compilador sea mucho más agresivo.

Answer 3

esto no va a vectorizar como es, debido a la doble vía indirecta de los índices de matriz. Como trabajas con dobles, no se gana mucho o nada de SSE, especialmente porque la mayoría de las CPU modernas tienen 2 FPU de todos modos.

Answer 4

Me gustaría empezar por desenrollar el bucle. Algo así como

double myNum1 = 0, myNum2=0; 
for(int i=0;i<n;i+=2) 
{ 
    myNum1 += a[b[ i ]] * c[ i ]; 
    myNum2 += a[b[i+1]] * c[i+1]; 
} 
// ...extra code to handle the remainder when n isn't a multiple of 2... 
double myNum = myNum1 + myNum2; 

Esperemos que permita al compilador intercalar las cargas con la aritmética; perfil y mire la asamblea para ver si hay una mejora. Idealmente, el compilador generará instrucciones de SSE, pero no estoy si eso sucede en la práctica.

Desenrollar más podría dejar de hacer esto:

__m128d sum0, sum1; 
// ...initialize to zero... 
for(int i=0;i<n;i+=4) 
{ 
    double temp0 = a[b[ i ]] * c[ i ]; 
    double temp1 = a[b[i+1]] * c[i+1]; 
    double temp2 = a[b[i+2]] * c[i+2]; 
    double temp3 = a[b[i+3]] * c[i+3]; 
    __m128d pair0 = _mm_set_pd(temp0, temp1); 
    __m128d pair1 = _mm_set_pd(temp2, temp3); 
    sum0 = _mm_add_pd(sum0, pair0); 
    sum1 = _mm_add_pd(sum1, pair1); 
} 
// ...extra code to handle the remainder when n isn't a multiple of 4... 
// ...add sum0 and sum1, then add the result's components... 

(disculpas por el pseudocódigo al comienzo y al final, supongo que la parte importante era el bucle). No estoy seguro si eso será más rápido; depende de las diversas latencias y de qué tan bien el compilador puede reorganizar todo. Asegúrese de que su perfil antes y después para ver si hubo una mejora real.

Espero que ayude.

Answer 5

procesadores Intel puede emitir dos operaciones de punto flotante, pero una carga por ciclo, por lo que los accesos de memoria son la restricción más apretado.Con esto en mente, primero apunté a usar cargas empaquetadas para reducir el número de instrucciones de carga, y usé la aritmética compacta solo porque era conveniente. Desde entonces, me he dado cuenta de que saturar el ancho de banda de la memoria puede ser el mayor problema, y ​​todo el juego con las instrucciones SSE podría haber sido una optimización prematura si el objetivo era hacer que el código fuera rápido en lugar de aprender a vectorizar.

SSE

Las cargas menor número posible con ninguna suposición sobre los índices en b requiere desenrollar el bucle de cuatro veces. Una carga de 128 bits obtiene cuatro índices de b, dos cargas de 128 bits cada una obtiene un par de dobles adyacentes desde c, y la recopilación de a requiere cargas independientes de 64 bits. Eso es un piso de 7 ciclos por cada cuatro iteraciones para el código de serie. (suficiente para saturar el ancho de banda de mi memoria si el acceso a a no guarda en caché). He dejado de lado algunas cosas molestas como manejar un número de iteraciones que no es un múltiplo de 4.

entry: ; (rdi,rsi,rdx,rcx) are (n,a,b,c) 
    xorpd xmm0, xmm0 
    xor r8, r8 
loop: 
    movdqa xmm1, [rdx+4*r8] 
    movapd xmm2, [rcx+8*r8] 
    movapd xmm3, [rcx+8*r8+8] 
    movd r9, xmm1 
    movq r10, xmm1 
    movsd xmm4, [rsi+8*r9] 
    shr r10, 32 
    movhpd xmm4, [rsi+8*r10] 
    punpckhqdq xmm1, xmm1 
    movd r9, xmm1 
    movq r10, xmm1 
    movsd xmm5, [rsi+8*r9] 
    shr r10, 32 
    movhpd xmm5, [rsi+8*r10] 
    add r8, 4 
    cmp r8, rdi 
    mulpd xmm2, xmm4 
    mulpd xmm3, xmm5 
    addpd xmm0, xmm2 
    addpd xmm0, xmm3 
    jl loop 

Obtención de los índices a cabo es la parte más complicada. movdqa carga 128 bits de datos enteros desde una dirección alineada de 16 bytes (Nehalem tiene penalizaciones de latencia para mezclar las instrucciones de SSE "entero" y "flotante"). punpckhqdq mueve alto 64 bits a 64 bits bajo, pero en modo entero a diferencia del nombre más simple movhlpd. Los cambios de 32 bits se realizan en los registros de propósito general. movhpd carga una doble en la parte superior de un registro xmm sin molestar a la parte inferior; esto se usa para cargar los elementos de a directamente en los registros empaquetados.

Este código es claramente más rápido que el código anterior, que a su vez es más rápido que el código simple, y en cada patrón de acceso, excepto en el caso simple B[i] = i donde el lazo ingenuo es realmente el más rápido. También intenté algunas cosas como una función alrededor de SUM(A(B(:)),C(:)) en Fortran que terminó básicamente equivalente al bucle simple.

Probé en un Q6600 (65 nm Core 2 a 2.4Ghz) con 4 GB de memoria DDR2-667, en 4 módulos. Probar el ancho de banda de la memoria da unos 5333 MB/s, por lo que parece que solo veo un solo canal. Estoy compilando con gcc 4.3.2-1.1, -O3 -Fast-math -msse2 -Ftree-vectorize -std = gnu99 de Debian.

Para la prueba de que estoy dejando n ser un millón, inicializar las matrices de manera a[b[i]] y c[i] ambos iguales 1.0/(i+1), con unos patrones diferentes de índices. Una asigna a con un millón de elementos y conjuntos de b a una permutación aleatoria, otro asigna a con elementos 10M y utiliza cada décima, y ​​los últimos asigna a con elementos 10M y configura b[i+1] mediante la adición de un número al azar de 1 a 9 para b[i]. Estoy cronometrando cuánto tarda una llamada con gettimeofday, borrando las cachés llamando al clflush en las matrices y midiendo 1000 pruebas de cada función. Tracé distribuciones de tiempo de ejecución suavizadas usando algún código de las entrañas de criterion (en particular, el estimador de densidad de kernel en el paquete statistics).

ancho de banda

Ahora, para la nota real importante acerca de ancho de banda. 5333MB/s con reloj de 2.4Ghz es un poco más de dos bytes por ciclo. Mis datos son lo suficientemente largos como para que nada pueda almacenarse en caché, y multiplicar el tiempo de ejecución de mi ciclo por (16 + 2 * 16 + 4 * 64) bytes cargados por iteración si todo falla me da casi exactamente el ancho de banda de ~ 5333MB/s que mi sistema tiene . Debería ser bastante fácil saturar ese ancho de banda sin SSE.Aun suponiendo a se almacenan en caché por completo, sólo la lectura y bc para una iteración mueve 12 bytes de datos, y la ingenua puede iniciar una nueva iteración de cada tercer ciclo con la canalización.

Suponiendo que algo menos que el almacenamiento en caché completo en a hace que la aritmética y las instrucciones cuenten incluso menos de un cuello de botella. No me sorprendería que la mayor parte del aumento de velocidad en mi código proviene de la emisión de un menor número de cargas a bc y espacio para que más está libre para realizar un seguimiento y especular últimos fallos de caché en a.

de hardware más amplia podría tener más diferencia. Un sistema Nehalem con tres canales de DDR3-1333 necesitaría mover 3 * 10667/2.66 = 12.6 bytes por ciclo para saturar el ancho de banda de la memoria. Eso sería imposible para un único hilo si a encaja en el caché, pero a 64 bytes carece de un caché de línea en el vector se agrega rápidamente, solo una de las cuatro cargas en mi bucle falta en cachés trae el ancho de banda promedio requerido a 16 bytes /ciclo.