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

218

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^nes más lenta que la transposición de una de tamaño2^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 <time.h>
#include <iostream>
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;
}

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

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?

Luchian Grigore
fuente
99
Tu código me parece poco amigable.
CodesInChaos
77
Es casi el mismo problema que esta pregunta: stackoverflow.com/questions/7905760/…
Mysticial
¿Te gustaría elavorar, @CodesInChaos? (O cualquier otra persona)
corazza
@Bane ¿Qué tal leer la respuesta aceptada?
CodesInChaos
44
@nzomkxia No tiene sentido medir nada sin optimizaciones. Con las optimizaciones deshabilitadas, el código generado estará lleno de basura extraña que ocultará otros cuellos de botella. (como la memoria)
Mysticial

Respuestas:

197

La explicación proviene de Agner Fog en Optimización de software en C ++ y se reduce a cómo se accede a los datos y se almacenan en el caché.

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

Un caché se organiza en conjuntos y líneas . A la vez, solo se usa un conjunto, del cual se puede usar cualquiera de las líneas que contiene. La memoria que una línea puede duplicar por la cantidad de líneas nos da el tamaño de caché.

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

set = ( address / lineSize ) % numberOfsets

Este tipo de fórmula idealmente proporciona una distribución uniforme entre los conjuntos, porque es probable que cada dirección de memoria sea leída (dije idealmente ).

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

Intentaré seguir el ejemplo de Agner:

Suponga que cada conjunto tiene 4 líneas, cada una con 64 bytes. Primero intentamos leer la dirección 0x2710, que va en conjunto 28. Y entonces también pretende hacer una lectura direcciones 0x2F00, 0x3700, 0x3F00y 0x4700. Todos estos pertenecen al mismo conjunto. Antes de leer 0x4700, todas las líneas en el conjunto habrían sido ocupadas. Leer esa memoria desaloja una línea existente en el conjunto, la línea que inicialmente estaba retenida 0x2710. El problema radica en el hecho de que leemos direcciones que están (para este ejemplo) 0x800separadas. Este es el paso crítico (nuevamente, para este ejemplo).

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

criticalStride = numberOfSets * lineSize

Las variables espaciadas criticalStrideo una separación múltiple compiten por 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 64x64 (recuerde, los efectos varían según el caché) con un caché de 8 kb, 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 (64 bits int).

La zancada crítica sería de 2048 bytes, que corresponden a 4 filas de la matriz (que es continua en la memoria).

Supongamos que estamos procesando la fila 28. Intentamos tomar los elementos de esta fila e intercambiarlos con los elementos de la columna 28. Los primeros 8 elementos de la fila forman una línea de caché, pero irán a 8 diferentes líneas de caché en la columna 28. Recuerde, el paso crítico está separado por 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á expulsado de la caché. Cuando llegamos al final de la columna, todas las líneas de caché anteriores se habrían perdido y habría que recargar al acceder al siguiente elemento (se sobrescribe toda la línea).

Tener un tamaño que no es 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 el número de recargas de caché se reduce drásticamente.

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

Luchian Grigore
fuente
Ah, y la próxima vez Solo hazme un ping directamente a través del Lounge . No encuentro cada instancia de nombre en SO. :) Solo vi esto a través de las notificaciones periódicas por correo electrónico.
Mysticial
@Mysticial @Luchian Grigore Uno de mis amigos me dice que su Intel core i3PC en ejecución Ubuntu 11.04 i386muestra casi el mismo rendimiento con gcc 4.6 . Y lo mismo ocurreIntel Core 2 Duo con mi computadora con mingw gcc4.4 , que está funcionando windows 7(32). Muestra una gran diferencia cuando Compilé este segmento con una PC un poco más antigua intel centrinocon gcc 4.6 , que se está ejecutando ubuntu 12.04 i386.
Hongxu Chen el
También tenga en cuenta que el acceso a la memoria donde las direcciones difieren en un múltiplo de 4096 tiene una dependencia falsa de las CPU de la familia Intel SnB. (es decir, el mismo desplazamiento dentro de una página). Esto puede reducir el rendimiento cuando algunas de las operaciones son tiendas, especialmente. Una mezcla de cargas y tiendas.
Peter Cordes
which goes in set 24¿quiso decir "en el set 28 " en su lugar? ¿Y asumes 32 sets?
Ruslan
Tienes razón, es 28. :) También verifiqué dos veces el documento vinculado, para la explicación original puedes navegar a 9.2 Organización de caché
Luchian Grigore
78

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 sobre los algoritmos ajenos a la memoria caché.

