Multiplicación de matriz: pequeña diferencia en el tamaño de la matriz, gran diferencia en los tiempos

Tengo un código de matriz multiplicada que se ve así:

for(i = 0; i < dimension; i++) for(j = 0; j < dimension; j++) for(k = 0; k < dimension; k++) C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j]; 

Aquí, el tamaño de la matriz está representado por dimension . Ahora, si el tamaño de las matrices es 2000, lleva 147 segundos ejecutar este fragmento de código, mientras que si el tamaño de las matrices es 2048, tardará 447 segundos. Entonces, mientras la diferencia en no. de multiplicaciones es (2048 * 2048 * 2048) / (2000 * 2000 * 2000) = 1.073, la diferencia en los tiempos es 447/147 = 3. ¿Puede alguien explicar por qué sucede esto? Esperaba que escalara linealmente, lo que no sucede. No estoy tratando de hacer el código de multiplicación de matriz más rápido, simplemente tratando de entender por qué sucede.

Especificaciones: nodo AMD Opteron de doble núcleo (2.2GHz), 2G RAM, gcc v 4.5.0

Progtwig comstackdo como gcc -O3 simple.c

También he ejecutado esto en el comstackdor icc de Intel y he visto resultados similares.

EDITAR:

Como se sugirió en los comentarios / respuestas, ejecuté el código con dimensión = 2060 y tarda 145 segundos.

Aquí está el progtwig completo:

 #include  #include  #include  /* change dimension size as needed */ const int dimension = 2048; struct timeval tv; double timestamp() { double t; gettimeofday(&tv, NULL); t = tv.tv_sec + (tv.tv_usec/1000000.0); return t; } int main(int argc, char *argv[]) { int i, j, k; double *A, *B, *C, start, end; A = (double*)malloc(dimension*dimension*sizeof(double)); B = (double*)malloc(dimension*dimension*sizeof(double)); C = (double*)malloc(dimension*dimension*sizeof(double)); srand(292); for(i = 0; i < dimension; i++) for(j = 0; j < dimension; j++) { A[dimension*i+j] = (rand()/(RAND_MAX + 1.0)); B[dimension*i+j] = (rand()/(RAND_MAX + 1.0)); C[dimension*i+j] = 0.0; } start = timestamp(); for(i = 0; i < dimension; i++) for(j = 0; j < dimension; j++) for(k = 0; k < dimension; k++) C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j]; end = timestamp(); printf("\nsecs:%f\n", end-start); free(A); free(B); free(C); return 0; } 

Aquí está mi conjetura salvaje: caché

Podría ser que pueda caber 2 filas de 2000 double en el caché. Que es un poco menos que el caché L1 de 32kb. (mientras deja espacio otras cosas necesarias)

Pero cuando llegas a 2048, usa todo el caché (y dertwigs un poco porque necesitas espacio para otras cosas)

Suponiendo que la política de caché es LRU, el derrame de la memoria caché solo un poco hará que toda la fila se vacíe repetidamente y se vuelva a cargar en la memoria caché L1.

La otra posibilidad es la asociatividad del caché debido a la potencia de dos. Aunque creo que el procesador es asociativo L1 bidireccional, así que no creo que importe en este caso. (Pero voy a tirar la idea por ahí de todos modos)

Posible explicación 2: falta de memoria caché de conflicto debido a la superalineación en la memoria caché L2.

Su matriz B se está iterando en la columna. Entonces, el acceso es strided. Su tamaño total de datos es 2k x 2k que es aproximadamente 32 MB por matriz. Eso es mucho más grande que tu caché L2.

Cuando los datos no están alineados perfectamente, tendrá una localidad espacial decente en B. Aunque esté saltando filas y solo usando un elemento por línea de caché, la línea de caché permanece en la caché L2 para ser reutilizada en la siguiente iteración del ciclo medio.

Sin embargo, cuando los datos se alinean perfectamente (2048), estos saltos aterrizarán en la misma “forma de caché” y superarán con creces su asociatividad de caché L2. Por lo tanto, las líneas de caché accedidas de B no permanecerán en la memoria caché para la siguiente iteración. En cambio, necesitarán ser arrastrados desde Ram.

Definitivamente está obteniendo lo que llamo una resonancia de caché. Esto es similar al aliasing , pero no exactamente lo mismo. Dejame explicar.

Los cachés son estructuras de datos de hardware que extraen una parte de la dirección y la utilizan como un índice en una tabla, similar a una matriz en el software. (De hecho, los llamamos matrices en hardware). La matriz de caché contiene líneas de datos y tags de caché, a veces una entrada por índice en la matriz (asignada directamente), a veces varias (asociatividad de conjuntos de N vías). Se extrae una segunda parte de la dirección y se compara con la etiqueta almacenada en la matriz. Juntos, el índice y la etiqueta identifican de forma única una dirección de memoria de línea de caché. Finalmente, el rest de los bits de dirección identifica a qué bytes de la línea de caché se dirigen, junto con el tamaño del acceso.

