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

Question

¿Qué características intrínsecas debo utilizar para vectorizar el siguiente(si es posible incluso para 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

Esta 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 hermosa usando gcc -O2 -msse2 (4.4.1).

Como se puede decir, que tienen un n incluso hará que este circuito sea más rápido, así como un c alineados. Si puede alinear c, el cambio _mm_loadu_pd a _mm_load_pd para una aún más rápido los tiempos de ejecución.

Answer 2

Yo empezaría 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;

Con suerte que le permite al compilador para intercalar las cargas con la aritmética; el perfil y la mirada a la asamblea para ver si hay una mejora. Lo ideal sería que el compilador generará instrucciones SSE, pero no soy si eso sucede en la práctica.

Unrolling 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 en el inicio y el final, I figura la parte importante fue el bucle). No sé a ciencia cierta si será más rápido; que depende de las distintas latencias y lo bien que el compilador puede reorganizar todo. Asegúrese de perfilar antes y después para ver si hubo una mejora real.

Espero que ayude.

Answer 3

Los procesadores Intel puede emitir dos operaciones de punto flotante, pero una carga por ciclo, así que los accesos a la memoria son los más estrechos restricción.Con eso en mente, me dirigido por primera vez a utilizar lleno de cargas para reducir el número de instrucciones de carga, y se utiliza lleno aritmética sólo porque era conveniente.Desde entonces me he dado cuenta de que saturando el ancho de banda de memoria puede ser el mayor problema, y todo el cachondeo con las instrucciones SSE podría haber sido prematuro de la optimización, si el punto era hacer que el código de ir rápido en lugar de aprender a vectorizar.

ESS

El menor número posible de cargas con ninguna suposición sobre los índices en b requiere de desenrollar el bucle de cuatro veces.Uno de 128 bits de carga consigue cuatro índices de b, dos de 128 bits cargas reciben cada uno un par de dobles adyacentes de c, y la recopilación de a se requiere independiente de 64 bits de carga.Que un piso de 7 ciclos por cuatro iteraciones para el código de serie.(lo suficiente como para saturar mi ancho de banda de memoria si el acceso a a no cache muy bien).Me he dejado algunas cosas molestas como el manejo de 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 es la parte más complicada. movdqa las cargas de 128 bits de datos enteros de 16 bytes alineado de dirección (Nehalem ha latencia de sanciones para la mezcla de la "entero" y "flotan" las instrucciones SSE). punpckhqdq movimientos de alta de 64 bits a la baja de 64 bits, pero en números enteros modo, a diferencia de la manera más simple nombre movhlpd.32 bits de los turnos de hecho en la general de registros de propósito. movhpd carga un doble en la parte superior de un xmm registro sin molestar a los de la parte inferior - este es usado para cargar los elementos de a directamente en el empacado de los registros.

Este código claramente más rápido que el código de arriba, la cual es más rápido que el código simple, y en cada acceso motivo, pero el caso simple B[i] = i donde el ingenuo bucle es en realidad más rápido.También probé con un par de cosa como una función alrededor de SUM(A(B(:)),C(:)) en Fortran que terminó básicamente equivalente a la de bucle simple.

He probado en un Q6600 (65 nm Core 2 a 2,4 Ghz) con 4 gb de memoria DDR2-667, en 4 módulos.Prueba de ancho de banda de memoria da sobre 5333 MB/s, así que parece que sólo estoy viendo un solo canal.Estoy compilando con Debian gcc 4.3.2-1.1, -O3 -Ffast-matemáticas -msse2 -Ftree-vectorización -std=gnu99.

Para las pruebas estoy dejando n ser uno de los millones, la inicialización de los arrays para a[b[i]] y c[i] los dos por igual 1.0/(i+1), con un par de diferentes patrones de índices.Uno asigna a con un millón de elementos y conjuntos de b al azar a una permutación, otro asigna a con 10 millones de elementos y usos de cada 10, y la última asigna a con 10 millones de elementos y conjuntos de hasta b[i+1] mediante la adición de un número aleatorio de 1 a 9 a b[i].Estoy de temporización de tiempo durante el que se toma con gettimeofday, borrar las cachés llamando clflush sobre las matrices, y la medición de 1000 ensayos de cada función.He trazado alisado de tiempo de ejecución de las distribuciones mediante el código de las entrañas de criterio (en particular, el estimador de densidad de kernel en la statistics el paquete).

El ancho de banda

Ahora, para la nota importante acerca de ancho de banda.5333MB/s, 2.4 Ghz de reloj es de poco más de dos bytes por ciclo.Mis datos es lo suficientemente larga que nada debe ser almacenable, y multiplicando el tiempo de ejecución de mi bucle (16+2*16+4*64) bytes cargados por iteración si todo se pierde me da casi exactamente el ~5333MB/s de ancho de banda tiene mi sistema.Debería ser bastante fácil para que saturar el ancho de banda sin ESS.Aun suponiendo que la a estaban completamente en caché, sólo lectura b y c para una iteración se mueve 12 bytes de datos, y el ingenuo puede iniciar una nueva iteración nunca tercer ciclo con la canalización.

Asumiendo que nada menos que completo el almacenamiento en caché a hace que la aritmética y el número de instrucciones, incluso menos de un cuello de botella.No me sorprendería si la mayoría de la aceleración en mi código proviene de la emisión de menos cargas para b y c así más espacio libre para el seguimiento y especular pasado caché se pierde en a.

Más amplia de hardware que podría hacer más diferencia.Un Nehalem sistema de ejecución de tres canales de memoria DDR3-1333 sería necesario para mover 3*10667/2.66 = 12.6 bytes por ciclo para saturar el ancho de banda de memoria.Que sería imposible para un solo hilo si a encaja en caché - pero a los 64 bytes de una línea de caché se pierde en el vector agregar rápidamente - sólo uno de los cuatro cargas en mi bucle que faltan en los cachés trae el promedio de ancho de banda requerido para 16 bytes/ciclo.

Answer 4

respuesta corta no. Respuesta larga sí, pero no de manera eficiente. Usted incurrirá en la pena para hacer cargas no alineados lo que eliminará cualquier tipo de beneficio. A menos que usted puede 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 usted sabe de antemano lo que los índices son, su mejor que es para desenrollar y especificar los índices explícitos. Hice algo similar usando especialización de plantilla y 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, de baja carga 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 función, lo desea, puede volver a comprobar. También, el uso de palabras clave con restringir los punteros si usted sabe que no puede haber problemas de aliasing. Esto permitirá que el compilador para ser mucho más agresivo.

Answer 5

Esto no va a vectorizar como es, debido a la doble indirección de los índices de matriz. Puesto que usted está trabajando con dobles hay poco o nada que ganar de SSE, especialmente en lo que la mayoría de las CPU modernas tienen 2 FPU de todos modos.