Frage Wie wird BLAS so extrem leistungsfähig?


Aus Neugier entschied ich mich, meine Matrixmultiplikationsfunktion mit der BLAS-Implementierung zu vergleichen ... Ich war am wenigsten überrascht über das Ergebnis:

Benutzerdefinierte Implementierung, 10 Versuche von   1000x1000 Matrixmultiplikation:

Took: 15.76542 seconds.

BLAS Implementierung, 10 Versuche von   1000x1000 Matrixmultiplikation:

Took: 1.32432 seconds.

Dies verwendet Gleitkommazahlen mit einfacher Genauigkeit.

Meine Implementierung:

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

Ich habe zwei Fragen:

  1. Angenommen, eine Matrix-Matrix-Multiplikation sagt: nxm * mxn erfordert n * n * m Multiplikationen, also im Fall über 1000 ^ 3 oder 1e9-Operationen. Wie ist es möglich, auf meinem 2.6Ghz Prozessor für BLAS 10 * 1e9 Operationen in 1.32 Sekunden zu machen? Selbst wenn Multiplikationen eine einzelne Operation waren und nichts anderes getan wurde, sollte es ~ 4 Sekunden dauern.
  2. Warum ist meine Implementierung so viel langsamer?

75
2017-08-19 23:30


Ursprung


Antworten:


Aus vielen Gründen.

Erstens, Fortran-Compiler sind hoch optimiert, und die Sprache erlaubt ihnen, so zu sein. C und C ++ sind in Bezug auf die Array-Handhabung sehr lose (z. B. im Falle von Zeigern, die sich auf denselben Speicherbereich beziehen). Dies bedeutet, dass der Compiler nicht im Voraus wissen kann, was zu tun ist, und gezwungen ist, generischen Code zu erstellen. In Fortran sind Ihre Fälle rationeller und der Compiler hat eine bessere Kontrolle darüber, was passiert, wodurch er mehr optimieren kann (z. B. mithilfe von Registern).

Eine andere Sache ist, dass Fortran Sachen spaltenweise speichert, während C Daten zeilenweise speichert. Ich habe Ihren Code nicht überprüft, aber seien Sie vorsichtig, wie Sie das Produkt ausführen. In C müssen Sie zeilenweise scannen: Auf diese Weise scannen Sie Ihr Array entlang des zusammenhängenden Speichers und reduzieren so die Cache-Fehler. Cache-Miss ist die erste Quelle von Ineffizienz.

Drittens hängt es von der blas-Implementierung ab, die Sie verwenden. Einige Implementierungen werden möglicherweise in Assembler geschrieben und für den spezifischen Prozessor optimiert, den Sie verwenden. Die Netlib-Version ist in Fortran 77 geschrieben.

Außerdem machen Sie viele Operationen, die meisten von ihnen wiederholt und redundant. Alle diese Multiplikationen, um den Index zu erhalten, sind für die Leistung nachteilig. Ich weiß nicht genau, wie das bei BLAS gemacht wird, aber es gibt viele Tricks, um teure Operationen zu verhindern.

Zum Beispiel könnten Sie Ihren Code auf diese Weise überarbeiten

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
    a1 = cc2*ADim2;
    a3 = cc2*BDim1
    for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
          a2=cc1*ADim1;
          ValT b = B[a3+cc1];
          for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
                    C[a1+cr1] += A[a2+cr1]*b;
           }
     }
  }
} 

Probieren Sie es aus, ich bin mir sicher, dass Sie etwas sparen werden.

Auf Ihrer # 1 Frage liegt der Grund darin, dass Matrixmultiplikation als O (n ^ 3) skaliert, wenn Sie einen trivialen Algorithmus verwenden. Es gibt Algorithmen, die viel besser skalieren.


-24
2017-08-19 23:36



Ein guter Ausgangspunkt ist das großartige Buch Die Wissenschaft der Programmierung von Matrixberechnungen von Robert A. van de Geijn und Enrique S. Quintana-Ortí. Sie bieten eine kostenlose Download-Version.

