punteros - suma de matrices c++

Medición del ancho de banda de memoria del producto de puntos de dos matrices

Question

El producto puntual de dos matrices.

for(int i=0; i<n; i++) {
sum += x[i]*y[i];
}

no reutiliza los datos, por lo que debe ser una operación vinculada a la memoria. Por lo tanto, debería poder medir el ancho de banda de memoria del producto punto.

Al usar el código en la why-vectorizing-the-loop-does-not-have-performance-improvement por la why-vectorizing-the-loop-does-not-have-performance-improvement obtengo un ancho de banda de 9.3 GB / s para mi sistema . Sin embargo, cuando intento calcular el ancho de banda utilizando el producto de puntos, obtengo más del doble de la velocidad para un solo hilo y más de tres veces la velocidad con varios subprocesos (mi sistema tiene cuatro núcleos / ocho hipervínculos). Esto no tiene sentido para mí, ya que una operación de enlace de memoria no debería beneficiarse de varios subprocesos. Aquí está la salida del código de abajo:

Xeon E5-1620, GCC 4.9.0, Linux kernel 3.13
dot 1 thread: 1.0 GB, sum 191054.81, time 4.98 s, 21.56 GB/s, 5.39 GFLOPS
dot_avx 1 thread 1.0 GB, sum 191043.33, time 5.16 s, 20.79 GB/s, 5.20 GFLOPS
dot_avx 2 threads: 1.0 GB, sum 191045.34, time 3.44 s, 31.24 GB/s, 7.81 GFLOPS
dot_avx 8 threads: 1.0 GB, sum 191043.34, time 3.26 s, 32.91 GB/s, 8.23 GFLOPS

¿Alguien puede explicarme por qué obtengo más del doble del ancho de banda para un hilo y más del triple del ancho de banda usando más de un hilo?

Aquí está el código que utilicé:

//g++ -O3 -fopenmp -mavx -ffast-math dot.cpp
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdint.h>
#include <x86intrin.h>
#include <omp.h>
extern "C" inline float horizontal_add(__m256 a) {
__m256 t1 = _mm256_hadd_ps(a,a);
__m256 t2 = _mm256_hadd_ps(t1,t1);
__m128 t3 = _mm256_extractf128_ps(t2,1);
__m128 t4 = _mm_add_ss(_mm256_castps256_ps128(t2),t3);
return _mm_cvtss_f32(t4);
}
extern "C" float dot_avx(float * __restrict x, float * __restrict y, const int n) {
x = (float*)__builtin_assume_aligned (x, 32);
y = (float*)__builtin_assume_aligned (y, 32);
float sum = 0;
#pragma omp parallel reduction(+:sum)
{
__m256 sum1 = _mm256_setzero_ps();
__m256 sum2 = _mm256_setzero_ps();
__m256 sum3 = _mm256_setzero_ps();
__m256 sum4 = _mm256_setzero_ps();
__m256 x8, y8;
#pragma omp for
for(int i=0; i<n; i+=32) {
x8 = _mm256_loadu_ps(&x[i]);
y8 = _mm256_loadu_ps(&y[i]);
sum1 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum1);
x8 = _mm256_loadu_ps(&x[i+8]);
y8 = _mm256_loadu_ps(&y[i+8]);
sum2 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum2);
x8 = _mm256_loadu_ps(&x[i+16]);
y8 = _mm256_loadu_ps(&y[i+16]);
sum3 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum3);
x8 = _mm256_loadu_ps(&x[i+24]);
y8 = _mm256_loadu_ps(&y[i+24]);
sum4 = _mm256_add_ps(_mm256_mul_ps(x8,y8),sum4);
}
sum += horizontal_add(_mm256_add_ps(_mm256_add_ps(sum1,sum2),_mm256_add_ps(sum3,sum4)));
}
return sum;
}
extern "C" float dot(float * __restrict x, float * __restrict y, const int n) {
x = (float*)__builtin_assume_aligned (x, 32);
y = (float*)__builtin_assume_aligned (y, 32);
float sum = 0;
for(int i=0; i<n; i++) {
sum += x[i]*y[i];
}
return sum;
}
int main(){
uint64_t LEN = 1 << 27;
float *x = (float*)_mm_malloc(sizeof(float)*LEN,64);
float *y = (float*)_mm_malloc(sizeof(float)*LEN,64);
for(uint64_t i=0; i<LEN; i++) { x[i] = 1.0*rand()/RAND_MAX - 0.5; y[i] = 1.0*rand()/RAND_MAX - 0.5;}
uint64_t size = 2*sizeof(float)*LEN;
volatile float sum = 0;
double dtime, rate, flops;
int repeat = 100;
dtime = omp_get_wtime();
for(int i=0; i<repeat; i++) sum += dot(x,y,LEN);
dtime = omp_get_wtime() - dtime;
rate = 1.0*repeat*size/dtime*1E-9;
flops = 2.0*repeat*LEN/dtime*1E-9;
printf("%f GB, sum %f, time %f s, %.2f GB/s, %.2f GFLOPS/n", 1.0*size/1024/1024/1024, sum, dtime, rate,flops);
sum = 0;
dtime = omp_get_wtime();
for(int i=0; i<repeat; i++) sum += dot_avx(x,y,LEN);
dtime = omp_get_wtime() - dtime;
rate = 1.0*repeat*size/dtime*1E-9;
flops = 2.0*repeat*LEN/dtime*1E-9;
printf("%f GB, sum %f, time %f s, %.2f GB/s, %.2f GFLOPS/n", 1.0*size/1024/1024/1024, sum, dtime, rate,flops);
}

