c optimization assembly x86-64 sse

Producto de punto rápido de un vector de bits y un vector de punto flotante



optimization assembly (7)

Estoy tratando de calcular el producto de puntos entre un vector flotante y un vector de bits de la manera más eficiente en un i7. En realidad, estoy haciendo esta operación en vectores de 128 o 256 dimensiones, pero para ilustrar, permítanme escribir el código para 64 dimensiones para ilustrar el problema:

// a has 64 elements. b is a bitvector of 64 dimensions. float dot(float *restrict a, uint64_t b) { float sum = 0; for(int i=0; b && i<64; i++, b>>=1) { if (b & 1) sum += a[i]; } return sum; }

Esto funciona, por supuesto, pero el problema es que este es el momento crítico de todo el programa (consume el 95% del tiempo de CPU de una ejecución de 50 minutos), así que necesito desesperadamente hacerlo más rápido.

Mi conjetura es que la ramificación de arriba es el asesino del juego (evita la ejecución fuera de orden, causa una mala predicción de la ramificación). No estoy seguro si las instrucciones vectoriales podrían ser usadas y útiles aquí. Al usar gcc 4.8 con -std = c99 -march = native -mtune = native -Ofast -funroll-loops, actualmente obtengo esta salida

movl $4660, %edx movl $5, %ecx xorps %xmm0, %xmm0 .p2align 4,,10 .p2align 3 .L4: testb $1, %cl je .L2 addss (%rdx), %xmm0 .L2: leaq 4(%rdx), %rax shrq %rcx testb $1, %cl je .L8 addss 4(%rdx), %xmm0 .L8: shrq %rcx testb $1, %cl je .L9 addss 4(%rax), %xmm0 .L9: shrq %rcx testb $1, %cl je .L10 addss 8(%rax), %xmm0 .L10: shrq %rcx testb $1, %cl je .L11 addss 12(%rax), %xmm0 .L11: shrq %rcx testb $1, %cl je .L12 addss 16(%rax), %xmm0 .L12: shrq %rcx testb $1, %cl je .L13 addss 20(%rax), %xmm0 .L13: shrq %rcx testb $1, %cl je .L14 addss 24(%rax), %xmm0 .L14: leaq 28(%rax), %rdx shrq %rcx cmpq $4916, %rdx jne .L4 ret

Editar Está bien permutar los datos (siempre que la permutación sea la misma para todos los parámetros), el orden no importa.

Me pregunto si hay algo que funcione a una velocidad> 3x del código SSE2 de Chris Dodd.

Nueva nota: ¡El código AVX / AVX2 también es bienvenido!

Edición 2 Dado un bitvector, tengo que multiplicarlo por 128 (o 256, si es de 256 bits) vectores de coma flotante diferentes (por lo que también está bien involucrar más que un solo vector de flotación a la vez). Este es todo el proceso. Cualquier cosa que acelere todo el proceso también es bienvenida!

Aquí hay algunas cosas para probar

Intenta que el compilador use CMOV lugar de una rama. ( Tenga en cuenta que usar una unión de esta manera está bien definido en C11 pero no definido en C ++ 11 ) .

union { int i; float f; } u; u.i = 0; if (b & 1) { u.f = a[i]; } sum += u.f;

Usa una multiplicación en lugar de una rama.

sum += (b & 1) * a[i];

Mantenga varias sumas y agréguelas al final para reducir las dependencias del flujo de datos. (Puede combinar cualquiera de las sugerencias anteriores con esta).

float sum0 = 0, sum1 = 0, sum2 = 0, sum3 = 0; for (int i = 0; i < 64; i += 4; b >>= 4) { if (b & 1) sum0 += a[i]; if (b & 2) sum1 += a[i+1]; if (b & 4) sum2 += a[i+2]; if (b & 8) sum3 += a[i+3]; } return sum0 + sum1 + sum2 + sum3;

Reduzca el número de ramas procesando varios bits a la vez:

