Warum ist das Transponieren einer Matrix von 512 x 512 viel langsamer als das Transponieren einer Matrix von 513 x 513?

Nach einigen Experimenten mit quadratischen Matrizen unterschiedlicher Größe wurde ein Muster erstellt. Das Transponieren einer Größenmatrix 2^nist immer langsamer als das Transponieren einer Größenmatrix2^n+1 . Für kleine Werte vonn ist der Unterschied nicht groß.

Große Unterschiede treten jedoch bei einem Wert von 512 auf. (Zumindest für mich)

Haftungsausschluss: Ich weiß, dass die Funktion die Matrix aufgrund des doppelten Austauschs von Elementen nicht wirklich transponiert, aber es macht keinen Unterschied.

Folgt dem 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;
}

Durch Ändern MATSIZEkönnen wir die Größe ändern (duh!). Ich habe zwei Versionen auf ideone gepostet:

• Größe 512 - Durchschnitt 2,46 ms - http://ideone.com/1PV7m
• Größe 513 - Durchschnitt 0,75 ms - http://ideone.com/NShpo

In meiner Umgebung (MSVS 2010, vollständige Optimierungen) ist der Unterschied ähnlich:

• Größe 512 - Durchschnitt 2,19 ms
• Größe 513 - Durchschnitt 0,57 ms

Warum passiert dies?

Die Erklärung stammt von Agner Fog in Optimizing Software in C ++ und reduziert sich darauf, wie auf Daten zugegriffen und im Cache gespeichert wird.

Begriffe und detaillierte Informationen finden Sie im Wiki-Eintrag zum Caching . Ich werde ihn hier eingrenzen.

Ein Cache ist in Mengen und Zeilen organisiert . Es wird jeweils nur ein Satz verwendet, von dem jede der darin enthaltenen Zeilen verwendet werden kann. Der Speicher, den eine Zeile mal die Anzahl der Zeilen spiegeln kann, gibt uns die Cache-Größe.

Für eine bestimmte Speicheradresse können wir mit der Formel berechnen, welche Menge sie spiegeln soll:

set = ( address / lineSize ) % numberOfsets

Diese Art von Formel ergibt idealerweise eine gleichmäßige Verteilung über die Sätze, da jede Speicheradresse genauso wahrscheinlich gelesen wird (ich sagte idealerweise ).

Es ist klar, dass Überlappungen auftreten können. Bei einem Cache-Fehler wird der Speicher im Cache gelesen und der alte Wert ersetzt. Denken Sie daran, dass jeder Satz eine Anzahl von Zeilen enthält, von denen die zuletzt verwendete mit dem neu gelesenen Speicher überschrieben wird.

Ich werde versuchen, dem Beispiel von Agner etwas zu folgen:

Angenommen, jeder Satz hat 4 Zeilen mit jeweils 64 Bytes. Wir versuchen zuerst, die Adresse zu lesen 0x2710, die im Set enthalten ist 28. Und dann versuchen wir auch Adressen zu lesen 0x2F00, 0x3700, 0x3F00und 0x4700. Alle diese gehören zum selben Set. Vor dem Lesen 0x4700wären alle Zeilen im Set belegt gewesen. Durch das Lesen dieses Speichers wird eine vorhandene Zeile im Satz entfernt, die ursprünglich gehalten wurde 0x2710. Das Problem liegt in der Tatsache, dass wir Adressen lesen, die (für dieses Beispiel) 0x800voneinander entfernt sind. Dies ist der entscheidende Schritt (auch in diesem Beispiel).

Der kritische Schritt kann auch berechnet werden:

criticalStride = numberOfSets * lineSize

Variablen mit Abstand criticalStrideoder mehreren Abständen konkurrieren um dieselben Cache-Zeilen.

Dies ist der theoretische Teil. Als nächstes die Erklärung (auch Agner, ich verfolge sie genau, um Fehler zu vermeiden):

Angenommen, eine Matrix von 64 x 64 (denken Sie daran, die Effekte variieren je nach Cache) mit einem 8-KB-Cache, 4 Zeilen pro Satz * Zeilengröße von 64 Byte. Jede Zeile kann 8 der Elemente in der Matrix enthalten (64-Bit int).

Der kritische Schritt wäre 2048 Bytes, was 4 Zeilen der Matrix entspricht (die im Speicher kontinuierlich ist).

Angenommen, wir verarbeiten Zeile 28. Wir versuchen, die Elemente dieser Zeile mit den Elementen aus Spalte 28 zu tauschen. Die ersten 8 Elemente der Zeile bilden eine Cache-Zeile, aber sie werden in 8 verschiedene Elemente unterteilt Cache-Zeilen in Spalte 28. Denken Sie daran, dass der kritische Schritt 4 Zeilen voneinander entfernt ist (4 aufeinanderfolgende Elemente in einer Spalte).

