c++ - ¿Por qué esta multiplicación SIMD no es más rápida que la multiplicación que no es SIMD?
performance assembly (3)
Quiero agregar otro punto de vista al problema. Las instrucciones SIMD dan un gran aumento de rendimiento si no hay restricciones de memoria. Pero hay demasiadas operaciones de carga y almacenamiento de memoria y muy pocos cálculos de CPU en el ejemplo actual. Así que la CPU está a tiempo para procesar los datos entrantes sin usar SIMD. Si usa datos de otro tipo (flotador de 32 bits, por ejemplo) o algoritmos más complejos, el rendimiento de la memoria no restringirá el rendimiento de la CPU y el uso de SIMD le dará más ventajas.
Supongamos que tenemos una función que multiplica dos matrices de 1000000 dobles cada una. En C / C ++, la función es así:
void mul_c(double* a, double* b)
{
for (int i = 0; i != 1000000; ++i)
{
a[i] = a[i] * b[i];
}
}
El compilador produce el siguiente ensamblaje con -O2 :
mul_c(double*, double*):
xor eax, eax
.L2:
movsd xmm0, QWORD PTR [rdi+rax]
mulsd xmm0, QWORD PTR [rsi+rax]
movsd QWORD PTR [rdi+rax], xmm0
add rax, 8
cmp rax, 8000000
jne .L2
rep ret
Del ensamblaje anterior parece que el compilador usa las instrucciones SIMD, pero solo multiplica un doble en cada iteración. Así que decidí escribir la misma función en ensamblaje en línea, donde aprovecho al máximo el registro xmm0 y multiplico dos dobles de una vez:
void mul_asm(double* a, double* b)
{
asm volatile
(
".intel_syntax noprefix /n/t"
"xor rax, rax /n/t"
"0: /n/t"
"movupd xmm0, xmmword ptr [rdi+rax] /n/t"
"mulpd xmm0, xmmword ptr [rsi+rax] /n/t"
"movupd xmmword ptr [rdi+rax], xmm0 /n/t"
"add rax, 16 /n/t"
"cmp rax, 8000000 /n/t"
"jne 0b /n/t"
".att_syntax noprefix /n/t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
}
Después de medir el tiempo de ejecución individualmente para estas dos funciones, parece que ambas tardan 1 ms en completarse:
> gcc -O2 main.cpp
> ./a.out < input
mul_c: 1 ms
mul_asm: 1 ms
[a lot of doubles...]
Esperaba que la implementación de SIMD fuera al menos dos veces más rápida (0 ms) ya que solo hay la mitad de la cantidad de instrucciones de multiplicaciones / memoria.
Entonces mi pregunta es: ¿por qué la implementación de SIMD no es más rápida que la implementación normal de C / C ++ cuando la implementación de SIMD solo hace la mitad de la cantidad de instrucciones de multiplicaciones / memoria?
Aquí está el programa completo:
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
void mul_c(double* a, double* b)
{
for (int i = 0; i != 1000000; ++i)
{
a[i] = a[i] * b[i];
}
}
void mul_asm(double* a, double* b)
{
asm volatile
(
".intel_syntax noprefix /n/t"
"xor rax, rax /n/t"
"0: /n/t"
"movupd xmm0, xmmword ptr [rdi+rax] /n/t"
"mulpd xmm0, xmmword ptr [rsi+rax] /n/t"
"movupd xmmword ptr [rdi+rax], xmm0 /n/t"
"add rax, 16 /n/t"
"cmp rax, 8000000 /n/t"
"jne 0b /n/t"
".att_syntax noprefix /n/t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
}
int main()
{
struct timeval t1;
struct timeval t2;
unsigned long long time;
double* a = (double*)malloc(sizeof(double) * 1000000);
double* b = (double*)malloc(sizeof(double) * 1000000);
double* c = (double*)malloc(sizeof(double) * 1000000);
for (int i = 0; i != 1000000; ++i)
{
double v;
scanf("%lf", &v);
a[i] = v;
b[i] = v;
c[i] = v;
}
gettimeofday(&t1, NULL);
mul_c(a, b);
gettimeofday(&t2, NULL);
time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
printf("mul_c: %llu ms/n", time);
gettimeofday(&t1, NULL);
mul_asm(b, c);
gettimeofday(&t2, NULL);
time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
printf("mul_asm: %llu ms/n/n", time);
for (int i = 0; i != 1000000; ++i)
{
printf("%lf/t/t/t%lf/n", a[i], b[i]);
}
return 0;
}
También traté de hacer uso de todos los registros xmm (0-7) y eliminar las dependencias de instrucciones para obtener una mejor computación en paralelo:
void mul_asm(double* a, double* b)
{
asm volatile
(
".intel_syntax noprefix /n/t"
"xor rax, rax /n/t"
"0: /n/t"
"movupd xmm0, xmmword ptr [rdi+rax] /n/t"
"movupd xmm1, xmmword ptr [rdi+rax+16] /n/t"
"movupd xmm2, xmmword ptr [rdi+rax+32] /n/t"
"movupd xmm3, xmmword ptr [rdi+rax+48] /n/t"
"movupd xmm4, xmmword ptr [rdi+rax+64] /n/t"
"movupd xmm5, xmmword ptr [rdi+rax+80] /n/t"
"movupd xmm6, xmmword ptr [rdi+rax+96] /n/t"
"movupd xmm7, xmmword ptr [rdi+rax+112] /n/t"
"mulpd xmm0, xmmword ptr [rsi+rax] /n/t"
"mulpd xmm1, xmmword ptr [rsi+rax+16] /n/t"
"mulpd xmm2, xmmword ptr [rsi+rax+32] /n/t"
"mulpd xmm3, xmmword ptr [rsi+rax+48] /n/t"
"mulpd xmm4, xmmword ptr [rsi+rax+64] /n/t"
"mulpd xmm5, xmmword ptr [rsi+rax+80] /n/t"
"mulpd xmm6, xmmword ptr [rsi+rax+96] /n/t"
"mulpd xmm7, xmmword ptr [rsi+rax+112] /n/t"
"movupd xmmword ptr [rdi+rax], xmm0 /n/t"
"movupd xmmword ptr [rdi+rax+16], xmm1 /n/t"
"movupd xmmword ptr [rdi+rax+32], xmm2 /n/t"
"movupd xmmword ptr [rdi+rax+48], xmm3 /n/t"
"movupd xmmword ptr [rdi+rax+64], xmm4 /n/t"
"movupd xmmword ptr [rdi+rax+80], xmm5 /n/t"
"movupd xmmword ptr [rdi+rax+96], xmm6 /n/t"
"movupd xmmword ptr [rdi+rax+112], xmm7 /n/t"
"add rax, 128 /n/t"
"cmp rax, 8000000 /n/t"
"jne 0b /n/t"
".att_syntax noprefix /n/t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
}
Pero todavía funciona a 1 ms, la misma velocidad que la implementación ordinaria de C / C ++.
ACTUALIZACIONES
Según lo sugerido por las respuestas / comentarios, he implementado otra forma de medir el tiempo de ejecución:
#include <stdio.h>
#include <stdlib.h>
void mul_c(double* a, double* b)
{
for (int i = 0; i != 1000000; ++i)
{
a[i] = a[i] * b[i];
}
}
void mul_asm(double* a, double* b)
{
asm volatile
(
".intel_syntax noprefix /n/t"
"xor rax, rax /n/t"
"0: /n/t"
"movupd xmm0, xmmword ptr [rdi+rax] /n/t"
"mulpd xmm0, xmmword ptr [rsi+rax] /n/t"
"movupd xmmword ptr [rdi+rax], xmm0 /n/t"
"add rax, 16 /n/t"
"cmp rax, 8000000 /n/t"
"jne 0b /n/t"
".att_syntax noprefix /n/t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
}
void mul_asm2(double* a, double* b)
{
asm volatile
(
".intel_syntax noprefix /n/t"
"xor rax, rax /n/t"
"0: /n/t"
"movupd xmm0, xmmword ptr [rdi+rax] /n/t"
"movupd xmm1, xmmword ptr [rdi+rax+16] /n/t"
"movupd xmm2, xmmword ptr [rdi+rax+32] /n/t"
"movupd xmm3, xmmword ptr [rdi+rax+48] /n/t"
"movupd xmm4, xmmword ptr [rdi+rax+64] /n/t"
"movupd xmm5, xmmword ptr [rdi+rax+80] /n/t"
"movupd xmm6, xmmword ptr [rdi+rax+96] /n/t"
"movupd xmm7, xmmword ptr [rdi+rax+112] /n/t"
"mulpd xmm0, xmmword ptr [rsi+rax] /n/t"
"mulpd xmm1, xmmword ptr [rsi+rax+16] /n/t"
"mulpd xmm2, xmmword ptr [rsi+rax+32] /n/t"
"mulpd xmm3, xmmword ptr [rsi+rax+48] /n/t"
"mulpd xmm4, xmmword ptr [rsi+rax+64] /n/t"
"mulpd xmm5, xmmword ptr [rsi+rax+80] /n/t"
"mulpd xmm6, xmmword ptr [rsi+rax+96] /n/t"
"mulpd xmm7, xmmword ptr [rsi+rax+112] /n/t"
"movupd xmmword ptr [rdi+rax], xmm0 /n/t"
"movupd xmmword ptr [rdi+rax+16], xmm1 /n/t"
"movupd xmmword ptr [rdi+rax+32], xmm2 /n/t"
"movupd xmmword ptr [rdi+rax+48], xmm3 /n/t"
"movupd xmmword ptr [rdi+rax+64], xmm4 /n/t"
"movupd xmmword ptr [rdi+rax+80], xmm5 /n/t"
"movupd xmmword ptr [rdi+rax+96], xmm6 /n/t"
"movupd xmmword ptr [rdi+rax+112], xmm7 /n/t"
"add rax, 128 /n/t"
"cmp rax, 8000000 /n/t"
"jne 0b /n/t"
".att_syntax noprefix /n/t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
}
unsigned long timestamp()
{
unsigned long a;
asm volatile
(
".intel_syntax noprefix /n/t"
"xor rax, rax /n/t"
"xor rdx, rdx /n/t"
"RDTSCP /n/t"
"shl rdx, 32 /n/t"
"or rax, rdx /n/t"
".att_syntax noprefix /n/t"
: "=a" (a)
:
: "memory", "cc"
);
return a;
}
int main()
{
unsigned long t1;
unsigned long t2;
double* a;
double* b;
a = (double*)malloc(sizeof(double) * 1000000);
b = (double*)malloc(sizeof(double) * 1000000);
for (int i = 0; i != 1000000; ++i)
{
double v;
scanf("%lf", &v);
a[i] = v;
b[i] = v;
}
t1 = timestamp();
mul_c(a, b);
//mul_asm(a, b);
//mul_asm2(a, b);
t2 = timestamp();
printf("mul_c: %lu cycles/n/n", t2 - t1);
for (int i = 0; i != 1000000; ++i)
{
printf("%lf/t/t/t%lf/n", a[i], b[i]);
}
return 0;
}
Cuando ejecuto el programa con esta medida, obtengo este resultado:
mul_c: ~2163971628 cycles
mul_asm: ~2532045184 cycles
mul_asm2: ~5230488 cycles <-- what???
Hay dos cosas que vale la pena tener en cuenta aquí, en primer lugar, el recuento de ciclos varía mucho, y supongo que es porque el sistema operativo permite que otros procesos se ejecuten en medio. ¿Hay alguna manera de evitar eso o solo contar los ciclos mientras se ejecuta mi programa? Además, mul_asm2 produce resultados idénticos en comparación con los otros dos, pero es mucho más rápido, ¿cómo?
Probé el programa de Z boson en mi sistema junto con mis 2 implementaciones y obtuve el siguiente resultado:
> g++ -O2 -fopenmp main.cpp
> ./a.out
mul time 1.33, 18.08 GB/s
mul_SSE time 1.13, 21.24 GB/s
mul_SSE_NT time 1.51, 15.88 GB/s
mul_SSE_OMP time 0.79, 30.28 GB/s
mul_SSE_v2 time 1.12, 21.49 GB/s
mul_v2 time 1.26, 18.99 GB/s
mul_asm time 1.12, 21.50 GB/s
mul_asm2 time 1.09, 22.08 GB/s
Su código asm está realmente bien. Lo que no es es la forma en que lo mides. Como señalé en los comentarios, debes:
a) usa mucho más iteraciones: 1 millón no es nada para la CPU moderna
b) use HPT para la medición
c) use RDTSC o RDTSCP para contar relojes reales de CPU
Además, ¿por qué tiene miedo de optar -O3? No olvides crear código para tu plataforma, así que usa -march = native. Si su CPU es compatible con AVX o AVX2, el compilador aprovechará la oportunidad para producir un código aún mejor.
Lo siguiente es darle al compilador algunas pistas sobre alias y alineación si conoce su código.
Esta es mi versión de tu mul_c : sí, es específica de GCC, pero has demostrado que mul_c GCC
void mul_c(double* restrict a, double* restrict b)
{
a = __builtin_assume_aligned (a, 16);
b = __builtin_assume_aligned (b, 16);
for (int i = 0; i != 1000000; ++i)
{
a[i] = a[i] * b[i];
}
}
Producirá:
mul_c(double*, double*):
xor eax, eax
.L2:
movapd xmm0, XMMWORD PTR [rdi+rax]
mulpd xmm0, XMMWORD PTR [rsi+rax]
movaps XMMWORD PTR [rdi+rax], xmm0
add rax, 16
cmp rax, 8000000
jne .L2
rep ret
Si tiene AVX2 y se asegura de que los datos estén alineados a 32 bytes, se convertirá en
mul_c(double*, double*):
xor eax, eax
.L2:
vmovapd ymm0, YMMWORD PTR [rdi+rax]
vmulpd ymm0, ymm0, YMMWORD PTR [rsi+rax]
vmovapd YMMWORD PTR [rdi+rax], ymm0
add rax, 32
cmp rax, 8000000
jne .L2
vzeroupper
ret
Así que no hay necesidad de un asm hecho a mano si el compilador puede hacerlo por ti;)
Hubo un error importante en la función de tiempo que utilicé para los puntos de referencia anteriores. Esto subestimó groseramente el ancho de banda sin vectorización, así como otras medidas. Además, había otro problema que sobreestimaba el ancho de banda debido a COW en la matriz que se leyó pero no se escribió. Finalmente, el ancho de banda máximo que utilicé fue incorrecto. He actualizado mi respuesta con las correcciones y he dejado la respuesta anterior al final de esta respuesta.
Su operación está vinculada al ancho de banda de la memoria. Esto significa que la CPU está pasando la mayor parte del tiempo esperando lecturas y escrituras lentas. Una excelente explicación para esto se puede encontrar aquí: ¿Por qué vectorizar el bucle no tiene una mejora en el rendimiento ?
Sin embargo, tengo que estar un poco en desacuerdo con una afirmación en esa respuesta.
Por lo tanto, independientemente de cómo esté optimizado (vectorizado, desenrollado, etc.), no será mucho más rápido.
De hecho, la vectorización , desenrollamiento y múltiples hilos pueden aumentar significativamente el ancho de banda incluso en operaciones de límite de ancho de banda de memoria. La razón es que es difícil obtener el ancho de banda de memoria máximo. Una buena explicación para esto se puede encontrar aquí: https://.com/a/25187492/2542702 .
El resto de mi respuesta mostrará cómo la vectorización y los subprocesos múltiples pueden acercarse al ancho de banda de memoria máximo.
Mi sistema de prueba: Ubuntu 16.10, Skylake ([email protected]), 32GB RAM, DDR4 de doble canal a 2400 GHz. El ancho de banda máximo de mi sistema es de 38,4 GB / s.
Del código a continuación, produzco las siguientes tablas. Configuré el número de subprocesos usando OMP_NUM_THREADS, por ejemplo export OMP_NUM_THREADS=4 . La eficiencia es bandwidth/max_bandwidth .
-O2 -march=native -fopenmp
Threads Efficiency
1 59.2%
2 76.6%
4 74.3%
8 70.7%
-O2 -march=native -fopenmp -funroll-loops
1 55.8%
2 76.5%
4 72.1%
8 72.2%
-O3 -march=native -fopenmp
1 63.9%
2 74.6%
4 63.9%
8 63.2%
-O3 -march=native -fopenmp -mprefer-avx128
1 67.8%
2 76.0%
4 63.9%
8 63.2%
-O3 -march=native -fopenmp -mprefer-avx128 -funroll-loops
1 68.8%
2 73.9%
4 69.0%
8 66.8%
Después de varias iteraciones de ejecución debido a las incertidumbres en las mediciones, formé las siguientes conclusiones:
- las operaciones escalares de un solo hilo obtienen más del 50% del ancho de banda.
- dos operaciones escalares con hilos obtienen el mayor ancho de banda.
- las operaciones de vector con un solo hilo son más rápidas que las operaciones escalares de un solo hilo.
- las operaciones SSE de una sola rosca son más rápidas que las operaciones AVX de rosca única.
- desenrollarse no es útil.
- desenrollar las operaciones de un solo hilo es más lento que sin desenrollar.
- más hilos que núcleos (Hyper-Threading) da un ancho de banda más bajo.
La solución que ofrece el mejor ancho de banda son las operaciones escalares con dos hilos.
El código que utilicé para comparar:
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <omp.h>
#define N 10000000
#define R 100
void mul(double *a, double *b) {
#pragma omp parallel for
for (int i = 0; i<N; i++) a[i] *= b[i];
}
int main() {
double maxbw = 2.4*2*8; // 2.4GHz * 2-channels * 64-bits * 1-byte/8-bits
double mem = 3*sizeof(double)*N*R*1E-9; // GB
double *a = (double*)malloc(sizeof *a * N);
double *b = (double*)malloc(sizeof *b * N);
//due to copy-on-write b must be initialized to get the correct bandwidth
//also, GCC will convert malloc + memset(0) to calloc so use memset(1)
memset(b, 1, sizeof *b * N);
double dtime = -omp_get_wtime();
for(int i=0; i<R; i++) mul(a,b);
dtime += omp_get_wtime();
printf("%.2f s, %.1f GB/s, %.1f%%/n", dtime, mem/dtime, 100*mem/dtime/maxbw);
free(a), free(b);
}
La vieja solución con el error de tiempo
La solución moderna para el ensamblaje en línea es usar elementos intrínsecos. Todavía hay casos en los que se necesita ensamblaje en línea, pero este no es uno de ellos.
Una solución intrínseca para su enfoque de ensamblaje en línea es simple:
void mul_SSE(double* a, double* b) {
for (int i = 0; i<N/2; i++)
_mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}
Déjame definir algunos códigos de prueba
#include <x86intrin.h>
#include <string.h>
#include <stdio.h>
#include <x86intrin.h>
#include <omp.h>
#define N 1000000
#define R 1000
typedef __attribute__(( aligned(32))) double aligned_double;
void (*fp)(aligned_double *a, aligned_double *b);
void mul(aligned_double* __restrict a, aligned_double* __restrict b) {
for (int i = 0; i<N; i++) a[i] *= b[i];
}
void mul_SSE(double* a, double* b) {
for (int i = 0; i<N/2; i++) _mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}
void mul_SSE_NT(double* a, double* b) {
for (int i = 0; i<N/2; i++) _mm_stream_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
}
void mul_SSE_OMP(double* a, double* b) {
#pragma omp parallel for
for (int i = 0; i<N; i++) a[i] *= b[i];
}
void test(aligned_double *a, aligned_double *b, const char *name) {
double dtime;
const double mem = 3*sizeof(double)*N*R/1024/1024/1024;
const double maxbw = 34.1;
dtime = -omp_get_wtime();
for(int i=0; i<R; i++) fp(a,b);
dtime += omp_get_wtime();
printf("%s /t time %.2f s, %.1f GB/s, efficency %.1f%%/n", name, dtime, mem/dtime, 100*mem/dtime/maxbw);
}
int main() {
double *a = (double*)_mm_malloc(sizeof *a * N, 32);
double *b = (double*)_mm_malloc(sizeof *b * N, 32);
//b must be initialized to get the correct bandwidth!!!
memset(a, 1, sizeof *a * N);
memset(b, 1, sizeof *a * N);
fp = mul, test(a,b, "mul ");
fp = mul_SSE, test(a,b, "mul_SSE ");
fp = mul_SSE_NT, test(a,b, "mul_SSE_NT ");
fp = mul_SSE_OMP, test(a,b, "mul_SSE_OMP");
_mm_free(a), _mm_free(b);
}
Ahora la primera prueba
g++ -O2 -fopenmp test.cpp
./a.out
mul time 1.67 s, 13.1 GB/s, efficiency 38.5%
mul_SSE time 1.00 s, 21.9 GB/s, efficiency 64.3%
mul_SSE_NT time 1.05 s, 20.9 GB/s, efficiency 61.4%
mul_SSE_OMP time 0.74 s, 29.7 GB/s, efficiency 87.0%
Entonces, con -O2 que no vectoriza los bucles, vemos que la versión SSE intrínseca es mucho más rápida que la solución simple C mul . efficiency = bandwith_measured/max_bandwidth donde max es 34.1 GB / s para mi sistema.
Segunda prueba
g++ -O3 -fopenmp test.cpp
./a.out
mul time 1.05 s, 20.9 GB/s, efficiency 61.2%
mul_SSE time 0.99 s, 22.3 GB/s, efficiency 65.3%
mul_SSE_NT time 1.01 s, 21.7 GB/s, efficiency 63.7%
mul_SSE_OMP time 0.68 s, 32.5 GB/s, efficiency 95.2%
Con -O3 vectoriza el bucle y la función intrínseca no ofrece esencialmente ninguna ventaja.
Tercera prueba
g++ -O3 -fopenmp -funroll-loops test.cpp
./a.out
mul time 0.85 s, 25.9 GB/s, efficency 76.1%
mul_SSE time 0.84 s, 26.2 GB/s, efficency 76.7%
mul_SSE_NT time 1.06 s, 20.8 GB/s, efficency 61.0%
mul_SSE_OMP time 0.76 s, 29.0 GB/s, efficency 85.0%
Con -funroll-loops GCC desenrolla los bucles ocho veces y vemos una mejora significativa a excepción de la solución de almacenamiento no temporal y no una ventaja real para la solución OpenMP.
Antes de desenrollar el lazo, el conjunto para mul wiht -O3 es
xor eax, eax
.L2:
movupd xmm0, XMMWORD PTR [rsi+rax]
mulpd xmm0, XMMWORD PTR [rdi+rax]
movaps XMMWORD PTR [rdi+rax], xmm0
add rax, 16
cmp rax, 8000000
jne .L2
rep ret
Con -O3 -funroll-loops el conjunto para mul es:
xor eax, eax
.L2:
movupd xmm0, XMMWORD PTR [rsi+rax]
movupd xmm1, XMMWORD PTR [rsi+16+rax]
mulpd xmm0, XMMWORD PTR [rdi+rax]
movupd xmm2, XMMWORD PTR [rsi+32+rax]
mulpd xmm1, XMMWORD PTR [rdi+16+rax]
movupd xmm3, XMMWORD PTR [rsi+48+rax]
mulpd xmm2, XMMWORD PTR [rdi+32+rax]
movupd xmm4, XMMWORD PTR [rsi+64+rax]
mulpd xmm3, XMMWORD PTR [rdi+48+rax]
movupd xmm5, XMMWORD PTR [rsi+80+rax]
mulpd xmm4, XMMWORD PTR [rdi+64+rax]
movupd xmm6, XMMWORD PTR [rsi+96+rax]
mulpd xmm5, XMMWORD PTR [rdi+80+rax]
movupd xmm7, XMMWORD PTR [rsi+112+rax]
mulpd xmm6, XMMWORD PTR [rdi+96+rax]
movaps XMMWORD PTR [rdi+rax], xmm0
mulpd xmm7, XMMWORD PTR [rdi+112+rax]
movaps XMMWORD PTR [rdi+16+rax], xmm1
movaps XMMWORD PTR [rdi+32+rax], xmm2
movaps XMMWORD PTR [rdi+48+rax], xmm3
movaps XMMWORD PTR [rdi+64+rax], xmm4
movaps XMMWORD PTR [rdi+80+rax], xmm5
movaps XMMWORD PTR [rdi+96+rax], xmm6
movaps XMMWORD PTR [rdi+112+rax], xmm7
sub rax, -128
cmp rax, 8000000
jne .L2
rep ret
Cuarta prueba
g++ -O3 -fopenmp -mavx test.cpp
./a.out
mul time 0.87 s, 25.3 GB/s, efficiency 74.3%
mul_SSE time 0.88 s, 24.9 GB/s, efficiency 73.0%
mul_SSE_NT time 1.07 s, 20.6 GB/s, efficiency 60.5%
mul_SSE_OMP time 0.76 s, 29.0 GB/s, efficiency 85.2%
Ahora la función no intrínseca es la más rápida (sin incluir la versión OpenMP).
Por lo tanto, no hay ninguna razón para utilizar intrínsecos o ensamblaje en línea en este caso porque podemos obtener el mejor rendimiento con las opciones de compilación adecuadas (por ejemplo, -funroll-loops , -funroll-loops , -mavx ).
Sistema de prueba: Ubuntu 16.10, Skylake ([email protected]), 32GB RAM. Ancho de banda máximo de memoria (34.1 GB / s) https://ark.intel.com/products/88967/Intel-Core-i7-6700HQ-Processor-6M-Cache-up-to-3_50-GHz
Aquí hay otra solución que vale la pena considerar. La instrucción cmp no es necesaria si contamos desde -N hasta cero y accedemos a las matrices como N+i . GCC debería haber arreglado esto hace mucho tiempo. Elimina una instrucción (aunque debido a la fusión macro-operativa el cmp y el jmp a menudo cuentan como una microoperación).
void mul_SSE_v2(double* a, double* b) {
for (ptrdiff_t i = -N; i<0; i+=2)
_mm_store_pd(&a[N + i], _mm_mul_pd(_mm_load_pd(&a[N + i]),_mm_load_pd(&b[N + i])));
Asamblea con -O3
mul_SSE_v2(double*, double*):
mov rax, -1000000
.L9:
movapd xmm0, XMMWORD PTR [rdi+8000000+rax*8]
mulpd xmm0, XMMWORD PTR [rsi+8000000+rax*8]
movaps XMMWORD PTR [rdi+8000000+rax*8], xmm0
add rax, 2
jne .L9
rep ret
}
Esta optimización solo será útil si las matrices encajan, por ejemplo, la memoria caché L1, es decir, no está leyendo desde la memoria principal.
Finalmente encontré una forma de obtener la solución C simple para no generar la instrucción cmp .
void mul_v2(aligned_double* __restrict a, aligned_double* __restrict b) {
for (int i = -N; i<0; i++) a[i] *= b[i];
}
Y luego llame a la función desde un archivo de objeto separado como este mul_v2(&a[N],&b[N]) así que esta es quizás la mejor solución. Sin embargo, si llama a la función desde el mismo archivo objeto (unidad de traducción) como el que está definido en el GCC genera la instrucción cmp nuevamente.
También,
void mul_v3(aligned_double* __restrict a, aligned_double* __restrict b) {
for (int i = -N; i<0; i++) a[N+i] *= b[N+i];
}
todavía genera la instrucción cmp y genera el mismo ensamblaje que la función mul .
La función mul_SSE_NT es tonta. Utiliza almacenes no temporales que solo son útiles cuando solo se escribe en la memoria, pero como la función lee y escribe en la misma dirección, los almacenes no temporales no solo son inútiles sino que dan resultados inferiores.
Las versiones anteriores de esta respuesta obtenían un ancho de banda incorrecto. La razón fue cuando las matrices no se inicializaron.