Su algoritmo básicamente hace:

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

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

¿Podemos hacerlo mejor que eso? Sí, podemos: Un enfoque general para este problema son los algoritmos ajenos a la caché que, como su nombre lo indica, evitan depender de tamaños de caché específicos [1]

La solución se verí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];
    }
}

Ligeramente más complejo, pero una breve prueba muestra algo bastante interesante en mi antiguo e8400 con la versión VS2010 x64, código de prueba 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 pronunciado, aunque todavía se nota en cierto grado, porque estamos usando 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, el caché no tiene influencia para mí [ 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 me sentiría demasiado cómodo si quisiéramos valores completamente precisos])

[1] Fuentes de estas cosas: Bueno, si no puede obtener una conferencia de alguien que trabajó con Leiserson y compañía 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 de página sobre ellos. Aún así, es una excelente manera de sorprender a la gente.


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

template<class InIt, class OutIt>
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];
            }
        }
    }
}
Voo
fuente
2
Esto sería relevante si comparara los tiempos entre matrices de diferentes tamaños, no recursivas e iterativas. Pruebe la solución recursiva en una matriz de los tamaños especificados.
Luchian Grigore
@Luchian Como ya explicaste por qué está viendo el comportamiento, pensé que era bastante interesante presentar una solución a este problema en general.
Voo
Porque, me pregunto por qué una matriz más grande tarda menos tiempo en procesarse, no busca un algoritmo más rápido ...
Luchian Grigore
@Luchian Las diferencias entre 16383 y 16384 son ... 28 vs 27ms para mí aquí, o alrededor de 3.5% - no es realmente significativo. Y me sorprendería si lo fuera.
Voo
3
Puede ser interesante explicar qué recursiveTransposehace, es decir, que no llena tanto el caché al operar en pequeños mosaicos (de LEAFSIZE x LEAFSIZEdimensión).
Matthieu M.
60

Como ilustración de la explicación en la respuesta de Luchian Grigore , así es como se ve la presencia de la memoria caché de matriz para los dos casos de matrices de 64x64 y 65x65 (consulte el enlace de arriba para obtener detalles sobre los números).

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

  • blanco - no en caché,
  • verde claro - en caché,
  • verde brillante - golpe de caché,
  • naranja - Acabo de leer de RAM,
  • rojo - caché señorita.

El caso de 64x64:

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

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

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

Aquí puede ver que la mayoría de los accesos después del calentamiento inicial son aciertos de caché. Así es como se pretende que el caché de la CPU funcione en general.


El código que generó marcos para las animaciones anteriores se puede ver aquí .

Ruslan
fuente
¿Por qué los hits de caché de escaneo vertical no se guardan en el primer caso, pero sí en el segundo caso? Parece que se accede a un bloque determinado exactamente una vez para la mayoría de los bloques en ambos ejemplos.
Josiah Yoder
De la respuesta de @ LuchianGrigore puedo ver que todas las líneas de la columna pertenecen al mismo conjunto.
Josiah Yoder
Sí, gran ilustración. Veo que están a la misma velocidad. Pero en realidad, no lo son, ¿no?
kelalaka
@kelalaka sí, la animación FPS es la misma. No simulé la desaceleración, solo los colores son importantes aquí.
Ruslan
Sería interesante tener dos imágenes estáticas que ilustren los diferentes conjuntos de caché.
Josiah Yoder