BLAS gliedert sich in drei Ebenen:

  • Ebene 1 definiert eine Reihe linearer Algebrafunktionen, die nur mit Vektoren arbeiten. Diese Funktionen profitieren von der Vektorisierung (z.B. von der Verwendung von SSE).

  • Level 2-Funktionen sind Matrix-Vektor-Operationen, z. etwas Matrix-Vektor-Produkt. Diese Funktionen könnten in Bezug auf Level1-Funktionen implementiert werden. Sie können jedoch die Leistung dieser Funktionen steigern, wenn Sie eine dedizierte Implementierung bereitstellen können, die eine Multiprozessorarchitektur mit gemeinsam genutztem Speicher verwendet.

  • Level-3-Funktionen sind Operationen wie das Matrix-Matrix-Produkt. Wiederum könnten Sie sie in Bezug auf Level2-Funktionen implementieren. Aber Level3-Funktionen führen O (N ^ 3) -Operationen an O (N ^ 2) -Daten durch. Wenn Ihre Plattform über eine Cache-Hierarchie verfügt, können Sie die Leistung steigern, wenn Sie eine dedizierte Implementierung bereitstellen Cache-optimiert / Cache-freundlich. Dies ist in dem Buch gut beschrieben. Der Hauptschub von Level3-Funktionen kommt von der Cache-Optimierung. Dieser Boost übersteigt deutlich den zweiten Schub von Parallelität und anderen Hardware-Optimierungen.

Übrigens sind die meisten (oder sogar alle) der Hochleistungs-BLAS-Implementierungen NICHT in Fortran implementiert. ATLAS ist in C. implementiert. GotoBLAS / OpenBLAS ist in C implementiert und seine leistungskritischen Teile in Assembler. In Fortran ist nur die Referenzimplementierung von BLAS implementiert. Alle diese BLAS-Implementierungen bieten jedoch eine Fortran-Schnittstelle, so dass sie mit LAPACK verknüpft werden kann (LAPACK erhält seine gesamte Leistung von BLAS).

Optimierte Compiler spielen dabei eine untergeordnete Rolle (und für GotoBLAS / OpenBLAS spielt der Compiler keine Rolle).

IMHO no BLAS-Implementierung verwendet Algorithmen wie den Coppersmith-Winograd-Algorithmus oder den Strassen-Algorithmus. Ich bin mir über den Grund nicht ganz sicher, aber das ist meine Vermutung:

  • Vielleicht ist es nicht möglich, eine Cache-optimierte Implementierung dieser Algorithmen bereitzustellen (d. H. Sie würden mehr verlieren als Sie gewinnen würden)
  • Diese Algorithmen sind numerisch nicht stabil. Da BLAS der Rechenkern von LAPACK ist, ist dies ein No-Go.

Bearbeiten / Aktualisieren:

Das neue und bahnbrechende Papier zu diesem Thema sind die BLIS Papiere. Sie sind außergewöhnlich gut geschrieben. Für meinen Vortrag "Software Basics for High Performance Computing" habe ich das Matrix-Matrix-Produkt nach der Arbeit implementiert. Tatsächlich habe ich mehrere Varianten des Matrix-Matrix-Produkts implementiert. Die einfachste Variante ist vollständig in C geschrieben und hat weniger als 450 Zeilen Code. Alle anderen Varianten optimieren lediglich die Schleifen

    for (l=0; l<MR*NR; ++l) {
        AB[l] = 0;
    }
    for (l=0; l<kc; ++l) {
        for (j=0; j<NR; ++j) {
            for (i=0; i<MR; ++i) {
                AB[i+j*MR] += A[i]*B[j];
            }
        }
        A += MR;
        B += NR;
    }

Die Gesamtleistung des Matrix-Matrix-Produkts nur hängt von diesen Schleifen ab. Ungefähr 99,9% der Zeit werden hier verbracht. In den anderen Varianten nutzte ich Intrinsics und Assembler-Code, um die Leistung zu verbessern. Sie können das Tutorial hier durch alle Varianten sehen:

ulmBLAS: Tutorial zu GEMM (Matrix-Matrix-Produkt)

Zusammen mit den BLIS-Papieren wird es relativ einfach zu verstehen, wie Bibliotheken wie Intel MKL eine solche Leistung erzielen können. Und warum ist es egal, ob Sie Zeilen- oder Spaltenspeicher verwenden!

Die letzten Benchmarks sind hier (wir nannten unser Projekt ulmBLAS):

Benchmarks für ulmBLAS, BLIS, MKL, openBLAS und Eigen

Ein weiteres Edit / Update:

Ich habe auch ein Tutorial geschrieben, wie BLAS für numerische lineare Algebra-Probleme wie das Lösen eines linearen Gleichungssystems verwendet wird:

Hochleistungs-LU-Faktorisierung

(Diese LU-Faktorisierung wird beispielsweise von Matlab zur Lösung eines linearen Gleichungssystems verwendet.)

