Matrixmultiplikation: Kleiner Unterschied in der Matrixgröße, großer Unterschied in den Timings

Ich habe einen Matrix-Multiplikationscode, der so aussieht:

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

Hier wird die Größe der Matrix durch dargestellt dimension. Wenn die Größe der Matrizen 2000 beträgt, dauert es 147 Sekunden, um diesen Code auszuführen. Wenn die Größe der Matrizen 2048 beträgt, dauert es 447 Sekunden. Also, während der Unterschied in Nr. von Multiplikationen ist (2048 * 2048 * 2048) / (2000 * 2000 * 2000) = 1,073, der Unterschied in den Timings ist 447/147 = 3. Kann jemand erklären, warum dies passiert? Ich habe erwartet, dass es linear skaliert, was nicht der Fall ist. Ich versuche nicht, den schnellsten Matrix-Multiplikationscode zu erstellen, sondern nur zu verstehen, warum dies geschieht.

Technische Daten: AMD Opteron Dual Core-Knoten (2,2 GHz), 2 G RAM, gcc v 4.5.0

Programm kompiliert als gcc -O3 simple.c

Ich habe dies auch auf Intels icc-Compiler ausgeführt und ähnliche Ergebnisse gesehen.

BEARBEITEN:

Wie in den Kommentaren / Antworten vorgeschlagen, habe ich den Code mit dimension = 2060 ausgeführt und es dauert 145 Sekunden.

Hier ist das komplette Programm:

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>

/* change dimension size as needed */
const int dimension = 2048;
struct timeval tv; 

double timestamp()
{
        double t;
        gettimeofday(&tv, NULL);
        t = tv.tv_sec + (tv.tv_usec/1000000.0);
        return t;
}