Por lo general, el índice y la etiqueta son simples bitfields. Entonces una dirección de memoria se ve como

  ...Tag... | ...Index... | Offset_within_Cache_Line 

(A veces, el índice y la etiqueta son hash, por ejemplo, algunos XOR de otros bits en los bits de rango medio que son el índice. Mucho más raramente, a veces el índice y más raramente la etiqueta, son cosas como tomar un módulo de dirección de línea de caché número primo. Estos cálculos de índices más complicados son bashs de combatir el problema de la resonancia, que explico aquí. Todos sufren alguna forma de resonancia, pero los esquemas de extracción de bitfield más simples sufren resonancia en los patrones de acceso comunes, como usted ha encontrado.

Entonces, los valores típicos … hay muchos modelos diferentes de “Opteron Dual Core”, y no veo nada aquí que especifique cuál tiene. Escogiendo uno al azar, el manual más reciente que veo en el sitio web de AMD, Bios y Kernel Developer’s Guide (BKDG) para AMD Family 15h Models 00h-0Fh , 12 de marzo de 2012.

(Familia 15h = Familia Bulldozer, el procesador de gama alta más reciente – el BKDG menciona doble núcleo, aunque no conozco el número de producto que es exactamente lo que describes. Pero, de todos modos, la misma idea de resonancia se aplica a todos los procesadores, es solo que los parámetros como el tamaño de la caché y la asociatividad pueden variar un poco).

De p.33:

El procesador AMD Family 15h contiene una caché de datos L1 predicha de 4 vías de 16 kbytes con dos puertos de 128 bits. Esta es una memoria caché de escritura simultánea que admite hasta dos cargas de 128 bytes por ciclo. Está dividido en 16 bancos, cada uno de 16 bytes de ancho. […] Solo se puede realizar una carga desde un banco determinado de la memoria caché L1 en un solo ciclo.

Para resumir:

  • Línea de 64 bytes de caché => 6 bits de desplazamiento dentro de la línea de caché

  • 16KB / 4 vías => la resonancia es 4KB.

    Es decir, los bits de dirección 0-5 son el desplazamiento de la línea de caché.

  • 16KB / 64B líneas de caché => 2 ^ 14/2 ^ 6 = 2 ^ 8 = 256 líneas de caché en el caché.
    (Corrección de errores: Originalmente calculé incorrectamente esto como 128. que he arreglado todas las dependencias).

  • 4 vías asociativas => 256/4 = 64 índices en la matriz de caché. Yo (Intel) llamo a estos “conjuntos”.

    es decir, puede considerar que la memoria caché es una matriz de 32 entradas o conjuntos, cada entrada contiene 4 líneas de memoria caché y sus tags. (Es más complicado que esto, pero está bien).

(Por cierto, los términos “establecer” y “forma” tienen varias definiciones ).

  • hay 6 bits de índice, bits 6-11 en el esquema más simple.

    Esto significa que cualquier línea de caché que tenga exactamente los mismos valores en los bits de índice, bits 6-11, se correlacionará con el mismo conjunto de caché.

Ahora mira tu progtwig.

 C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j]; 

Loop k es el lazo más interno. El tipo base es doble, 8 bytes. Si dimensión = 2048, es decir 2K, entonces los elementos sucesivos de B[dimension*k+j] accedidos por el bucle serán 2048 * 8 = 16K bytes de separación. Todos se asignarán al mismo conjunto de la memoria caché L1; todos tendrán el mismo índice en la memoria caché. Lo que significa que, en lugar de tener 256 líneas de caché en el caché disponibles, solo habrá 4: la “asociatividad de 4 vías” del caché.

Es decir, es probable que obtenga un error de caché cada 4 iteraciones en este ciclo. No está bien.

(En realidad, las cosas son un poco más complicadas. Pero lo anterior es una buena primera comprensión. Las direcciones de las entradas de B mencionadas anteriormente son una dirección virtual. Por lo tanto, podría haber direcciones físicas ligeramente diferentes. Además, Bulldozer tiene un caché de predicción de manera, probablemente usando bits de direcciones virtuales para que no tenga que esperar una traducción de direcciones virtuales a físicas. Pero, en cualquier caso, su código tiene una “resonancia” de 16 K. La memoria caché de datos L1 tiene una resonancia de 16 K. No es bueno .)]

Si cambia la dimensión solo un poco, por ejemplo, a 2048 + 1, las direcciones de la matriz B se distribuirán en todos los conjuntos de la caché. Y obtendrá significativamente menos errores de caché.