Acabo de descargar, cumplir y ejecutar STREAM como lo sugirió Jonathan Dursi y aquí están los resultados:

Un hilo

Function Rate (MB/s) Avg time Min time Max time
Copy: 14292.1657 0.0023 0.0022 0.0023
Scale: 14286.0807 0.0023 0.0022 0.0023
Add: 14724.3906 0.0033 0.0033 0.0033
Triad: 15224.3339 0.0032 0.0032 0.0032

Ocho hilos

Function Rate (MB/s) Avg time Min time Max time
Copy: 24501.2282 0.0014 0.0013 0.0021
Scale: 23121.0556 0.0014 0.0014 0.0015
Add: 25263.7209 0.0024 0.0019 0.0056
Triad: 25817.7215 0.0020 0.0019 0.0027

Answer 1

Hay algunas cosas que suceden aquí, que se reducen a:

• Debe trabajar bastante duro para obtener hasta el último bit de rendimiento del subsistema de memoria; y
• Diferentes puntos de referencia miden diferentes cosas.

La primera ayuda a explicar por qué necesita varios subprocesos para saturar el ancho de banda de memoria disponible. Hay mucha concurrencia en el sistema de memoria, y aprovechando eso a menudo requerirá cierta concurrencia en su código de CPU. Una razón importante por la que la ayuda de múltiples subprocesos de ejecución es la ocultación de la latencia : mientras que un subproceso está bloqueado esperando que lleguen los datos, otro subproceso puede aprovechar algunos otros datos que están disponibles.

En este caso, el hardware lo ayuda mucho en un solo subproceso: como el acceso a la memoria es tan predecible, el hardware puede obtener los datos antes de que los necesite, lo que le brinda algunas de las ventajas de la ocultación de latencia incluso con un solo subproceso; pero hay límites a lo que puede hacer la captación previa. El prefetcher no se encargará de cruzar los límites de la página, por ejemplo. La referencia canónica para gran parte de esto es Lo que todo programador debería saber sobre la memoria de Ulrich Drepper , que ahora tiene la suficiente antigüedad como para que algunos espacios comiencen a mostrarse (la descripción general de chips de Intel de su procesador Sandy Bridge está here , en particular la integración más estrecha del hardware de gestión de memoria con la CPU).

