¿Por qué la transposición de una matriz de 512×512 es mucho más lenta que la transposición de una matriz de 513×513?

Después de realizar algunos experimentos en matrices cuadradas de diferentes tamaños, surgió un patrón. Invariablemente, la transposición de una matriz de tamaño 2^n es más lenta que la transposición de uno de tamaño 2^n+1 . Para valores pequeños de n , la diferencia no es mayor.

Sin embargo, se producen grandes diferencias sobre un valor de 512. (al menos para mí)

Descargo de responsabilidad: Sé que la función en realidad no transpone la matriz debido al doble intercambio de elementos, pero no hace ninguna diferencia.

Sigue el código:

 #define SAMPLES 1000 #define MATSIZE 512 #include  #include  int mat[MATSIZE][MATSIZE]; void transpose() { for ( int i = 0 ; i < MATSIZE ; i++ ) for ( int j = 0 ; j < MATSIZE ; j++ ) { int aux = mat[i][j]; mat[i][j] = mat[j][i]; mat[j][i] = aux; } } int main() { //initialize matrix for ( int i = 0 ; i < MATSIZE ; i++ ) for ( int j = 0 ; j < MATSIZE ; j++ ) mat[i][j] = i+j; int t = clock(); for ( int i = 0 ; i < SAMPLES ; i++ ) transpose(); int elapsed = clock() - t; std::cout << "Average for a matrix of " << MATSIZE << ": " << elapsed / SAMPLES; } 

Cambiar MATSIZE nos permite alterar el tamaño (¡duh!). Publiqué dos versiones en ideone:

  • tamaño 512 – promedio 2.46 ms – http://ideone.com/1PV7m
  • tamaño 513 – promedio 0.75 ms – http://ideone.com/NShpo

En mi entorno (MSVS 2010, optimizaciones completas), la diferencia es similar:

  • tamaño 512 – promedio 2,19 ms
  • tamaño 513 – promedio 0,57 ms

¿Por qué está pasando esto?

La explicación proviene de Agner Fog en Optimizing software in C ++ y se reduce a la forma en que se accede a los datos y se almacenan en la memoria caché.

Para conocer los términos y la información detallada, consulte la entrada del wiki sobre el almacenamiento en caché , voy a limitarlo aquí.

Un caché está organizado en conjuntos y líneas . En un momento, solo se usa un juego, del cual se puede usar cualquiera de las líneas que contiene. La memoria que una línea puede reflejar multiplicada por el número de líneas nos da el tamaño de la caché.

Para una dirección de memoria particular, podemos calcular en qué conjunto debe reflejarse con la fórmula:

 set = ( address / lineSize ) % numberOfsets 

Este tipo de fórmula ofrece una distribución idealmente uniforme en todos los conjuntos, ya que cada dirección de memoria tiene la misma probabilidad de leerse (lo dije de manera ideal ).

Está claro que las superposiciones pueden ocurrir. En caso de falta de caché, la memoria se lee en el caché y se reemplaza el valor anterior. Recuerde que cada conjunto tiene varias líneas, de las cuales la que se usó menos recientemente se sobrescribe con la memoria recién leída.

Trataré de seguir de alguna manera el ejemplo de Agner:

Supongamos que cada conjunto tiene 4 líneas, cada una con 64 bytes. Primero intentamos leer la dirección 0x2710 , que va en el conjunto 28 . Y luego también intentamos leer las direcciones 0x2F00 , 0x3700 , 0x3F00 y 0x4700 . Todos estos pertenecen al mismo conjunto. Antes de leer 0x4700 , todas las líneas del conjunto se habrían ocupado. Leer esa memoria desaloja una línea existente en el conjunto, la línea que inicialmente contenía 0x2710 . El problema radica en el hecho de que leemos las direcciones que están (para este ejemplo) separadas por 0x800 . Este es el paso crítico (nuevamente, para este ejemplo).

El paso crítico también se puede calcular:

 criticaStride = numberOfSets * lineSize 

Variables espaciadas criticalStride o una competencia múltiple aparte para las mismas líneas de caché.

Esta es la parte de la teoría. A continuación, la explicación (también Agner, la sigo de cerca para evitar cometer errores):

Suponga una matriz de 64×64 (recuerde, los efectos varían de acuerdo con la memoria caché) con un caché de 8kb, 4 líneas por conjunto * tamaño de línea de 64 bytes. Cada línea puede contener 8 de los elementos en la matriz ( int 64 bits).

El paso crítico sería 2048 bytes, que corresponden a 4 filas de la matriz (que es continua en la memoria).

Supongamos que estamos procesando la fila 28. Estamos intentando tomar los elementos de esta fila y cambiarlos por los elementos de la columna 28. Los primeros 8 elementos de la fila forman una línea de caché, pero entrarán en 8 diferentes líneas de caché en la columna 28. Recuerde, el paso crítico está a 4 filas (4 elementos consecutivos en una columna).

Cuando se alcanza el elemento 16 en la columna (4 líneas de caché por conjunto y 4 filas separadas = problema), el elemento ex-0 será desalojado de la memoria caché. Cuando lleguemos al final de la columna, todas las líneas de caché anteriores se habrían perdido y habría que volver a cargarlas al acceder al siguiente elemento (toda la línea se sobrescribe).

Tener un tamaño que no sea un múltiplo de la zancada crítica arruina este escenario perfecto para el desastre, ya que ya no estamos lidiando con elementos que son críticos en la vertical, por lo que la cantidad de recarga de caché se reduce drásticamente.