Ich hoffe, Zeit zu finden das Tutorial zu erweitern, um zu beschreiben und zu demonstrieren, wie eine hochskalierbare parallele Implementierung der LU - Faktorisierung wie in realisiert werden kann PLASMA.

Ok, hier hast du es: Codieren eines Cache Optimierte parallele LU-Faktorisierung

P.S .: Ich habe auch einige Experimente gemacht, um die Leistung von uBLAS zu verbessern. Es ist eigentlich ziemlich einfach, die Leistung von uBLAS zu steigern (yeah, play on words :))

Experimente auf uBLAS.

Hier ein ähnliches Projekt mit BLASEN:

Experimente auf BLAZE.


98
2017-07-10 20:23



BLAS ist also zunächst nur eine Schnittstelle von etwa 50 Funktionen. Es gibt viele konkurrierende Implementierungen der Schnittstelle.

Zunächst werde ich Dinge erwähnen, die weitgehend unabhängig sind:

  • Fortran vs C, macht keinen Unterschied
  • Fortgeschrittene Matrixalgorithmen wie Strassen, Implementierungen benutzen sie nicht, da sie in der Praxis nicht helfen

Die meisten Implementierungen unterbrechen jede Operation in mehrdimensionaler Matrix- oder Vektoroperation auf mehr oder weniger offensichtliche Weise. Zum Beispiel kann eine große 1000 × 1000 Matrixmultiplikation in eine Sequenz von 50 × 50 Matrixmultiplikationen aufgeteilt werden.

Diese kleinformatigen Operationen mit fester Größe (die als Kernel bezeichnet werden) sind im CPU-spezifischen Assemblercode unter Verwendung mehrerer CPU-Funktionen ihres Ziels fest codiert:

  • SIMD-ähnliche Anweisungen
  • Parallelität der Anweisungsstufe
  • Cache-Bewusstsein

Darüber hinaus können diese Kernel parallel zueinander unter Verwendung mehrerer Threads (CPU-Kerne) in dem typischen Map-Reduced-Entwurfsmuster ausgeführt werden.

Sehen Sie sich ATLAS an, die am häufigsten verwendete Open-Source-BLAS-Implementierung. Es hat viele verschiedene konkurrierende Kernel und während des ATLAS-Library-Build-Prozesses läuft ein Wettbewerb zwischen ihnen (einige sind sogar parametrisiert, so dass der gleiche Kernel verschiedene Einstellungen haben kann). Es versucht verschiedene Konfigurationen und wählt dann das Beste für das jeweilige Zielsystem aus.

(Tipp: Wenn Sie ATLAS verwenden, sollten Sie daher die Bibliothek für Ihre spezielle Maschine manuell erstellen und abstimmen, bevor Sie eine vordefinierte Version verwenden.)


19
2018-06-26 14:10



Erstens gibt es effizientere Algorithmen für die Matrixmultiplikation als die, die Sie verwenden.

Zweitens kann Ihre CPU viel mehr als eine Anweisung gleichzeitig ausführen.

Ihre CPU führt 3-4 Befehle pro Zyklus aus, und wenn die SIMD-Einheiten verwendet werden, verarbeitet jede Anweisung 4 Floats oder 2 Doubles. (Natürlich ist diese Zahl auch nicht genau, da die CPU typischerweise nur einen SIMD-Befehl pro Zyklus verarbeiten kann)

Drittens ist Ihr Code weit davon entfernt, optimal zu sein:

  • Sie verwenden rohe Zeiger, was bedeutet, dass der Compiler davon ausgehen muss, dass sie einen Aliasnamen haben. Es gibt compilerspezifische Schlüsselwörter oder Flags, die Sie angeben können, um dem Compiler mitzuteilen, dass sie keinen Aliasnamen haben. Alternativ sollten Sie andere Typen als rohe Zeiger verwenden, die sich um das Problem kümmern.
  • Sie dreschen den Cache durch eine naive Durchquerung jeder Zeile / Spalte der Eingabematrizen. Sie können die Blockierung verwenden, um auf einem kleineren Block der Matrix, der in den CPU-Cache passt, so viel Arbeit wie möglich auszuführen, bevor Sie mit dem nächsten Block fortfahren.
  • Für rein numerische Aufgaben ist Fortran ziemlich unschlagbar, und C ++ braucht viel Überredungskunst, um eine ähnliche Geschwindigkeit zu erreichen. Es kann gemacht werden, und es gibt ein paar Bibliotheken, die es demonstrieren (normalerweise mit Expression-Templates), aber es ist nicht trivial und tut es auch nicht gerade geschehen.

12
2017-08-20 12:12