En cuanto a la pregunta sobre la comparación con memset, mbw o STREAM , la comparación entre los puntos de referencia siempre causará dolores de cabeza, incluso los puntos de referencia que afirman que miden lo mismo. En particular, el "ancho de banda de la memoria" no es un solo número, el rendimiento varía bastante dependiendo de las operaciones. Tanto mbw como Stream realizan alguna versión de una operación de copia, con las operaciones de STREAM que se detallan aquí (tomadas directamente de la página web, todos los operandos son puntos flotantes de doble precisión):

------------------------------------------------------------------
name kernel bytes/iter FLOPS/iter
------------------------------------------------------------------
COPY: a(i) = b(i) 16 0
SCALE: a(i) = q*b(i) 16 1
SUM: a(i) = b(i) + c(i) 24 1
TRIAD: a(i) = b(i) + q*c(i) 24 2
------------------------------------------------------------------

así que aproximadamente 1 / 2-1 / 3 de las operaciones de memoria en estos casos se escriben (y todo es una escritura en el caso de memset). Si bien las escrituras individuales pueden ser un poco más lentas que las lecturas, el problema más grande es que es mucho más difícil saturar el subsistema de memoria con escrituras porque, por supuesto, no puede hacer el equivalente de realizar una escritura previa. La intercalación de las lecturas y escrituras ayuda, pero su ejemplo de producto de puntos, que es esencialmente todas las lecturas, será el mejor caso posible para fijar la aguja en el ancho de banda de la memoria.

Además, el índice de referencia STREAM está (intencionalmente) escrito de manera completamente portátil, con solo algunos pragmas del compilador para sugerir la vectorización, por lo que superar el índice de referencia STREAM no es necesariamente una señal de advertencia, especialmente cuando lo que estás haciendo son dos lecturas de transmisión.

Answer 2

Hice mi propio código de referencia de memoria github.com/zboson/bandwidth

Aquí están los resultados actuales para ocho hilos:

write: 0.5 GB, time 2.96e-01 s, 18.11 GB/s
copy: 1 GB, time 4.50e-01 s, 23.85 GB/s
scale: 1 GB, time 4.50e-01 s, 23.85 GB/s
add: 1.5 GB, time 6.59e-01 s, 24.45 GB/s
mul: 1.5 GB, time 6.56e-01 s, 24.57 GB/s
triad: 1.5 GB, time 6.61e-01 s, 24.37 GB/s
vsum: 0.5 GB, time 1.49e-01 s, 36.09 GB/s, sum -8.986818e+03
vmul: 0.5 GB, time 9.00e-05 s, 59635.10 GB/s, sum 0.000000e+00
vmul_sum: 1 GB, time 3.25e-01 s, 33.06 GB/s, sum 1.910421e+04

Aquí están los resultados de corrientes para 1 hilo:

write: 0.5 GB, time 4.65e-01 s, 11.54 GB/s
copy: 1 GB, time 7.51e-01 s, 14.30 GB/s
scale: 1 GB, time 7.45e-01 s, 14.41 GB/s
add: 1.5 GB, time 1.02e+00 s, 15.80 GB/s
mul: 1.5 GB, time 1.07e+00 s, 15.08 GB/s
triad: 1.5 GB, time 1.02e+00 s, 15.76 GB/s
vsum: 0.5 GB, time 2.78e-01 s, 19.29 GB/s, sum -8.990941e+03
vmul: 0.5 GB, time 1.15e-05 s, 468719.08 GB/s, sum 0.000000e+00
vmul_sum: 1 GB, time 5.72e-01 s, 18.78 GB/s, sum 1.910549e+04

1. escribir: escribe una constante (3.14159) en una matriz. Esto debería ser como memset .
2. Copiar, escalar, agregar y tríada se definen de la misma manera que en STREAM
3. mul: a(i) = b(i) * c(i)
4. vsum: sum += a(i)
5. vmul: sum *= a(i)
6. vmul_sum: sum += a(i)*b(i) // el producto punto

Mis resultados son consistentes con STREAM. Tengo el mayor ancho de banda para vsum . El método vmul no funciona actualmente (una vez que el valor es cero, termina antes). Puedo obtener resultados ligeramente mejores (alrededor del 10%) usando intrínsecos y desenrollando el bucle que agregaré más adelante.