Pourquoi la transposition d'une matrice de 512x512 est-elle beaucoup plus lente que la transposition d'une matrice de 513x513?

Luchian donne une explication de pourquoi ce comportement se produit, mais j'ai pensé que ce serait une bonne idée de montrer une solution possible à ce problème et en même temps montrer un peu sur les algorithmes inconscients de cache.

votre algorithme fait essentiellement:

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

qui est tout simplement horrible pour un CPU moderne. Une solution consiste à connaître les détails de votre système de cache et à modifier l'algorithme pour éviter ces problèmes. Ça marche bien tant que tu connais ces détails.. pas spécialement portable.

Peut-on faire mieux que ça? Oui nous le pouvons: une approche générale de ce problème est algorithmes de cache inconscients qui, comme le nom le dit, évite de dépendre de tailles de cache spécifiques [1]

la solution ressemblerait à ceci:

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];
    }
}

légèrement plus complexe, mais un court test montre quelque chose de très intéressant sur mon ancien e8400 avec version VS2010 x64, code test pour 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: à propos de l'influence de la taille: elle est beaucoup moins prononcée bien qu'encore perceptible à un certain degré, c'est parce que nous utilisons la solution itérative comme un noeud de feuille au lieu de revenir à 1 (l'optimisation habituelle pour les algorithmes récursifs). Si nous définissons LEAFSIZE = 1, le cache n'a aucune influence pour moi [ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms - c'est à l'intérieur de la marge d'erreur, les fluctuations sont dans les 100ms zone; cette "référence" n'est pas quelque chose que je serais trop à l'aise avec, si nous voulions tout à fait précis valeurs])

[1] Sources pour ce genre de choses: eh Bien, si vous ne pouvez pas obtenir une conférence de quelqu'un qui a travaillé avec Leiserson et la coopération sur ce.. Je suppose que leurs papiers un bon point de départ. Ces algorithmes sont encore assez rarement décrits - CLR a une seule note de bas de page à leur sujet. C'est quand même une bonne façon de surprendre les gens.

Edit (note: Je ne suis pas celui qui a posté cette réponse; je voulais juste ajouter ceci):

Voici une version C++ complète du code ci-dessus:

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];
            }
        }
    }
}

comme illustration de l'explication dans réponse de Luchian Grigore , voici à quoi ressemble la présence de cache matricielle pour les deux cas de matrices 64x64 et 65x65 (voir le lien ci-dessus pour plus de détails sur les nombres).

les couleurs dans les animations ci-dessous signifient ce qui suit:

• – pas en cache,
• "151920920 – - en cache,
• – en cache,
• –
• – cache miss.

La 64x64 cas:

remarquez comment presque chaque l'accès à une nouvelle rangée entraîne une erreur de cache. Et maintenant comment il regarde pour le cas normal, une matrice 65x65:

ici vous pouvez voir que la plupart des accès après le réchauffement initial sont des hits de cache. C'est ainsi que le cache CPU est censé fonctionner en général.