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

Nachdem einige Experimente mit quadratischen Matrizen unterschiedlicher Größe durchgeführt wurden, kam ein Muster zum Vorschein. Ausnahmslos, Transponieren einer Größenmatrix 2^n ist langsamer als das Transponieren einer Größe 2^n+1. Für kleine Werte von nder Unterschied ist nicht groß.

Große Unterschiede treten jedoch ab einem Wert von 512 auf. (zumindest bei mir)

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

Ändern MATSIZE lassen Sie uns 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 kommt von Agner Fog in Optimierung von Software in C++ und es reduziert sich darauf, wie auf Daten zugegriffen und im Cache gespeichert wird.

Bedingungen und detaillierte Informationen finden Sie unter Wiki-Eintrag zum Cachingich werde es hier eingrenzen.

Ein Cache ist in organisiert setzt und Linien. Es wird jeweils nur ein Satz verwendet, aus dem alle darin enthaltenen Zeilen verwendet werden können. Der Speicher, den eine Zeile spiegeln kann, multipliziert mit der Anzahl der Zeilen ergibt die Cache-Größe.

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

set = ( address / lineSize ) % numberOfsets

Diese Art von Formel ergibt im Idealfall eine gleichmäßige Verteilung über die Sätze, da jede Speicheradresse mit gleicher Wahrscheinlichkeit gelesen wird (ich sagte im Idealfall).

Es ist klar, dass es zu Überschneidungen kommen kann. Bei einem Cache-Miss wird der Speicher im Cache gelesen und der alte Wert ersetzt. Denken Sie daran, dass jeder Satz eine Reihe von Zeilen hat, von denen die am längsten 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 0x2710die in den Satz geht 28. Und dann versuchen wir auch, Adressen zu lesen 0x2F00, 0x3700, 0x3F00 und 0x4700. All dies gehört zu demselben Satz. Vor dem Lesen 0x4700, wären alle Leitungen im Set besetzt gewesen. Das Lesen dieses Speichers entfernt eine vorhandene Zeile in der Menge, die Zeile, die ursprünglich gehalten wurde 0x2710. Das Problem liegt darin, dass wir Adressen lesen, die (für dieses Beispiel) 0x800 ein Teil. Dies ist das kritischer Schritt (wieder für dieses Beispiel).

Die kritische Schrittlänge kann auch berechnet werden:

criticalStride = numberOfSets * lineSize

Variablen beabstandet criticalStride oder mehrere auseinander konkurrieren um die gleichen Cache-Zeilen.

Dies ist der Theorieteil. Als nächstes die Erklärung (auch Agner, ich folge ihr genau, um Fehler zu vermeiden):

Nehmen Sie eine Matrix von 64 x 64 an (denken Sie daran, dass die Auswirkungen je nach Cache variieren) mit einem 8-KB-Cache, 4 Zeilen pro Satz * Zeilengröße von 64 Bytes. 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 zu nehmen und sie mit den Elementen aus Spalte 28 auszutauschen. Die ersten 8 Elemente der Zeile bilden eine Cache-Zeile, aber sie gehen in 8 verschiedene 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 wird (4 Cache-Zeilen pro Satz und 4 Zeilen auseinander = Problem), wird das ex-0-Element aus dem Cache entfernt. Wenn wir das Ende der Spalte erreichen, wären alle vorherigen Cache-Zeilen verloren gegangen und müssten beim Zugriff auf das nächste Element neu geladen werden (die ganze Zeile wird überschrieben).

Eine Größe zu haben, die kein Vielfaches der kritischen Schrittlänge ist, bringt dies durcheinander perfektes Szenario für eine Katastrophe, da wir uns nicht mehr mit Elementen befassen, die in der Vertikalen kritisch auseinander gehen, sodass die Anzahl der Cache-Neuladungen stark reduziert wird.

Ein weiterer Haftungsausschluss – Ich habe mir gerade die Erklärung einfallen lassen und hoffe, dass ich es getroffen habe, aber ich könnte mich irren. Wie auch immer, ich warte auf eine Antwort (oder Bestätigung) von Mystcial. 🙂

Zur Veranschaulichung der Erklärung in Luchian Grigores Antwort sehen Sie hier, wie das Vorhandensein des Matrix-Cache für die beiden Fälle von 64 x 64- und 65 x 65-Matrizen aussieht (Einzelheiten zu Zahlen finden Sie unter dem obigen Link).

Farben in den Animationen unten bedeuten Folgendes:

• – nicht im Cache,
• – im Cache,
• – Cache-Treffer,
• – nur aus RAM lesen,
• – Cache-Miss.

Der 64×64-Fall:

Beachte wie Fast jede Der Zugriff auf eine neue Zeile führt zu einem Cache-Miss. Und wie sieht es nun für den Normalfall, eine 65×65-Matrix aus:

Hier sieht man, dass die meisten Zugriffe nach dem anfänglichen Aufwärmen Cache-Hits sind. So soll der CPU-Cache im Allgemeinen funktionieren.

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

Luchian gibt eine Erklärung warum Dieses Verhalten tritt auf, aber ich dachte, es wäre eine nette Idee, eine mögliche Lösung für dieses Problem zu zeigen und gleichzeitig etwas über Cache-vergessene Algorithmen zu zeigen.

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 großartig, solange Sie diese Details kennen. Nicht besonders portabel.

Können wir das besser? Yes we can: Eine allgemeine Herangehensweise an dieses Problem sind Cache-vergessene Algorithmen das vermeidet, wie der Name schon sagt, die Abhängigkeit von bestimmten Cache-Größen [1]

Die Lösung sähe so aus:

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 bemerkbar ist, das liegt daran, dass wir die iterative Lösung als Blattknoten verwenden, anstatt auf 1 zu rekursiv (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 – that’s inside the margin of error, the fluctuations are in the 100ms area; this “benchmark” isn’t something that I’d be too comfortable with if we wanted completely accurate values])

[1] Quellen für dieses Zeug: Nun, wenn Sie keinen Vortrag von jemandem bekommen können, der mit Leiserson und Co. an diesem Thema gearbeitet hat, nehme ich an, dass ihre Papiere ein guter Ausgangspunkt sind. Diese Algorithmen werden immer noch ziemlich selten beschrieben – CLR hat eine einzige Fußnote über sie. Trotzdem ist es eine großartige Möglichkeit, Leute 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];
            }
        }
    }
}