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

Après avoir mené quelques expériences sur des matrices carrées de différentes tailles, un modèle est apparu. Invariablement, la transposition d'une matrice de taille 2^nest plus lente que la transposition d'une matrice de taille2^n+1 . Pour les petites valeurs de n, la différence n'est pas majeure.

De grandes différences se produisent cependant sur une valeur de 512. (au moins pour moi)

Avertissement: je sais que la fonction ne transpose pas réellement la matrice à cause du double échange d'éléments, mais cela ne fait aucune différence.

Suit le code:

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

Changer MATSIZEnous permet de modifier la taille (duh!). J'ai posté deux versions sur ideone:

• taille 512 - moyenne 2,46 ms - http://ideone.com/1PV7m
• taille 513 - moyenne 0,75 ms - http://ideone.com/NShpo

Dans mon environnement (MSVS 2010, optimisations complètes), la différence est similaire:

• taille 512 - moyenne 2,19 ms
• taille 513 - moyenne 0,57 ms

Pourquoi cela arrive-t-il?

L'explication vient d'Agner Fog in Optimizing software in C ++ et elle se réduit à la façon dont les données sont accessibles et stockées dans le cache.

Pour les termes et les informations détaillées, voir l' entrée wiki sur la mise en cache , je vais l'affiner ici.

Un cache est organisé en ensembles et en lignes . À la fois, un seul ensemble est utilisé, parmi lequel n'importe laquelle des lignes qu'il contient peut être utilisée. La mémoire qu'une ligne peut refléter multipliée par le nombre de lignes nous donne la taille du cache.

Pour une adresse mémoire particulière, nous pouvons calculer quel ensemble doit la refléter avec la formule:

set = ( address / lineSize ) % numberOfsets

Ce type de formule donne idéalement une distribution uniforme à travers les ensembles, car chaque adresse mémoire est aussi susceptible d'être lue (je l'ai dit idéalement ).

Il est clair que des chevauchements peuvent se produire. En cas de manque de cache, la mémoire est lue dans le cache et l'ancienne valeur est remplacée. N'oubliez pas que chaque ensemble a un certain nombre de lignes, dont la moins récemment utilisée est remplacée par la mémoire nouvellement lue.

Je vais essayer de suivre un peu l'exemple d'Agner:

Supposons que chaque ensemble comporte 4 lignes, chacune contenant 64 octets. Nous essayons d'abord de lire l'adresse 0x2710, qui va dans l'ensemble 28. Et puis nous essayons également de lire les adresses 0x2F00, 0x3700, 0x3F00et 0x4700. Tous ces éléments appartiennent au même ensemble. Avant la lecture 0x4700, toutes les lignes de l'ensemble auraient été occupées. La lecture de cette mémoire évite une ligne existante dans l'ensemble, la ligne qui tenait initialement 0x2710. Le problème réside dans le fait que nous lisons des adresses qui sont (pour cet exemple) 0x800séparées. Il s'agit de la foulée critique (encore une fois, pour cet exemple).

La foulée critique peut également être calculée:

criticalStride = numberOfSets * lineSize

Les variables espacées criticalStrideou multiples se disputent les mêmes lignes de cache.

Ceci est la partie théorique. Ensuite, l'explication (aussi Agner, je la suis de près pour éviter de faire des erreurs):

Supposons une matrice de 64x64 (rappelez-vous, les effets varient en fonction du cache) avec un cache de 8 Ko, 4 lignes par jeu * taille de ligne de 64 octets. Chaque ligne peut contenir 8 des éléments de la matrice (64 bits int).

La foulée critique serait de 2048 octets, ce qui correspond à 4 lignes de la matrice (qui est continue en mémoire).

Supposons que nous traitons la ligne 28. Nous essayons de prendre les éléments de cette ligne et de les échanger avec les éléments de la colonne 28. Les 8 premiers éléments de la ligne constituent une ligne de cache, mais ils iront dans 8 différents lignes de cache dans la colonne 28. N'oubliez pas que la foulée critique est distante de 4 lignes (4 éléments consécutifs dans une colonne).

Lorsque l'élément 16 est atteint dans la colonne (4 lignes de cache par ensemble et 4 lignes séparées = problème), l'élément ex-0 sera expulsé du cache. Lorsque nous atteignons la fin de la colonne, toutes les lignes de cache précédentes auraient été perdues et auraient dû être rechargées lors de l'accès à l'élément suivant (la ligne entière est remplacée).

Avoir une taille qui n'est pas un multiple de la foulée critique perturbe ce scénario parfait de catastrophe, car nous ne traitons plus d'éléments qui se démarquent sur la verticale, de sorte que le nombre de rechargements de cache est fortement réduit.

Un autre avertissement - je viens de comprendre ma explication et j'espère que je l'ai cloué, mais je peux me tromper. Quoi qu'il en soit, j'attends une réponse (ou une confirmation) de Mysticial . :)

Luchian explique 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 de montrer un peu les algorithmes de cache inconscients.

Votre algorithme fait essentiellement:

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

ce qui est tout simplement horrible pour un processeur moderne. Une solution consiste à connaître les détails de votre système de cache et à modifier l'algorithme pour éviter ces problèmes. Fonctionne très bien tant que vous connaissez ces détails .. pas spécialement portable.

Pouvons-nous faire mieux que ça? Oui, nous pouvons: Une approche générale de ce problème consiste à utiliser des algorithmes inconscients du cache qui, comme son nom l'indique, évitent 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];
    }
}

Un peu plus complexe, mais un court test montre quelque chose de très intéressant sur mon ancien e8400 avec la version VS2010 x64, testcode 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: A propos de l'influence de la taille: elle est beaucoup moins prononcée bien qu'elle soit encore perceptible dans une certaine mesure, c'est parce que nous utilisons la solution itérative comme nœud 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 la zone de 100 ms ; cette "référence" n'est pas quelque chose avec laquelle je serais trop à l'aise si nous voulions des valeurs complètement précises])

[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 ses collègues à ce sujet .. Je suppose que leurs papiers sont 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 un excellent moyen de surprendre les gens.

Edit (note: ce n'est pas moi qui ai posté cette réponse; je voulais juste l'ajouter):
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];
            }
        }
    }
}

Pour illustrer l'explication de la réponse de Luchian Grigore , voici à quoi ressemble la présence du cache de matrice 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 dans le cache,
• - en cache,
• - cache hit,
• - il suffit de lire à partir de la RAM,
• - cache miss.

Le boîtier 64x64:

Remarquez comment presque chaque accès à une nouvelle ligne entraîne un échec de cache. Et maintenant à quoi cela ressemble dans le cas normal, une matrice 65x65:

Ici, vous pouvez voir que la plupart des accès après l'échauffement initial sont des accès au cache. C'est ainsi que le cache CPU est censé fonctionner en général.

_{Le code qui a généré des images pour les animations ci-dessus peut être vu ici .}