int main(int argc, char *argv[])
{
        int i, j, k;
        double *A, *B, *C, start, end;

        A = (double*)malloc(dimension*dimension*sizeof(double));
        B = (double*)malloc(dimension*dimension*sizeof(double));
        C = (double*)malloc(dimension*dimension*sizeof(double));

        srand(292);

        for(i = 0; i < dimension; i++)
                for(j = 0; j < dimension; j++)
                {   
                        A[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                        B[dimension*i+j] = (rand()/(RAND_MAX + 1.0));
                        C[dimension*i+j] = 0.0;
                }   

        start = timestamp();
        for(i = 0; i < dimension; i++)
                for(j = 0; j < dimension; j++)
                        for(k = 0; k < dimension; k++)
                                C[dimension*i+j] += A[dimension*i+k] *
                                        B[dimension*k+j];

        end = timestamp();
        printf("\nsecs:%f\n", end-start);

        free(A);
        free(B);
        free(C);

        return 0;
}

Hier ist meine wilde Vermutung: Cache

Es kann sein, dass Sie 2 Zeilen von 2000 doubles in den Cache einfügen können . Das ist etwas weniger als der 32kb L1-Cache. (beim Verlassen des Raumes andere notwendige Dinge)

Wenn Sie es jedoch auf 2048 erhöhen, wird der gesamte Cache verwendet (und Sie verschütten einige, weil Sie Platz für andere Dinge benötigen).

Angenommen, die Cache-Richtlinie ist LRU. Wenn Sie den Cache nur geringfügig verschütten, wird die gesamte Zeile wiederholt geleert und erneut in den L1-Cache geladen.

Die andere Möglichkeit ist die Cache-Assoziativität aufgrund der Zweierpotenz. Obwohl ich denke, dass der Prozessor 2-Wege-L1-assoziativ ist, denke ich nicht, dass es in diesem Fall wichtig ist. (aber ich werde die Idee trotzdem rauswerfen)

Mögliche Erklärung 2: Konflikt-Cache-Fehler aufgrund von Super-Alignment im L2-Cache.

Ihr BArray wird in der Spalte iteriert. Der Zugang ist also geschritten. Ihre Gesamtdatengröße 2k x 2kbeträgt ca. 32 MB pro Matrix. Das ist viel größer als Ihr L2-Cache.

Wenn die Daten nicht perfekt ausgerichtet sind, haben Sie eine anständige räumliche Lokalität auf B. Obwohl Sie Zeilen hüpfen und nur ein Element pro Cacheline verwenden, verbleibt die Cacheline im L2-Cache, um bei der nächsten Iteration der mittleren Schleife wiederverwendet zu werden.

Wenn die Daten jedoch perfekt ausgerichtet sind (2048), landen diese Hops alle auf demselben "Cache-Weg" und überschreiten Ihre L2-Cache-Assoziativität bei weitem. Daher Bbleiben die Cache-Zeilen, auf die zugegriffen wird, für die nächste Iteration nicht im Cache. Stattdessen müssen sie vollständig vom Widder eingezogen werden.

Sie bekommen definitiv eine sogenannte Cache- Resonanz . Dies ähnelt dem Aliasing , ist jedoch nicht genau dasselbe. Lassen Sie mich erklären.

Caches sind Hardwaredatenstrukturen, die einen Teil der Adresse extrahieren und als Index in einer Tabelle verwenden, ähnlich wie ein Array in der Software. (Tatsächlich nennen wir sie Arrays in der Hardware.) Das Cache-Array enthält Cache-Zeilen mit Daten und Tags - manchmal einen solchen Eintrag pro Index im Array (direkt zugeordnet), manchmal mehrere solcher (N-Wege-Satzassoziativität). Ein zweiter Teil der Adresse wird extrahiert und mit dem im Array gespeicherten Tag verglichen. Index und Tag identifizieren zusammen eine Cache-Zeilen-Speicheradresse eindeutig. Schließlich identifiziert der Rest der Adressbits zusammen mit der Größe des Zugriffs, welche Bytes in der Cache-Zeile adressiert sind.

Normalerweise sind Index und Tag einfache Bitfelder. Eine Speicheradresse sieht also so aus

  ...Tag... | ...Index... | Offset_within_Cache_Line

(Manchmal sind der Index und das Tag Hashes, z. B. einige XORs anderer Bits in den Mittelbereichsbits, die der Index sind. Viel seltener, manchmal sind der Index und seltener das Tag Dinge wie das Aufnehmen der Cache-Zeilenadresse modulo a Primzahl. Diese komplizierteren Indexberechnungen sind Versuche, das hier erläuterte Resonanzproblem zu bekämpfen. Alle leiden unter irgendeiner Form von Resonanz, aber die einfachsten Bitfeldextraktionsschemata leiden unter Resonanz bei allgemeinen Zugriffsmustern, wie Sie festgestellt haben.)

Also, typische Werte ... es gibt viele verschiedene Modelle von "Opteron Dual Core", und ich sehe hier nichts, was angibt, welches Sie haben. Wählen Sie eines nach dem Zufallsprinzip aus, das neueste Handbuch, das ich auf der AMD-Website, im Bios and Kernel Developer's Guide (BKDG) für 15-Stunden-Modelle der AMD-Familie 00h-0Fh , 12. März 2012, sehe .

(Familie 15h = Bulldozer-Familie, der jüngste High-End-Prozessor - die BKDG erwähnt Dual Core, obwohl ich die Produktnummer nicht genau kenne, die Sie beschreiben. Trotzdem gilt für alle Prozessoren dieselbe Resonanzidee. Es ist nur so, dass die Parameter wie Cache-Größe und Assoziativität etwas variieren können.)

Ab S.33:

Der 15-Stunden-Prozessor der AMD-Familie enthält einen 16-KByte-4-Wege-L1-Datencache mit zwei 128-Bit-Ports. Dies ist ein Durchschreibcache, der bis zu zwei 128-Byte-Ladevorgänge pro Zyklus unterstützt. Es ist in 16 Bänke mit einer Breite von jeweils 16 Bytes unterteilt. [...] In einem Zyklus kann nur ein Ladevorgang von einer bestimmten Bank des L1-Cache ausgeführt werden.

Um zusammenzufassen:

• 64-Byte-Cache-Zeile => 6 Offset-Bits innerhalb der Cache-Zeile

• 16KB / 4-Wege => Die Resonanz beträgt 4KB.

  Das heißt, die Adressbits 0-5 sind der Cache-Zeilenversatz.

• 16 KB / 64 KB Cache-Zeilen => 2 ^ 14/2 ^ 6 = 2 ^ 8 = 256 Cache-Zeilen im Cache.
  (Bugfix: Ich habe dies ursprünglich als 128 falsch berechnet, dass ich alle Abhängigkeiten behoben habe.)

• 4-Wege-Assoziativ => 256/4 = 64 Indizes im Cache-Array. Ich (Intel) nenne diese "Sets".

  Das heißt, Sie können den Cache als ein Array von 32 Einträgen oder Sätzen betrachten, wobei jeder Eintrag 4 Cache-Zeilen und ihre Tags enthält. (Es ist komplizierter als das, aber das ist okay).

(Übrigens haben die Begriffe "set" und "way" unterschiedliche Definitionen .)

• Es gibt 6 Indexbits, Bits 6-11 im einfachsten Schema.

  Dies bedeutet, dass alle Cache-Zeilen, die genau dieselben Werte in den Indexbits 6 bis 11 haben, demselben Satz des Caches zugeordnet werden.

Schauen Sie sich jetzt Ihr Programm an.

C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

Schleife k ist die innerste Schleife. Der Basistyp ist doppelt, 8 Bytes. Wenn dimension = 2048, dh 2 KB, sind aufeinanderfolgende Elemente, auf B[dimension*k+j]die die Schleife zugreift, 2048 * 8 = 16 KB voneinander entfernt. Sie werden alle demselben Satz des L1-Cache zugeordnet - sie haben alle denselben Index im Cache. Dies bedeutet, dass anstelle von 256 Cache-Zeilen im Cache nur 4 zur Verfügung stehen - die "4-Wege-Assoziativität" des Caches.

Dh Sie werden wahrscheinlich alle 4 Iterationen um diese Schleife einen Cache-Fehler bekommen. Nicht gut.

(Eigentlich sind die Dinge etwas komplizierter. Aber das Obige ist ein gutes erstes Verständnis. Die Adressen der Einträge von B, die oben erwähnt wurden, sind eine virtuelle Adresse. Es kann also leicht unterschiedliche physikalische Adressen geben. Darüber hinaus verfügt Bulldozer über einen prädiktiven Cache. Wahrscheinlich werden virtuelle Adressbits verwendet, damit nicht auf eine Übersetzung von virtuellen zu physischen Adressen gewartet werden muss. Auf jeden Fall: Ihr Code hat eine "Resonanz" von 16 KB. Der L1-Datencache hat eine Resonanz von 16 KB. Nicht gut .)]

