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

suma de matrices transpuestas

Question

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 de una matriz 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 <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;
}

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

• tamaño 512 - promedio 2.46 ms - http://ideone.com/1PV7m
• tamaño 513 - promedio de 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?

Answer 1

A modo de ilustración de la explicación de la respuesta de Luchian Grigore , a continuación se muestra cómo se ve la presencia del 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 de abajo significan lo siguiente:

• - no en caché,
• - en el caché,
• - golpe de caché,
• - acaba de leer de la memoria RAM,
• - señorita caché.

El caso 64x64:

Observe cómo casi todos los accesos a una nueva fila producen una falta de caché. Y ahora cómo busca el caso normal, una matriz de 65x65:

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

Answer 2

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

Para conocer los términos y la información detallada, consulte la entrada de la wiki sobre almacenamiento en caché . Lo limitaré 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 reflejar por el número de líneas nos da el tamaño del caché.

Para una dirección de memoria 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, ya que cada dirección de memoria es tan probable que se lea (dije idealmente ).

Está claro que pueden ocurrir solapamientos. En caso de 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 utiliza menos recientemente se sobrescribe con la memoria que se acaba de leer.

Intentaré seguir un poco 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 habrían estado ocupadas. La lectura de 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 separadas (para este ejemplo) 0x800 . Este es el paso crítico (nuevamente, para este ejemplo).

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

criticalStride = numberOfSets * lineSize

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

Esta es la parte teórica. A continuación, la explicación (también Agner, la estoy siguiendo de cerca para evitar cometer errores):

Suponga una matriz de 64x64 (recuerde, los efectos varían según el 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 e intercambiarlos con 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 guarde las líneas 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 se desalojará del caché. Cuando llegamos al final de la columna, todas las líneas de caché anteriores se habrían perdido y sería necesario volver a cargar en el acceso al siguiente elemento (se sobrescribe toda la línea).

Tener un tamaño que no sea un múltiplo de la zancada crítica confunde este escenario perfecto para el desastre, ya que ya no estamos tratando con elementos que son cruciales en la vertical, por lo que el número de recargas de caché se reduce considerablemente.

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

Answer 3

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 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 que es simplemente 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 especialmente portátil.

¿Podemos hacerlo mejor que eso? Sí, podemos: un enfoque general de este problema son los algoritmos de memoria caché ajenos que, como su nombre 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];
    }
}

Un poco más complejo, pero una breve prueba muestra algo bastante interesante en mi antiguo 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

Edit: Acerca de la influencia del tamaño: es mucho menos pronunciado, aunque aún se nota en cierta medida, porque estamos usando la solución iterativa como un nodo de hoja en lugar de recurrir a 1 (la optimización habitual para los algoritmos recursivos). Si establecemos LEAFSIZE = 1, el 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 100ms; este "punto de referencia" no es algo con lo que me sentirí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 colaboró ​​en esto ... supongo que sus trabajos son un buen punto de partida. Esos algoritmos aún son raramente descritos: CLR tiene una única 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 en 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];
            }
        }
    }
}