Otro descargo de responsabilidad : acabo de entender la explicación y espero haberlo descifrado, pero podría estar equivocado. De todos modos, estoy esperando una respuesta (o confirmación) de Mysticial . 🙂

Luchian da una explicación de por qué ocurre este comportamiento, pero pensé que sería una buena idea mostrar una posible solución a este problema y al mismo tiempo mostrar un poco acerca de los algoritmos de memoria caché.

Tu algoritmo básicamente hace:

 for (int i = 0; i < N; i++) for (int j = 0; j < N; j++) A[j][i] = A[i][j]; 

que es simplemente horrible para una CPU moderna. Una solución es conocer los detalles sobre su sistema de caché y modificar el algoritmo para evitar esos problemas. Funciona muy bien siempre que sepa esos detalles ... no especialmente portátil.

¿Podemos hacer algo mejor que eso? Sí, podemos: un enfoque general a este problema son los algoritmos de memoria caché que, como su nombre lo dice, evita depender de tamaños de caché específicos [1]

La solución sería así:

 void recursiveTranspose(int i0, int i1, int j0, int j1) { int di = i1 - i0, dj = j1 - j0; const int LEAFSIZE = 32; // well ok caching still affects this one here if (di >= dj && di > LEAFSIZE) { int im = (i0 + i1) / 2; recursiveTranspose(i0, im, j0, j1); recursiveTranspose(im, i1, j0, j1); } else if (dj > LEAFSIZE) { int jm = (j0 + j1) / 2; recursiveTranspose(i0, i1, j0, jm); recursiveTranspose(i0, i1, jm, j1); } else { for (int i = i0; i < i1; i++ ) for (int j = j0; j < j1; j++ ) mat[j][i] = mat[i][j]; } } 

Un poco más complejo, pero una breve prueba muestra algo bastante interesante en mi antigua e8400 con la versión VS2010 x64, testcode para MATSIZE 8192

 int main() { LARGE_INTEGER start, end, freq; QueryPerformanceFrequency(&freq); QueryPerformanceCounter(&start); recursiveTranspose(0, MATSIZE, 0, MATSIZE); QueryPerformanceCounter(&end); printf("recursive: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000)); QueryPerformanceCounter(&start); transpose(); QueryPerformanceCounter(&end); printf("iterative: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000)); return 0; } results: recursive: 480.58ms iterative: 3678.46ms 

Editar: Acerca de la influencia del tamaño: es mucho menos pronunciada, aunque se nota hasta cierto punto, porque usamos la solución iterativa como un nodo hoja en lugar de recurrir a 1 (la optimización habitual para algoritmos recursivos). Si establecemos LEAFSIZE = 1, la memoria caché no tiene ninguna influencia para mí [ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms - eso está dentro del margen de error, las fluctuaciones están en el área de 100 ms; este "punto de referencia" no es algo con lo que estaría demasiado cómodo si quisiéramos valores completamente precisos])

[1] Fuentes para esto: Bueno, si no puede obtener una conferencia de alguien que trabajó con Leiserson y colabore en esto ... supongo que sus documentos son un buen punto de partida. Esos algoritmos todavía se describen muy raramente: CLR tiene una sola nota al pie sobre ellos. Aún así es una gran manera de sorprender a la gente.


Editar (nota: no soy el que publicó esta respuesta; solo quería agregar esto):
Aquí hay una versión completa de C ++ del código anterior:

 template void transpose(InIt const input, OutIt const output, size_t const rows, size_t const columns, size_t const r1 = 0, size_t const c1 = 0, size_t r2 = ~(size_t) 0, size_t c2 = ~(size_t) 0, size_t const leaf = 0x20) { if (!~c2) { c2 = columns - c1; } if (!~r2) { r2 = rows - r1; } size_t const di = r2 - r1, dj = c2 - c1; if (di >= dj && di > leaf) { transpose(input, output, rows, columns, r1, c1, (r1 + r2) / 2, c2); transpose(input, output, rows, columns, (r1 + r2) / 2, c1, r2, c2); } else if (dj > leaf) { transpose(input, output, rows, columns, r1, c1, r2, (c1 + c2) / 2); transpose(input, output, rows, columns, r1, (c1 + c2) / 2, r2, c2); } else { for (ptrdiff_t i1 = (ptrdiff_t) r1, i2 = (ptrdiff_t) (i1 * columns); i1 < (ptrdiff_t) r2; ++i1, i2 += (ptrdiff_t) columns) { for (ptrdiff_t j1 = (ptrdiff_t) c1, j2 = (ptrdiff_t) (j1 * rows); j1 < (ptrdiff_t) c2; ++j1, j2 += (ptrdiff_t) rows) { output[j2 + i1] = input[i2 + j1]; } } } } 

Como ilustración de la explicación en la respuesta de Luchian Grigore , aquí se muestra la apariencia de la memoria caché de la matriz para los dos casos de matrices de 64×64 y 65×65 (consulte el enlace anterior para obtener detalles sobre los números).

Los colores en las animaciones a continuación significan lo siguiente:

  • blanco – No en el caché,
  • verde claro – en caché,
  • verde brillante – golpe de caché,
  • naranja – solo lea de RAM,
  • rojo – caché señorita.

El caso 64×64:

animación de presencia de caché para matriz 64x64

Observe cómo casi cada acceso a una nueva fila da como resultado una falta de caché. Y ahora cómo se ve el caso normal, una matriz de 65×65:

animación de presencia de caché para la matriz 65x65

Aquí puede ver que la mayoría de los accesos después del calentamiento inicial son éxitos de caché. Así es como la memoria caché de la CPU está diseñada para funcionar en general.