Ich weiß nicht genau über BLAS-Implementierung, aber es gibt effizientere Algorithmen für Matrix-Multiplikation, die besser als O (n3) -Komplexität ist. Ein bekannter ist es Straßenalgorithmus 


9
2017-08-20 00:15



Die meisten Argumente für die zweite Frage - Assembler, Aufspalten in Blöcke usw. (aber nicht die weniger als N ^ 3 Algorithmen, sie sind wirklich überentwickelt) - spielen eine Rolle. Aber die niedrige Geschwindigkeit Ihres Algorithmus wird im Wesentlichen durch die Matrixgröße und die unglückliche Anordnung der drei verschachtelten Schleifen verursacht. Ihre Matrizen sind so groß, dass sie nicht sofort in den Cache-Speicher passen. Sie können die Schleifen so anordnen, dass so viel wie möglich in einer Zeile im Cache ausgeführt wird, wodurch Cache-Aktualisierungen erheblich reduziert werden (die Aufteilung in kleine Blöcke hat einen analogen Effekt, am besten, wenn Schleifen über die Blöcke ähnlich angeordnet sind). Es folgt eine Modellimplementierung für quadratische Matrizen. Auf meinem Computer war der Zeitaufwand etwa 1:10 gegenüber der Standardimplementierung (wie bei Ihnen). Mit anderen Worten: programmieren Sie niemals eine Matrixmultiplikation entlang des Schemas "row times column", das wir in der Schule gelernt haben. Nachdem die Schleifen neu angeordnet wurden, werden weitere Verbesserungen durch Abwickeln von Schleifen, Assemblercode usw. erreicht.

    void vector(int m, double ** a, double ** b, double ** c) {
      int i, j, k;
      for (i=0; i<m; i++) {
        double * ci = c[i];
        for (k=0; k<m; k++) ci[k] = 0.;
        for (j=0; j<m; j++) {
          double aij = a[i][j];
          double * bj = b[j];
          for (k=0; k<m; k++)  ci[k] += aij*bj[k];
        }
      }
    }

Eine weitere Bemerkung: Diese Implementierung ist auf meinem Computer sogar noch besser, als alle durch die BLAS-Routine cblas_dgemm zu ersetzen (versuchen Sie es auf Ihrem Computer!). Aber viel schneller (1: 4) ruft dgemm_ der Fortran-Bibliothek direkt auf. Ich denke, dass diese Routine in der Tat nicht Fortran, sondern Assembler-Code ist (ich weiß nicht, was in der Bibliothek ist, ich habe keine Quellen). Völlig unklar für mich ist, warum cblas_dgemm nicht so schnell ist, da es meines Wissens nur ein Wrapper für dgemm_ ist.


4
2017-11-30 20:11



Das ist eine realistische Beschleunigung. Ein Beispiel dafür, was mit SIMD Assembler über C ++ Code gemacht werden kann, finden Sie in einem Beispiel iPhone Matrixfunktionen - Diese waren über 8x schneller als die C-Version und sind noch nicht einmal "optimierte" Assembly - es gibt noch kein Pipeline-Lining und es gibt unnötige Stack-Operationen.

Auch dein Code ist nicht "richtig beschränken"- Wie weiß der Compiler, dass er, wenn er C modifiziert, A und B nicht modifiziert?


3
2017-08-20 00:10



In Bezug auf den ursprünglichen Code in MM multiplizieren, Speicherreferenz für die meisten Operationen ist der Hauptgrund für schlechte Leistung. Der Speicher läuft 100-1000 mal langsamer als der Cache.

Die meiste Geschwindigkeit kommt von der Verwendung von Schleifenoptimierungstechniken für diese Dreifachschleifenfunktion in MM multiplizieren. Zwei Hauptschleifenoptimierungstechniken werden verwendet; entrollen und blockieren. In Bezug auf das Abrollen rollen wir die äußersten zwei Schleifen ab und blockieren sie für die Datenwiederverwendung im Cache. Durch das Abwickeln der äußeren Schleife wird der Datenzugriff zeitlich optimiert, indem die Anzahl der Speicherreferenzen für die gleichen Daten zu verschiedenen Zeiten während des gesamten Vorgangs reduziert wird. Das Blockieren des Schleifenindex um eine bestimmte Nummer hilft dabei, die Daten im Cache zu behalten. Sie können wählen, ob Sie für L2-Cache oder L3-Cache optimieren möchten.

https://en.wikipedia.org/wiki/Loop_nest_optimization


1
2018-05-02 12:07