for (int i = 0; i < 64; i += 4, b >>= 4) { switch (b & 0xf) { case 0: break; case 1: sum += a[i]; break; case 2: sum += a[i + 1]; break; case 3: sum += a[i] + a[i+1]; break; case 4: sum += a[i+2]; break; // etc. for cases up to and including 15 } }

Podría mantener varias sumas y, para cada suma, procesar varios bits a la vez. En ese caso, probablemente querrá usar una macro o una función en línea e invocarla cuatro veces.

Descubrí que el ensamblado generado para el código de Chris Dodd es muy dependiente del compilador; clang convierte en un bucle mientras que gcc (4.6 y 4.7) e Intel icc (12.xy 13.x) desenrollan el bucle. Aún así, uno puede reducir las dependencias (la necesidad de esperar la sum += anterior sum += ) convirtiéndola en un mapa para reducir,

float dot(__m128 *a, uint64_t b) { __m128 sum[8]; int i; for (i = 0; i < 8; i++) { sum[i] = _mm_add_ps( _mm_and_ps(a[2*i], mask[b & 0xf].xmm), _mm_and_ps(a[2*i+1], mask[(b & 0xf0) >> 4].xmm)); b >>= 8; } for (i = 0; i < 4; i++) { sum[i] = _mm_add_ps(sum[2*i], sum[2*i+1]); } sum[0] = _mm_add_ps(sum[0], sum[1]); sum[2] = _mm_add_ps(sum[2], sum[3]); sum[0] = _mm_add_ps(sum[0], sum[2]); sum[0] = _mm_hadd_ps(sum[0], sum[0]); sum[0] = _mm_hadd_ps(sum[0], sum[0]); i = _mm_extract_ps(sum[0], 0); return *((float*)(&i)); }

Esto crea un ensamblaje claramente inferior con clang (que almacena la sum[] en la pila) pero mejor código (sin dependencias en addps posteriores) con gcc e icc . De manera bastante interesante, solo gcc tiene la idea al final de que la float inferior en sum[0] puede devolverse en el lugar ...

Buen ejercicio sobre cómo ajustar para compiladores específicos ...

La mejor opción es utilizar instrucciones SSE ps que operen en 4 flotadores a la vez. Puede aprovechar el hecho de que un flotante 0.0 es de 0 bits para usar una instrucción andps para enmascarar los elementos no deseados:

#include <stdint.h> #include <xmmintrin.h> union { uint32_t i[4]; __m128 xmm; } mask[16] = { { 0, 0, 0, 0 }, { ~0, 0, 0, 0 }, { 0, ~0, 0, 0 }, { ~0, ~0, 0, 0 }, { 0, 0, ~0, 0 }, { ~0, 0, ~0, 0 }, { 0, ~0, ~0, 0 }, { ~0, ~0, ~0, 0 }, { 0, 0, 0, ~0 }, { ~0, 0, 0, ~0 }, { 0, ~0, 0, ~0 }, { ~0, ~0, 0, ~0 }, { 0, 0, ~0, ~0 }, { ~0, 0, ~0, ~0 }, { 0, ~0, ~0, ~0 }, { ~0, ~0, ~0, ~0 }, }; float dot(__m128 *a, uint64_t b) { __m128 sum = { 0.0 }; for (int i = 0; i < 16; i++, b>>=4) sum += _mm_and_ps(a[i], mask[b&0xf].xmm); return sum[0] + sum[1] + sum[2] + sum[3]; }

Si esperas que haya una gran cantidad de 0 en la máscara, podría ser más rápido hacer un corto de selección de los 0:

for (int i = 0; b; i++, b >>= 4) if (b & 0xf) sum += _mm_and_ps(a[i], mask[b&0xf].xmm);

pero si b es aleatorio, esto será más lento.

Para ampliar la respuesta de Aki Suihkonen, remodelar la cadena de bits es útil para mover flotadores condicionalmente. En la solución a continuación, una permutación de bits de dos etapas con las instrucciones SSE PMOVMASKB y PSHUFB, más la instrucción BLENDVPS se ha utilizado para lograr 1.25 elementos manejados / ciclos en un Core 2 Duo 2.26GHz, que es 20 veces la velocidad de mi referencia Código C

[EDIT: Se agregó una implementación AVX2. El rendimiento es desconocido porque no puedo probarlo yo mismo, pero se espera que duplique la velocidad. ]

Aquí está mi implementación y testbench, la explicación sigue.

SSE4.1 (antiguo, para AVX2 ver más abajo)

Código

/* Includes */ #include <stdio.h> #include <stdint.h> #include <stdlib.h> #include <smmintrin.h> /* SSE 4.1 */ #include <time.h> /* Defines */ #define ALIGNTO(n) __attribute__((aligned(n))) #define USE_PINSRW 1 #define NUM_ITERS 2260000 /** * Bit mask shuffle. * * This version uses a loop to store eight u16 and reloads them as one __m128i. */ __m128 bitMaskShuffleStoreAndReload(__m128i mask){ const __m128i perm ALIGNTO(16) = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4, 11, 3, 10, 2, 9, 1, 8, 0); int i; uint16_t interMask[8] ALIGNTO(16); /* Shuffle bitmask */ /* Stage 1 */ for(i=7;i>=0;i--){ interMask[i] = _mm_movemask_epi8(mask); mask = _mm_slli_epi32(mask, 1); } /* Stage 2 */ return _mm_castsi128_ps( _mm_shuffle_epi8( _mm_load_si128((const __m128i*)interMask), perm) ); } /** * Bit mask shuffle. * * This version uses the PINSTRW instruction. */ __m128 bitMaskShufflePINSRW(__m128i mask){ const __m128i perm ALIGNTO(16) = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4, 11, 3, 10, 2, 9, 1, 8, 0); __m128i imask ALIGNTO(16); /* Shuffle bitmask */ /* Stage 1 */ imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 7); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 6); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 5); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 4); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 3); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 2); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 1); mask = _mm_slli_epi16(mask, 1); imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 0); /* Stage 2 */ return _mm_castsi128_ps( _mm_shuffle_epi8( imask, perm) ); } /** * SSE 4.1 implementation. */ float dotSSE41(__m128 f[32], unsigned char maskArg[16]){ int i, j, k; __m128i mask ALIGNTO(16) = _mm_load_si128((const __m128i*)maskArg); __m128 shufdMask ALIGNTO(16); __m128 zblended ALIGNTO(16); __m128 sums ALIGNTO(16) = _mm_setzero_ps(); float sumsf[4] ALIGNTO(16); /* Shuffle bitmask */ #if USE_PINSRW shufdMask = bitMaskShufflePINSRW(mask); #else shufdMask = bitMaskShuffleStoreAndReload(mask); #endif /* Dot product */ for(i=1;i>=0;i--){ for(j=1;j>=0;j--){ for(k=7;k>=0;k--){ zblended = _mm_setzero_ps(); zblended = _mm_blendv_ps(zblended, f[i*16+j+k*2], shufdMask); sums = _mm_add_ps(sums, zblended); shufdMask = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(shufdMask), 1)); } } } /* Final Summation */ _mm_store_ps(sumsf, sums); return sumsf[0] + sumsf[1] + sumsf[2] + sumsf[3]; } /** * Reference C implementation */ float dotRefC(float f[128], unsigned char mask[16]){ float sum = 0.0; int i; for(i=0;i<128;i++){ sum += ((mask[i>>3]>>(i&7))&1) ? f[i] : 0.0; } return sum; } /** * Main */ int main(void){ /* Variables */ /* Loop Counter */ int i; /* Data to process */ float data[128] ALIGNTO(16); unsigned char mask[16] ALIGNTO(16); float refCSum, sseSum; /* Time tracking */ clock_t t1, t2, t3; double refCTime, sseTime; /* Initialize mask and float arrays with some random data. */ for(i=0;i<128;i++){ if(i<16) mask[i]=rand(); data[i] = random(); } /* RUN TESTS */ t1 = clock(); for(i=0;i<NUM_ITERS;i++){ refCSum = dotRefC(data, mask); } t2 = clock(); for(i=0;i<NUM_ITERS;i++){ sseSum = dotSSE41((__m128*)data, mask); } t3 = clock(); /* Compute time taken */ refCTime = (double)(t2-t1)/CLOCKS_PER_SEC; sseTime = (double)(t3-t2)/CLOCKS_PER_SEC; /* Print out results */ printf("Results:/n" "RefC: Time: %f Value: %f/n" "SSE: Time: %f Value: %f/n", refCTime, refCSum, sseTime, sseSum); return 0; }

Explicación

BLENDVPS usa el bit superior en los cuatro carriles de 32 bits del registro XMM0 de 128 bits para determinar si se mueve o no para mover el valor en el carril correspondiente de su operando de origen a su operando de destino. Al cargar datos con MOVAPS, se obtienen 4 flotadores consecutivos: por ejemplo, los flotadores 8º, 9º, 10º y 11º. Por supuesto, su selección o deselección debe ser controlada por el conjunto correspondiente de bits: por ejemplo, los bits 8, 9, 10 y 11 en la cadena de bits.

El problema es que cuando se carga la máscara por primera vez, los bits de estos conjuntos están uno junto al otro (en las posiciones 8, 9, 10 y 11), cuando en realidad deben estar separados por 32 bits; Recuerde, en algún punto tendrán que ocupar la posición del bit 31 de cada carril (las posiciones 31, 63, 95 y 127 en el registro XMM0).

Lo que idealmente sucedería es una transposición de bits que traiga los bits 0, 4, 8, 12, ... en el carril cero, los bits 1, 5, 9, 13, ... en el carril uno, los bits 2, 6, 10, 14 , ... en el carril dos y los bits 3, 7, 11, 15, ... en el carril tres. Por lo tanto, todos los conjuntos de 4 bits que antes eran contiguos ahora están separados por 32 bits, uno en cada uno de los cuatro carriles de 32 bits. Luego, todo lo que necesita es un bucle que se repite 32 veces, cada vez cambiando a la posición superior de bits de cada carril un nuevo conjunto de 4 bits.

Desafortunadamente, x86 no está dotado con buenas instrucciones de manipulación de bits, por lo que, a falta de una forma limpia de hacer una transposición perfecta, el compromiso aquí es el razonable.

En la máscara, los 128 bits.

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127

son permutados, por ocho PMOVMASKB y ocho instrucciones PSLLW, primero a

0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 1 9 17 25 33 41 49 57 65 73 81 89 97 105 113 121 2 10 18 26 34 42 50 58 66 74 82 90 98 106 114 122 3 11 19 27 35 43 51 59 67 75 83 91 99 107 115 123 4 12 20 28 36 44 52 60 68 76 84 92 100 108 116 124 5 13 21 29 37 45 53 61 69 77 85 93 101 109 117 125 6 14 22 30 38 46 54 62 70 78 86 94 102 110 118 126 7 15 23 31 39 47 55 63 71 79 87 95 103 111 119 127

y por una sola instrucción PSHUFB para

0 8 16 24 32 40 48 56 4 12 20 28 36 44 52 60 64 72 80 88 96 104 112 120 68 76 84 92 100 108 116 124 1 9 17 25 33 41 49 57 5 13 21 29 37 45 53 61 65 73 81 89 97 105 113 121 69 77 85 93 101 109 117 125 2 10 18 26 34 42 50 58 6 14 22 30 38 46 54 62 66 74 82 90 98 106 114 122 70 78 86 94 102 110 118 126 3 11 19 27 35 43 51 59 7 15 23 31 39 47 55 63 67 75 83 91 99 107 115 123 71 79 87 95 103 111 119 127

. Ahora iteramos en cuatro "corridas", cada una de las cuales contiene ocho conjuntos de cuatro bits distribuidos a intervalos de 32 bits (como deseamos), utilizando estos conjuntos como el control de máscara para BLENDVPS. La incomodidad inherente de la mezcla de bits conduce al bucle anidado triplemente de aspecto incómodo en dotSSE41() , pero con

clang -Ofast -ftree-vectorize -finline-functions -funroll-loops -msse4.1 -mssse3 dot.c -o dottest

Los bucles se desenrollan de todos modos. Las iteraciones del bucle interno consisten en 16 repeticiones de

blendvps 0x90(%rsi),%xmm1 addps %xmm4,%xmm1 pslld $0x1,%xmm2 movdqa %xmm2,%xmm0 xorps %xmm4,%xmm4

.

Dejando de lado, no pude precisar cuál de mis versiones aleatorias de dos bits fue la más rápida, así que di ambas implementaciones en mi respuesta.

AVX2 (nuevo, pero no probado)

Código

/* Includes */ #include <stdio.h> #include <stdint.h> #include <stdlib.h> #include <immintrin.h> /* AVX2 */ #include <time.h> /* Defines */ #define ALIGNTO(n) __attribute__((aligned(n))) #define NUM_ITERS 2260000 /** * Bit mask shuffle. * * This version uses the PINSTRW instruction. */ __m256 bitMaskShufflePINSRW(__m256i mask){ __m256i imask ALIGNTO(32); /* Shuffle bitmask */ imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 7); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 6); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 5); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 4); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 3); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 2); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 1); mask = _mm256_slli_epi32(mask, 1); imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 0); /* Return bitmask */ return _mm256_castsi256_ps(imask); } /** * AVX2 implementation. */ float dotAVX2(__m256 f[16], unsigned char maskArg[16]){ int i, j, k; /* Use _mm_loadu_si128 */ __m256i mask ALIGNTO(32) = _mm256_castsi128_si256(_mm_load_si128((const __m128i*)maskArg)); __m256 shufdMask ALIGNTO(32); __m256 zblended ALIGNTO(32); __m256 sums ALIGNTO(32) = _mm256_setzero_ps(); float sumsf[8] ALIGNTO(32); /* Shuffle bitmask */ shufdMask = bitMaskShufflePINSRW(mask); shufdMask = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(shufdMask), 16)); /* Dot product */ for(i=15;i>=0;i--){ zblended = _mm256_setzero_ps(); /* Replace f[i] with _mm256_loadu_ps((float*)&f[i]) */ zblended = _mm256_blendv_ps(zblended, f[i], shufdMask); sums = _mm256_add_ps(sums, zblended); shufdMask = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(shufdMask), 1)); } /* Final Summation */ _mm256_store_ps(sumsf, sums); return sumsf[0] + sumsf[1] + sumsf[2] + sumsf[3] + sumsf[4] + sumsf[5] + sumsf[6] + sumsf[7]; } /** * Reference C implementation */ float dotRefC(float f[128], unsigned char mask[16]){ float sum = 0.0; int i; for(i=0;i<128;i++){ sum += ((mask[i>>3]>>(i&7))&1) ? f[i] : 0.0; } return sum; } /** * Main */ int main(void){ /* Variables */ /* Loop Counter */ int i; /* Data to process */ float data[128] ALIGNTO(32); unsigned char mask[16] ALIGNTO(32); float refCSum, sseSum; /* Time tracking */ clock_t t1, t2, t3; double refCTime, sseTime; /* Initialize mask and float arrays with some random data. */ for(i=0;i<128;i++){ if(i<16) mask[i]=rand(); data[i] = random(); } /* RUN TESTS */ t1 = clock(); for(i=0;i<NUM_ITERS;i++){ refCSum = dotRefC(data, mask); } t2 = clock(); for(i=0;i<NUM_ITERS;i++){ sseSum = dotAVX2((__m256*)data, mask); } t3 = clock(); /* Compute time taken */ refCTime = (double)(t2-t1)/CLOCKS_PER_SEC; sseTime = (double)(t3-t2)/CLOCKS_PER_SEC; /* Print out results */ printf("Results:/n" "RefC: Time: %f Value: %f/n" "SSE: Time: %f Value: %f/n", refCTime, refCSum, sseTime, sseSum); return 0; }

Explicación

Se utiliza el mismo concepto que para SSE4.1. La diferencia es que ahora estamos procesando 8 flotantes a la vez y haciendo uso de los registros de 256 bits de AVX2 y PMOVMASKB de los registros ymm (que reúnen 256/8 = 32 bits). Debido a esto, ahora tenemos una mezcla de máscara de bits más simple y un bucle mucho más simple.

En la máscara, los 256 bits.

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255

se permutan, usando 8 PMOVMASKB y 8 instrucciones PSLLW, para

0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120 128 136 144 152 160 168 176 184 192 200 208 216 224 232 240 248 1 9 17 25 33 41 49 57 65 73 81 89 97 105 113 121 129 137 145 153 161 169 177 185 193 201 209 217 225 233 241 249 2 10 18 26 34 42 50 58 66 74 82 90 98 106 114 122 130 138 146 154 162 170 178 186 194 202 210 218 226 234 242 250 3 11 19 27 35 43 51 59 67 75 83 91 99 107 115 123 131 139 147 155 163 171 179 187 195 203 211 219 227 235 243 251 4 12 20 28 36 44 52 60 68 76 84 92 100 108 116 124 132 140 148 156 164 172 180 188 196 204 212 220 228 236 244 252 5 13 21 29 37 45 53 61 69 77 85 93 101 109 117 125 133 141 149 157 165 173 181 189 197 205 213 221 229 237 245 253 6 14 22 30 38 46 54 62 70 78 86 94 102 110 118 126 134 142 150 158 166 174 182 190 198 206 214 222 230 238 246 254 7 15 23 31 39 47 55 63 71 79 87 95 103 111 119 127 135 143 151 159 167 175 183 191 199 207 215 223 231 239 247 255

. Para los productos de punto de bit con flotador de 128 elementos, luego iteramos en paralelo en ocho conjuntos de 16 elementos. Esta implementación se puede ampliar fácilmente para DP de 256 elementos iterando en conjuntos de 32 elementos. Ahora solo se requiere un bucle.

Específicamente, para cambiar esto para que funcione con productos de puntos de 256 elementos,

  • Doble el tamaño de los argumentos de la función. __m256 f[32], unsigned char maskArg[32] .
  • Intercambie la carga de la máscara ( = _mm256_castsi128_si256(_mm_load_si128((const __m128i*)maskArg)); ) con = _mm256_load_si256((const __m256i*)maskArg); .
  • Elimine el desplazamiento compensatorio a la izquierda por 16 justo debajo de la llamada a bitMaskShufflePINSRW .
  • Ejecutar el bucle hacia abajo desde el conjunto 31 en lugar de 15: for(i=31;i>=0;i--)

No puedo probar ni ejecutar el código ya que mi CPU es SSE4.1 , pero Clang con

clang -Ofast -ftree-vectorize -finline-functions -funroll-loops -mavx2 -msse4.1 -mssse3 dotavx2.c -o dottest

compilado limpiamente (puede obtener un código más rápido sin desenrollar), produciendo esto:

(gdb) disas dotAVX2 Dump of assembler code for function dotAVX2: 0x0000000100001730 <+0>: push %rbp 0x0000000100001731 <+1>: mov %rsp,%rbp 0x0000000100001734 <+4>: vmovdqa (%rsi),%xmm0 0x0000000100001738 <+8>: vpslld $0x1,%ymm0,%ymm1 0x000000010000173d <+13>: vpslld $0x1,%ymm1,%ymm2 0x0000000100001742 <+18>: vpmovmskb %ymm2,%eax 0x0000000100001746 <+22>: vpslld $0x1,%ymm2,%ymm2 0x000000010000174b <+27>: vpmovmskb %ymm2,%ecx 0x000000010000174f <+31>: vxorps %ymm3,%ymm3,%ymm3 0x0000000100001753 <+35>: vmovd %ecx,%xmm4 0x0000000100001757 <+39>: vpinsrd $0x1,%eax,%xmm4,%xmm4 0x000000010000175d <+45>: vpmovmskb %ymm1,%eax 0x0000000100001761 <+49>: vpinsrd $0x2,%eax,%xmm4,%xmm1 0x0000000100001767 <+55>: vpslld $0x1,%ymm2,%ymm2 0x000000010000176c <+60>: vpslld $0x1,%ymm2,%ymm4 0x0000000100001771 <+65>: vpslld $0x1,%ymm4,%ymm5 0x0000000100001776 <+70>: vpmovmskb %ymm0,%eax 0x000000010000177a <+74>: vpinsrd $0x3,%eax,%xmm1,%xmm0 0x0000000100001780 <+80>: vpmovmskb %ymm5,%eax 0x0000000100001784 <+84>: vpslld $0x1,%ymm5,%ymm1 0x0000000100001789 <+89>: vpmovmskb %ymm1,%ecx 0x000000010000178d <+93>: vmovd %ecx,%xmm1 0x0000000100001791 <+97>: vpinsrd $0x1,%eax,%xmm1,%xmm1 0x0000000100001797 <+103>: vpmovmskb %ymm4,%eax 0x000000010000179b <+107>: vpinsrd $0x2,%eax,%xmm1,%xmm1 0x00000001000017a1 <+113>: vpmovmskb %ymm2,%eax 0x00000001000017a5 <+117>: vpinsrd $0x3,%eax,%xmm1,%xmm1 0x00000001000017ab <+123>: vinserti128 $0x1,%xmm0,%ymm1,%ymm0 0x00000001000017b1 <+129>: vpslld $0x10,%ymm0,%ymm0 0x00000001000017b6 <+134>: vblendvps %ymm0,0x1e0(%rdi),%ymm3,%ymm1 0x00000001000017c0 <+144>: vpslld $0x1,%ymm0,%ymm0 0x00000001000017c5 <+149>: vblendvps %ymm0,0x1c0(%rdi),%ymm3,%ymm2 0x00000001000017cf <+159>: vaddps %ymm2,%ymm1,%ymm1 0x00000001000017d3 <+163>: vpslld $0x1,%ymm0,%ymm0 0x00000001000017d8 <+168>: vblendvps %ymm0,0x1a0(%rdi),%ymm3,%ymm2 0x00000001000017e2 <+178>: vaddps %ymm2,%ymm1,%ymm1 0x00000001000017e6 <+182>: vpslld $0x1,%ymm0,%ymm0 0x00000001000017eb <+187>: vblendvps %ymm0,0x180(%rdi),%ymm3,%ymm2 0x00000001000017f5 <+197>: vaddps %ymm2,%ymm1,%ymm1 0x00000001000017f9 <+201>: vpslld $0x1,%ymm0,%ymm0 0x00000001000017fe <+206>: vblendvps %ymm0,0x160(%rdi),%ymm3,%ymm2 0x0000000100001808 <+216>: vaddps %ymm2,%ymm1,%ymm1 0x000000010000180c <+220>: vpslld $0x1,%ymm0,%ymm0 0x0000000100001811 <+225>: vblendvps %ymm0,0x140(%rdi),%ymm3,%ymm2 0x000000010000181b <+235>: vaddps %ymm2,%ymm1,%ymm1 0x000000010000181f <+239>: vpslld $0x1,%ymm0,%ymm0 0x0000000100001824 <+244>: vblendvps %ymm0,0x120(%rdi),%ymm3,%ymm2 0x000000010000182e <+254>: vaddps %ymm2,%ymm1,%ymm1 0x0000000100001832 <+258>: vpslld $0x1,%ymm0,%ymm0 0x0000000100001837 <+263>: vblendvps %ymm0,0x100(%rdi),%ymm3,%ymm2 0x0000000100001841 <+273>: vaddps %ymm2,%ymm1,%ymm1 0x0000000100001845 <+277>: vpslld $0x1,%ymm0,%ymm0 0x000000010000184a <+282>: vblendvps %ymm0,0xe0(%rdi),%ymm3,%ymm2 0x0000000100001854 <+292>: vaddps %ymm2,%ymm1,%ymm1 0x0000000100001858 <+296>: vpslld $0x1,%ymm0,%ymm0 0x000000010000185d <+301>: vblendvps %ymm0,0xc0(%rdi),%ymm3,%ymm2 0x0000000100001867 <+311>: vaddps %ymm2,%ymm1,%ymm1 0x000000010000186b <+315>: vpslld $0x1,%ymm0,%ymm0 0x0000000100001870 <+320>: vblendvps %ymm0,0xa0(%rdi),%ymm3,%ymm2 0x000000010000187a <+330>: vaddps %ymm2,%ymm1,%ymm1 0x000000010000187e <+334>: vpslld $0x1,%ymm0,%ymm0 0x0000000100001883 <+339>: vblendvps %ymm0,0x80(%rdi),%ymm3,%ymm2 0x000000010000188d <+349>: vaddps %ymm2,%ymm1,%ymm1 0x0000000100001891 <+353>: vpslld $0x1,%ymm0,%ymm0 0x0000000100001896 <+358>: vblendvps %ymm0,0x60(%rdi),%ymm3,%ymm2 0x000000010000189d <+365>: vaddps %ymm2,%ymm1,%ymm1 0x00000001000018a1 <+369>: vpslld $0x1,%ymm0,%ymm0 0x00000001000018a6 <+374>: vblendvps %ymm0,0x40(%rdi),%ymm3,%ymm2 0x00000001000018ad <+381>: vaddps %ymm2,%ymm1,%ymm1 0x00000001000018b1 <+385>: vpslld $0x1,%ymm0,%ymm0 0x00000001000018b6 <+390>: vblendvps %ymm0,0x20(%rdi),%ymm3,%ymm2 0x00000001000018bd <+397>: vaddps %ymm2,%ymm1,%ymm1 0x00000001000018c1 <+401>: vpslld $0x1,%ymm0,%ymm0 0x00000001000018c6 <+406>: vblendvps %ymm0,(%rdi),%ymm3,%ymm0 0x00000001000018cc <+412>: vaddps %ymm0,%ymm1,%ymm0 0x00000001000018d0 <+416>: vpshufd $0x1,%xmm0,%xmm1 0x00000001000018d5 <+421>: vaddss %xmm1,%xmm0,%xmm1 0x00000001000018d9 <+425>: vmovhlps %xmm0,%xmm0,%xmm2 0x00000001000018dd <+429>: vaddss %xmm1,%xmm2,%xmm1 0x00000001000018e1 <+433>: vpshufd $0x3,%xmm0,%xmm2 0x00000001000018e6 <+438>: vaddss %xmm1,%xmm2,%xmm1 0x00000001000018ea <+442>: vextracti128 $0x1,%ymm0,%xmm0 0x00000001000018f0 <+448>: vaddss %xmm1,%xmm0,%xmm1 0x00000001000018f4 <+452>: vpshufd $0x1,%xmm0,%xmm2 0x00000001000018f9 <+457>: vaddss %xmm1,%xmm2,%xmm1 0x00000001000018fd <+461>: vpshufd $0x3,%xmm0,%xmm2 0x0000000100001902 <+466>: vmovhlps %xmm0,%xmm0,%xmm0 0x0000000100001906 <+470>: vaddss %xmm1,%xmm0,%xmm0 0x000000010000190a <+474>: vaddss %xmm0,%xmm2,%xmm0 0x000000010000190e <+478>: pop %rbp 0x000000010000190f <+479>: vzeroupper 0x0000000100001912 <+482>: retq End of assembler dump.

Como podemos ver, el núcleo tiene 3 instrucciones (vblendvps, vaddps, vpslld) ahora.

Si se permite una permutación ligeramente diferente en los datos de datos float data[128] , o se realiza la permutación correspondiente en la máscara de bits para la __m128 mask; , uno puede mejorar ligeramente el algoritmo sugerido por Chris Dodd arriba. (Sin contar el tiempo requerido para permutar la máscara, esta implementación (+ sobrecarga) es aproximadamente un 25% más rápida). Por supuesto, esto es solo un borrador rápido de mi idea que se proporciona en los comentarios.

union { unsigned int i[4]; float f[4]; __m128 xmm; } mask = { 0xFF00FF00, 0xF0F0F0F0, 0xCCCCCCCC, 0xAAAAAAAA }; float dot2(__m128 *a, __m128 mask); // 20M times 1.161s float dotref(__m128 *a, unsigned int *mask) // 20M times 8.174s { float z=0.0f; int i; for (i=0;i<32;i++) { if (mask[0] & (0x80000000U >> i)) z+= a[i][0]; if (mask[1] & (0x80000000U >> i)) z+= a[i][1]; if (mask[2] & (0x80000000U >> i)) z+= a[i][2]; if (mask[3] & (0x80000000U >> i)) z+= a[i][3]; } return z; }

La implementación del ensamblador correspondiente sería:

dot2: // warm up stage: fill in initial data and // set up registers pxor %xmm1, %xmm1 ;; // clear partial sum1 pxor %xmm2, %xmm2 ;; // clear partial sum2 movaps (%rdi), %xmm3 ;; // register warm up stage1 movaps 16(%rdi), %xmm4 ;; // next 4 values pxor %xmm5, %xmm5 pxor %xmm6, %xmm6 lea 32(%rdi), %rdi movl $16, %ecx ;; // process 2x4 items per iteration (total=128) a: ;; // inner loop -- 2 independent data paths blendvps %xmm3, %xmm5 pslld $1, %xmm0 movaps (%rdi), %xmm3 blendvps %xmm4, %xmm6 pslld $1, %xmm0 movaps 16(%rdi), %xmm4 addps %xmm5, %xmm1 pxor %xmm5, %xmm5 addps %xmm6, %xmm2 pxor %xmm6, %xmm6 lea 32(%rdi), %rdi loop a ;; // cool down stage: gather results (xmm0 = xmm1+xmm2) ;; // in beautiful world this stage is interleaved ;; // with the warm up stage of the next block addps %xmm2, %xmm1 movaps %xmm1, %xmm2 movaps %xmm1, %xmm0 shufps $85, %xmm1, %xmm2 addss %xmm2, %xmm0 movaps %xmm1, %xmm2 unpckhps %xmm1, %xmm2 shufps $255, %xmm1, %xmm1 addss %xmm2, %xmm0 addss %xmm1, %xmm0 ret

Si tiene i7, entonces tiene SSE4.1 y puede usar el _mm_dp_ps intrinsic, que hace un producto de puntos. Con tu código de ejemplo se vería como

#include <stdint.h> #include <immintrin.h> const float fltmask[][4] = {{0, 0, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}, {0, 0, 1, 1}, {0, 1, 0, 0}, {0, 1, 0, 1}, {0, 1, 1, 0}, {0, 1, 1, 1}, {1, 0, 0, 0}, {1, 0, 0, 1}, {1, 0, 1, 0}, {1, 0, 1, 1}, {1, 1, 0, 0}, {1, 1, 0, 1}, {1, 1, 1, 0}, {1, 1, 1, 1}}; // a has 64 elements. b is a bitvector of 64 dimensions. float dot(float * restrict a, uint64_t b) { int i; float sum = 0; for(i=0; b && i<64; i+=4,b>>=4) { __m128 t0 = _mm_load_ps (a); a += 4; __m128 t1 = _mm_load_ps (fltmask[b & 15]); sum += _mm_cvtss_f32 (_mm_dp_ps (t0, t1, 15)); } return sum; }

PD.Las matrices ay fltmaskmejor se alinean 16 bytes!

PPS.Cuando se compila con gcc -std=c99 -msse4 -O2el bucle se ve así:

.L3: movq %rdx, %rax movaps (%rcx), %xmm1 shrq $4, %rdx andl $15, %eax addq $16, %rcx addl $4, %r8d salq $4, %rax testq %rdx, %rdx dpps $15, (%r9,%rax), %xmm1 addss %xmm1, %xmm0 jne .L13

Y con -O3su desenrollado, por supuesto.

Puedes eliminar la rama así:

for(int i=0; b && i<64; i++, b>>=1) sum += a[i] * (b & 1);

Aunque esto agregará un mult adicional, por lo menos no arruinará tu canalización.

Otra forma de controlar la rama es usarla de la manera que lo eres, pero también usar una macro compiladora. Supongo que en gcc es la likely(if ...)macro. Utilizará la rama, pero de esa manera le está diciendo al compilador que la rama se ejecutará más a menudo, y gcc optimizará más.

Otra optimización que se puede hacer es "almacenar en caché" el producto de puntos. Entonces, en lugar de tener una función para calcular el producto de puntos, tendría una variable que mantiene el producto inicializado en 0 al principio. Y cada vez que inserta / elimina / actualiza un elemento del vector, también actualizaría la variable que contiene el resultado.