Why is transposing a matrix of <strong>512x512</strong> much slower than transposing a matrix of <strong>513x513</strong>? Pourquoi la transposition d'une matrice de <stro

Question

Après avoir mené quelques expériences sur des matrices carrées de tailles différentes, un schéma est apparu. Invariablement, transposer une matrice de taille 2^n est plus lent que de transposer une de taille 2^n+1. Pour de petites valeurs de n, la différence n'est pas majeure.

Cependant, de grandes différences surviennent au-dessus d'une valeur de 512. (du moins pour moi)

Avertissement : Je sais que la fonction ne transpose pas réellement la matrice en raison de l'échange double des éléments, mais cela ne change rien.

Voici le code :

#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()
{
   //initialisation de la matrice
   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 << "Moyenne pour une matrice de " << MATSIZE << ": " << elapsed / SAMPLES;
}

Changer MATSIZE nous permet de modifier la taille (évidemment!). 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 se produit-il?

Answer 1

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

Pour les termes et informations détaillées, consultez l' entrée wiki sur le caching, je vais simplifier ici.

Un cache est organisé en ensembles et lignes. À un moment donné, un seul ensemble est utilisé, parmi lesquels l'une 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 :

ensemble = ( adresse / tailleLigne ) % nombreEnsembles

Ce type de formule donne idéalement une distribution uniforme à travers les ensembles, car chaque adresse mémoire a autant de chances d'être lue (j'ai dit idéalement).

Il est clair que des chevauchements peuvent se produire. En cas de cache miss, la mémoire est lue dans le cache et la vieille valeur est remplacée. Rappelez-vous que chaque ensemble a un certain nombre de lignes, dont la moins récemment utilisée est écrasée par la mémoire nouvellement lue.

Je vais essayer de suivre en quelque sorte l'exemple d'Agner :

Supposons que chaque ensemble a 4 lignes, chacune contenant 64 octets. Nous tentons d'abord de lire l'adresse 0x2710, qui se trouve dans l'ensemble 28. Puis nous tentons également de lire les adresses 0x2F00, 0x3700, 0x3F00 et 0x4700. Toutes appartiennent au même ensemble. Avant de lire 0x4700, toutes les lignes de l'ensemble seront occupées. La lecture de cette mémoire évince une ligne existante dans l'ensemble, la ligne qui contenait initialement 0x2710. Le problème réside dans le fait que nous lisons des adresses qui sont (pour cet exemple) espacées de 0x800. C'est le pas critique (encore une fois, pour cet exemple).

Le pas critique peut également être calculé :

pasCritique = nombreEnsembles * tailleLigne

Les variables espacées de pasCritique ou d'un multiple concourent pour les mêmes lignes de cache.

C'est la partie théorique. Ensuite, l'explication (également d'Agner, je le suis de près pour éviter les erreurs) :

Supposons une matrice de 64x64 (rappelez-vous, les effets varient selon le cache) avec un cache de 8 ko, 4 lignes par ensemble * taille de ligne de 64 octets. Chaque ligne peut contenir 8 des éléments de la matrice (entier de 64 bits int).

Le pas critique serait de 2048 octets, correspondant à 4 rangées de la matrice (qui est continue en mémoire).

Supposons que nous traitions la rangée 28. Nous tentons de prendre les éléments de cette rangée et de les échanger avec les éléments de la colonne 28. Les 8 premiers éléments de la rangée forment une ligne de cache, mais ils iront dans 8 lignes de cache différentes dans la colonne 28. Rappelez-vous, le pas critique est de 4 rangées (4 éléments consécutifs dans une colonne).

Lorsque l'on atteint l'élément 16 dans la colonne (4 lignes de cache par ensemble & 4 rangées d'écart = problème), l'élément ex-0 sera évincé du cache. Lorsque nous atteignons la fin de la colonne, toutes les lignes de cache précédentes auront été perdues et devront être rechargées lors de l'accès à l'élément suivant (toute la ligne est écrasée).

Avoir une taille qui n'est pas un multiple du pas critique perturbe ce scénario parfait pour le désastre, car nous ne traitons plus d'éléments espacés de manière critique sur la verticale, le nombre de rechargements de cache est donc considérablement réduit.

Autre avertissement - Je viens de comprendre l'explication et j'espère l'avoir bien comprise, mais je pourrais me tromper. De toute façon, j'attends une réponse (ou une confirmation) de Mysticial. :)

Answer 2

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

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

• – pas dans le cache,
• – dans le cache,
• – cache hit,
• – juste lu depuis la RAM,
• – cache miss.

Le cas 64x64:

Remarquez comment presque chaque accès à une nouvelle ligne entraîne un cache miss. Et maintenant, voyez comment cela se passe pour le cas normal, une matrice 65x65:

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

---

_{Le code qui a généré les images pour les animations ci-dessus peut être consulté <a href="https://github.com/10110111/cpu-cache-visualization" rel="noreferrer">ici</a>.}

Answer 3

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 de montrer un peu sur les algorithmes sans 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];

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

Pouvons-nous faire mieux que cela ? Oui nous le pouvons: Une approche générale à ce problème est basée sur les algorithmes sans cache qui, comme le nom l'indique, évitent de dépendre des 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; // bon ok le caching affecte toujours celui-ci
    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 VS2010 x64 version finale, code de test pour MATSIZE 8192

int main() {
    LARGE_INTEGER start, end, freq;
    QueryPerformanceFrequency(&freq);
    QueryPerformanceCounter(&start);
    recursiveTranspose(0, MATSIZE, 0, MATSIZE);
    QueryPerformanceCounter(&end);
    printf("récursif: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));

    QueryPerformanceCounter(&start);
    transpose();
    QueryPerformanceCounter(&end);
    printf("itératif: %.2fms\n", (end.QuadPart - start.QuadPart) / (double(freq.QuadPart) / 1000));
    return 0;
}

résultats: 
récursif: 480.58ms
itératif: 3678.46ms

Édit: À propos de l'influence de la taille: Elle est beaucoup moins prononcée bien que toujours perceptible dans une certaine mesure, c'est parce que nous utilisons la solution itérative comme un nœud feuille au lieu de descendre en récursif jusqu'à 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 des 100ms; ce "test" n'est pas quelque chose avec lequel je serais trop à l'aise si nous voulions des valeurs complètement précises])

[1] Sources pour ces choses: Eh bien, si vous ne pouvez pas obtenir un cours magistral de quelqu'un qui a travaillé avec Leiserson et compagnie sur ce sujet.. je suppose que leurs documents sont un bon point de départ. Ces algorithmes sont encore assez rarement décrits - CLR a une note unique à leur sujet. Néanmoins, c'est un excellent moyen de surprendre les gens.

---

Édition (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
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];
            }
        }
    }
}