Wenn Sie die Dimension nur geringfügig ändern, z. B. auf 2048 + 1, werden die Adressen von Array B auf alle Sätze des Caches verteilt. Und Sie erhalten deutlich weniger Cache-Fehler.

Es ist eine ziemlich übliche Optimierung, Ihre Arrays aufzufüllen, z. B. 2048 auf 2049 zu ändern, um diese Resonanz zu vermeiden. "Cache-Blockierung ist jedoch eine noch wichtigere Optimierung. Http://suif.stanford.edu/papers/lam-asplos91.pdf

Neben der Cache-Line-Resonanz gibt es hier noch andere Dinge. Beispielsweise hat der L1-Cache 16 Bänke mit einer Breite von jeweils 16 Bytes. Mit der Dimension = 2048 werden aufeinanderfolgende B-Zugriffe in der inneren Schleife immer an dieselbe Bank gesendet. Sie können also nicht parallel geschaltet werden - und wenn der A-Zugriff zufällig auf dieselbe Bank geht, verlieren Sie.

Ich denke nicht, dass dies so groß ist wie die Cache-Resonanz.

Und ja, möglicherweise gibt es ein Aliasing. Beispielsweise vergleicht der STLF (Store To Load Forwarding-Puffer) möglicherweise nur ein kleines Bitfeld und erhält falsche Übereinstimmungen.

(Wenn Sie darüber nachdenken, ist Resonanz im Cache wie Aliasing, das mit der Verwendung von Bitfeldern zusammenhängt. Resonanz wird durch mehrere Cache-Zeilen verursacht, die denselben Satz abbilden und nicht über diese verteilt sind. Alisaing wird durch Matching basierend auf unvollständiger Adresse verursacht Bits.)

Insgesamt meine Empfehlung zur Abstimmung:

1. Versuchen Sie, den Cache ohne weitere Analyse zu blockieren. Ich sage das, weil das Blockieren des Caches einfach ist und es sehr wahrscheinlich ist, dass dies alles ist, was Sie tun müssten.