Es una optimización bastante común para rellenar sus matrices, por ejemplo, para cambiar 2048 a 2049, para evitar este srt de resonancia. Pero “el locking de caché es una optimización aún más importante. http://suif.stanford.edu/papers/lam-asplos91.pdf


Además de la resonancia de línea de caché, hay otras cosas sucediendo aquí. Por ejemplo, la memoria caché L1 tiene 16 bancos, cada uno de 16 bytes de ancho. Con dimensión = 2048, los accesos B sucesivos en el bucle interno siempre irán al mismo banco. Entonces no pueden ir en paralelo, y si el acceso A pasa al mismo banco, perderás.

No creo, mirándolo, que esto sea tan grande como la resonancia del caché.

Y, sí, posiblemente, puede haber aliasing en marcha. Por ejemplo, los almacenamientos intermedios STLF (Almacenar para cargar) pueden estar comparando solo con un pequeño campo de bits y obteniendo coincidencias falsas.

(En realidad, si lo piensas, la resonancia en el caché es como un alias, relacionado con el uso de los campos de bits. La resonancia es causada por múltiples líneas de caché que mapean el mismo conjunto, sin extenderse. bits)


En general, mi recomendación para el ajuste:

  1. Pruebe el locking de caché sin ningún análisis adicional. Digo esto porque el locking de caché es fácil, y es muy probable que esto sea todo lo que necesita hacer.

  2. Después de eso, use VTune u OProf. O Cachegrind. O …

  3. Mejor aún, use una rutina de biblioteca bien ajustada para hacer la multiplicación matricial.

hay varias explicaciones posibles. Una posible explicación es lo que Mysticial sugiere: agotamiento de un recurso limitado (ya sea caché o TLB). Otra posibilidad probable es un locking de alias falso, que puede ocurrir cuando los accesos de memoria consecutivos están separados por un múltiplo de alguna potencia de dos (a menudo 4 KB).

Puede comenzar a reducir lo que está trabajando al trazar el tiempo / dimensión ^ 3 para un rango de valores. Si ha volado una memoria caché o un scope de TLB agotado, verá una sección más o menos plana seguida de un aumento brusco entre 2000 y 2048, seguido de otra sección plana. Si está viendo puestos relacionados con aliasing, verá un gráfico más o menos plano con un pico estrecho hacia arriba en 2048.

Por supuesto, esto tiene poder de diagnóstico, pero no es concluyente. Si quiere saber de manera concluyente cuál es el origen de la ralentización, le conviene aprender sobre los contadores de rendimiento , que pueden responder a este tipo de preguntas de manera definitiva.

Sé que esto es demasiado viejo, pero tomaré un bocado. Es (como se ha dicho) un problema de caché lo que causa la desaceleración en torno a potencias de dos. Pero hay otro problema con esto: es demasiado lento. Si miras tu ciclo de cálculo.

 for(i = 0; i < dimension; i++) for(j = 0; j < dimension; j++) for(k = 0; k < dimension; k++) C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j]; 

El bucle más interno cambia k por 1 en cada iteración, lo que significa que accedes solo a 1 doble del último elemento que usaste de A, pero una "dimensión" completa se duplica con respecto al último elemento de B. Esto no aprovecha ninguna ventaja de el almacenamiento en caché de los elementos de B.

Si cambias esto a:

 for(i = 0; i < dimension; i++) for(j = 0; j < dimension; j++) for(k = 0; k < dimension; k++) C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k]; 

Obtiene los mismos resultados exactos (errores de asociatividad de sum doble de módulo), pero es mucho más compatible con el caché ( local ). Lo intenté y me da mejoras sustanciales. Esto se puede resumir como

No multiplique matrices por definición, sino más bien, por filas


Ejemplo de aceleración (cambié tu código para tomar la dimensión como argumento)

 $ diff ac bc 42c42 < C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j]; --- > C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k]; $ make a cc ac -oa $ make b cc bc -ob $ ./a 1024 secs:88.732918 $ ./b 1024 secs:12.116630 

Como una ventaja (y lo que hace que esto esté relacionado con esta pregunta) es que este ciclo no sufre del problema anterior.

Si ya sabías todo esto, ¡me disculpo!

Un par de respuestas mencionaron problemas de caché L2.

Usted puede verificar esto con una simulación de caché. La herramienta cachegrind de Valgrind puede hacer eso.

 valgrind --tool=cachegrind --cache-sim=yes your_executable 

Establezca los parámetros de línea de comando para que coincidan con los parámetros L2 de su CPU.

Pruébalo con diferentes tamaños de matriz, probablemente verás un aumento repentino en la relación de falta de L2.