Wenn Element 16 in der Spalte erreicht ist (4 Cache-Zeilen pro Satz und 4 Zeilen voneinander entfernt = Problem), wird das Ex-0-Element aus dem Cache entfernt. Wenn wir das Ende der Spalte erreicht haben, wären alle vorherigen Cache-Zeilen verloren gegangen und müssten beim Zugriff auf das nächste Element neu geladen werden (die gesamte Zeile wird überschrieben).

Eine Größe, die nicht ein Vielfaches des kritischen Schrittes ist, bringt dieses perfekte Szenario für eine Katastrophe durcheinander , da es sich nicht mehr um Elemente handelt, die in der Vertikalen einen kritischen Schritt voneinander entfernt sind, sodass die Anzahl der Cache-Neuladungen erheblich reduziert wird.

Noch ein Haftungsausschluss - ich habe mich gerade mit der Erklärung beschäftigt und hoffe, dass ich sie verstanden habe, aber ich könnte mich irren. Wie auch immer, ich warte auf eine Antwort (oder Bestätigung) von Mysticial . :) :)

Luchian gibt eine Erklärung dafür, warum dieses Verhalten auftritt, aber ich dachte, es wäre eine gute Idee, eine mögliche Lösung für dieses Problem aufzuzeigen und gleichzeitig ein wenig über Algorithmen zu zeigen, die den Cache nicht kennen.

Ihr Algorithmus macht im Grunde:

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

Das ist einfach schrecklich für eine moderne CPU. Eine Lösung besteht darin, die Details Ihres Cache-Systems zu kennen und den Algorithmus zu optimieren, um diese Probleme zu vermeiden. Funktioniert hervorragend, solange Sie diese Details kennen. Nicht besonders tragbar.

Können wir es besser machen? Ja, wir können: Ein allgemeiner Ansatz für dieses Problem sind Algorithmen , die den Cache nicht kennen und die , wie der Name schon sagt, vermeiden, von bestimmten Cache-Größen abhängig zu sein [1].

Die Lösung würde so aussehen:

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

Etwas komplexer, aber ein kurzer Test zeigt etwas ziemlich Interessantes auf meinem alten e8400 mit VS2010 x64 Release, Testcode für 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

Bearbeiten: Über den Einfluss der Größe: Es ist viel weniger ausgeprägt, obwohl es bis zu einem gewissen Grad immer noch spürbar ist. Dies liegt daran, dass wir die iterative Lösung als Blattknoten verwenden, anstatt auf 1 zu rekursieren (die übliche Optimierung für rekursive Algorithmen). Wenn wir LEAFSIZE = 1 setzen, hat der Cache für mich keinen Einfluss [ 8193: 1214.06; 8192: 1171.62ms, 8191: 1351.07ms- das liegt innerhalb der Fehlergrenze, die Schwankungen liegen im Bereich von 100 ms; Dieser "Benchmark" ist nichts, mit dem ich mich zu wohl fühlen würde, wenn wir völlig genaue Werte wollten])

[1] Quellen für dieses Zeug: Nun, wenn Sie keinen Vortrag von jemandem bekommen können, der mit Leiserson und Co. daran gearbeitet hat. Ich gehe davon aus, dass ihre Arbeiten ein guter Ausgangspunkt sind. Diese Algorithmen werden noch recht selten beschrieben - CLR enthält eine einzige Fußnote. Trotzdem ist es eine großartige Möglichkeit, Menschen zu überraschen.

Bearbeiten (Hinweis: Ich bin nicht derjenige, der diese Antwort gepostet hat; ich wollte dies nur hinzufügen):
Hier ist eine vollständige C ++ - Version des obigen Codes:

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

Zur Veranschaulichung der Erklärung in Luchian Grigores Antwort sehen Sie hier, wie die Matrix-Cache-Präsenz für die beiden Fälle von 64x64- und 65x65-Matrizen aussieht (Einzelheiten zu Zahlen finden Sie oben).

Die Farben in den folgenden Animationen bedeuten Folgendes:

• - nicht im Cache,
• - im Cache,
• - Cache getroffen,
• - einfach aus dem RAM lesen,
• - Cache-Miss.

Der 64x64-Fall:

Beachten Sie, dass fast jeder Zugriff auf eine neue Zeile zu einem Cache-Fehler führt. Und jetzt, wie es für den Normalfall aussieht, eine 65x65-Matrix:

Hier können Sie sehen, dass die meisten Zugriffe nach dem ersten Aufwärmen Cache-Treffer sind. So soll der CPU-Cache im Allgemeinen funktionieren.

_{Der Code, der Frames für die obigen Animationen generiert hat, ist hier zu sehen .}