2. Verwenden Sie danach VTune oder OProf. Oder Cachegrind. Oder ...

3. Besser noch, verwenden Sie eine gut abgestimmte Bibliotheksroutine, um die Matrixmultiplikation durchzuführen.

Es gibt mehrere mögliche Erklärungen. Eine wahrscheinliche Erklärung ist das, was Mysticial vorschlägt: Erschöpfung einer begrenzten Ressource (entweder Cache oder TLB). Eine andere wahrscheinliche Möglichkeit ist ein falscher Aliasing-Stillstand, der auftreten kann, wenn aufeinanderfolgende Speicherzugriffe durch ein Vielfaches einer Zweierpotenz (häufig 4 KB) getrennt sind.

Sie können beginnen, die Arbeit einzugrenzen, indem Sie Zeit / Dimension ^ 3 für einen Wertebereich zeichnen. Wenn Sie einen Cache oder eine erschöpfte TLB-Reichweite gesprengt haben, sehen Sie einen mehr oder weniger flachen Abschnitt, gefolgt von einem starken Anstieg zwischen 2000 und 2048, gefolgt von einem weiteren flachen Abschnitt. Wenn Sie Aliasing-bezogene Stände sehen, sehen Sie ein mehr oder weniger flaches Diagramm mit einer schmalen Spitze nach oben bei 2048.

Dies hat natürlich diagnostische Aussagekraft, ist aber nicht schlüssig. Wenn Sie abschließend wissen möchten , woher die Verlangsamung stammt, sollten Sie sich über Leistungsindikatoren informieren , die diese Art von Frage definitiv beantworten können.

Ich weiß, das ist viel zu alt, aber ich werde einen Bissen nehmen. Es ist (wie gesagt) ein Cache-Problem, das die Verlangsamung bei Zweierpotenzen verursacht. Aber es gibt noch ein anderes Problem: Es ist zu langsam. Wenn Sie sich Ihre Rechenschleife ansehen.

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];

Die innerste Schleife ändert k bei jeder Iteration um 1, was bedeutet, dass Sie nur 1 Doppel vom letzten Element, das Sie von A verwendet haben , zugreifen, aber eine ganze 'Dimension' verdoppelt sich vom letzten Element von B. Dies nutzt keinen Vorteil das Zwischenspeichern der Elemente von B.

Wenn Sie dies ändern in:

for(i = 0; i < dimension; i++)
    for(j = 0; j < dimension; j++)
        for(k = 0; k < dimension; k++)
            C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];

Sie erhalten genau die gleichen Ergebnisse (Modulo-Doppeladditions-Assoziativitätsfehler), aber es ist viel cachefreundlicher ( lokal ). Ich habe es versucht und es gibt wesentliche Verbesserungen. Dies kann zusammengefasst werden als

Multiplizieren Sie Matrizen nicht per Definition, sondern mit Zeilen

Beispiel für die Beschleunigung (Ich habe Ihren Code geändert, um die Dimension als Argument zu verwenden.)

$ diff a.c b.c
42c42
<               C[dimension*i+j] += A[dimension*i+k] * B[dimension*k+j];
---
>               C[dimension*i+k] += A[dimension*i+j] * B[dimension*j+k];
$ make a
cc     a.c   -o a
$ make b
cc     b.c   -o b
$ ./a 1024

secs:88.732918
$ ./b 1024

secs:12.116630

Als Bonus (und was dies mit dieser Frage zusammenhängt) ist, dass diese Schleife nicht unter dem vorherigen Problem leidet.

Wenn Sie das alles schon gewusst haben, dann entschuldige ich mich!

In einigen Antworten wurden Probleme mit dem L2-Cache erwähnt.

Sie können dies tatsächlich mit einer Cache- Simulation überprüfen . Valgrinds Cachegrind- Tool kann das.

valgrind --tool=cachegrind --cache-sim=yes your_executable

Stellen Sie die Befehlszeilenparameter so ein, dass sie mit den L2-Parametern Ihrer CPU übereinstimmen.

Wenn Sie es mit verschiedenen Matrixgrößen testen, werden Sie wahrscheinlich einen plötzlichen Anstieg der L2-Fehlerquote feststellen.