Wie kommt BLAS zu solch extremer Leistung?

108

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

Benutzerdefinierte Implementierung, 10 Versuche zur 1000x1000-Matrixmultiplikation:

Took: 15.76542 seconds.

BLAS-Implementierung, 10 Versuche zur 1000x1000-Matrixmultiplikation:

Took: 1.32432 seconds.

Hierbei werden Gleitkommazahlen mit einfacher Genauigkeit verwendet.

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. Vorausgesetzt, eine Matrix-Matrix-Multiplikation sagt: nxm * mxn erfordert n * n * m Multiplikationen, also im Fall über 1000 ^ 3 oder 1e9 Operationen. Wie ist es auf meinem 2,6-GHz-Prozessor möglich, dass BLAS 10 * 1e9-Operationen in 1,32 Sekunden ausführt? Selbst wenn Multiplikationen eine einzelne Operation wären und nichts anderes getan würde, sollte es ~ 4 Sekunden dauern.
  2. Warum ist meine Implementierung so viel langsamer?
DeusAduro
quelle
17
BLAS wurde von Fachleuten auf der einen Seite und auf der anderen Seite optimiert. Ich
gehe
3
Wie können Sie dennoch 1E10-Operationen auf einem Prozessor mit 2,63E9 Zyklen / Sekunde in 1,3 Sekunden ausführen?
DeusAduro
9
Mehrere Ausführungseinheiten, Pipe-Lining und Single Instruction Multiple Data (SIMD) (dh mehrere Operationen gleichzeitig mit mehr als einem Operandenpaar). Einige Compiler können die SIMD-Einheiten auf gängigen Chips ausrichten, aber Sie müssen sie fast immer explizit einschalten, und es hilft zu wissen, wie alles funktioniert ( en.wikipedia.org/wiki/SIMD ). Die Versicherung gegen Cache-Fehler ist mit ziemlicher Sicherheit der schwierige Teil.
dmckee --- Ex-Moderator Kätzchen
13
Annahme ist falsch. Es sind bessere Algorithmen bekannt, siehe Wikipedia.
MSalters
2
@DeusAduro: In meiner Antwort für Wie schreibe ich ein Matrixmatrixprodukt, das mit Eigen konkurrieren kann? Ich habe ein kleines Beispiel zur Implementierung eines cache-effizienten Matrix-Matrix-Produkts veröffentlicht.
Michael Lehn

Antworten:

141

Ein guter Ausgangspunkt ist das großartige Buch The Science of Programming Matrix Computations von Robert A. van de Geijn und Enrique S. Quintana-Ortí. Sie bieten eine kostenlose Download-Version.

BLAS ist in drei Ebenen unterteilt:

  • Stufe 1 definiert eine Reihe von linearen Algebra-Funktionen, die nur mit Vektoren arbeiten. Diese Funktionen profitieren von der Vektorisierung (z. B. von der Verwendung von SSE).

  • Funktionen der Ebene 2 sind Matrixvektoroperationen, z. B. ein Matrixvektorprodukt. Diese Funktionen könnten in Form von Level1-Funktionen implementiert werden. Sie können jedoch die Leistung dieser Funktionen steigern, wenn Sie eine dedizierte Implementierung bereitstellen, die eine Multiprozessorarchitektur mit gemeinsam genutztem Speicher verwendet.

  • Level 3-Funktionen sind Operationen wie das Matrix-Matrix-Produkt. Auch hier können Sie sie in Bezug auf Level2-Funktionen implementieren. Level3-Funktionen führen jedoch O (N ^ 3) -Operationen für O (N ^ 2) -Daten aus. Wenn Ihre Plattform über eine Cache-Hierarchie verfügt, können Sie die Leistung steigern, wenn Sie eine dedizierte Implementierung bereitstellen, die cache-optimiert / cache-freundlich ist . Dies ist im Buch gut beschrieben. Der Hauptschub der Level3-Funktionen liegt in der Cache-Optimierung. Dieser Boost übertrifft den zweiten Boost aufgrund von Parallelität und anderen Hardwareoptimierungen erheblich.

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

Optimierte Compiler spielen in dieser Hinsicht eine untergeordnete Rolle (und für GotoBLAS / OpenBLAS spielt der Compiler überhaupt keine Rolle).

IMHO keine 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:

  • Möglicherweise ist es nicht möglich, eine Cache-optimierte Implementierung dieser Algorithmen bereitzustellen (dh Sie würden mehr verlieren als gewinnen).
  • 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 meine Vorlesung "Software-Grundlagen für High Performance Computing" habe ich das Matrix-Matrix-Produkt nach deren Arbeit implementiert. Eigentlich habe ich mehrere Varianten des Matrix-Matrix-Produkts implementiert. Die einfachste Variante ist vollständig in einfachem C geschrieben und enthält weniger als 450 Codezeilen. 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 hängt nur von diesen Schleifen ab. Etwa 99,9% der Zeit wird hier verbracht. In den anderen Varianten habe ich Intrinsics und Assembler-Code verwendet, um die Leistung zu verbessern. Sie können das Tutorial sehen, das alle Varianten hier durchläuft:

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

Zusammen mit den BLIS-Papieren wird es ziemlich einfach zu verstehen, wie Bibliotheken wie Intel MKL eine solche Leistung erzielen können. Und warum spielt es keine Rolle, ob Sie Zeilen- oder Spalten-Hauptspeicher verwenden!

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

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

Noch eine Bearbeitung / Aktualisierung:

Ich habe auch ein Tutorial darüber 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 zum Lösen eines linearen Gleichungssystems verwendet.)

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

Ok, los geht's: Codieren einer Cache-optimierten parallelen LU-Faktorisierung

PS: Ich habe auch einige Experimente zur Verbesserung der Leistung von uBLAS durchgeführt. Es ist eigentlich ziemlich einfach, die Leistung von uBLAS zu steigern (ja, mit Worten spielen :)):

Experimente zu uBLAS .

Hier ein ähnliches Projekt mit BLAZE :

Experimente mit BLAZE .

Michael Lehn
quelle
3
Neuer Link zu "Benchmarks für ulmBLAS, BLIS, MKL, openBLAS und Eigen": apfel.mathematik.uni-ulm.de/~lehn/ulmBLAS/#toc3
Ahmed Fasih
Es stellt sich heraus, dass IBMs ESSL eine Variation des Strassen-Algorithmus verwendet - ibm.com/support/knowledgecenter/en/SSFHY8/essl_welcome.html
ben-albrecht
2
die meisten Links sind tot
Aurélien Pierre
Ein PDF von TSoPMC finden Sie auf der Seite des Autors unter cs.utexas.edu/users/rvdg/tmp/TSoPMC.pdf
Alex Shpilkin
Obwohl der Coppersmith-Winograd-Algorithmus auf dem Papier eine schöne zeitliche Komplexität aufweist, verbirgt die Big O-Notation eine sehr große Konstante, sodass sie erst für lächerlich große Matrizen realisierbar wird.
DiehardTheTryhard
26

Zuallererst ist BLAS also nur eine Schnittstelle von ungefähr 50 Funktionen. Es gibt viele konkurrierende Implementierungen der Schnittstelle.

Zunächst möchte ich Dinge erwähnen, die weitgehend unabhängig sind:

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

Die meisten Implementierungen unterteilen jede Operation auf mehr oder weniger offensichtliche Weise in Matrix- oder Vektoroperationen mit kleinen Dimensionen. Zum Beispiel kann eine große 1000x1000-Matrixmultiplikation in eine Folge von 50x50-Matrixmultiplikationen aufgeteilt werden.

Diese kleinen Operationen mit fester Größe (Kernel genannt) werden im CPU-spezifischen Assemblycode unter Verwendung mehrerer CPU-Funktionen ihres Ziels fest codiert:

  • Anweisungen im SIMD-Stil
  • Parallelität auf Befehlsebene
  • Cache-Bewusstsein

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

Schauen Sie sich ATLAS an, die am häufigsten verwendete Open-Source-BLAS-Implementierung. Es hat viele verschiedene konkurrierende Kernel und während des Erstellungsprozesses der ATLAS-Bibliothek wird ein Wettbewerb zwischen ihnen ausgeführt (einige sind sogar parametrisiert, sodass derselbe Kernel unterschiedliche Einstellungen haben kann). Es werden verschiedene Konfigurationen ausprobiert und dann die beste für das jeweilige Zielsystem ausgewählt.

(Tipp: Wenn Sie ATLAS verwenden, ist es daher besser, die Bibliothek von Hand für Ihren bestimmten Computer zu erstellen und zu optimieren, als einen vorgefertigten.)

Andrew Tomazos
quelle
ATLAS ist nicht mehr die am häufigsten verwendete Open-Source-BLAS-Implementierung. Es wurde von OpenBLAS (eine Gabelung des GotoBLAS) und BLIS (ein Refactoring des GotoBLAS) übertroffen.
Robert van de Geijn
1
@ ulaff.net: Das vielleicht. Dies wurde vor 6 Jahren geschrieben. Ich denke, die derzeit schnellste BLAS-Implementierung (natürlich auf Intel) ist Intel MKL, aber es ist kein Open Source.
Andrew Tomazos
14

Erstens gibt es effizientere Algorithmen für die Matrixmultiplikation als den von Ihnen verwendeten.

Zweitens kann Ihre CPU viel mehr als einen Befehl gleichzeitig ausführen.

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

Drittens ist Ihr Code alles andere als optimal:

  • Sie verwenden unformatierte Zeiger, was bedeutet, dass der Compiler davon ausgehen muss, dass es sich möglicherweise um Alias ​​handelt. Es gibt compilerspezifische Schlüsselwörter oder Flags, die Sie angeben können, um dem Compiler mitzuteilen, dass sie keinen Alias ​​haben. Alternativ sollten Sie andere Typen als Rohzeiger verwenden, um das Problem zu beheben.
  • Sie verdrängen den Cache, indem Sie eine naive Durchquerung jeder Zeile / Spalte der Eingabematrizen durchführen. Sie können das Blockieren verwenden, um so viel Arbeit wie möglich an einem kleineren Block der Matrix auszuführen, der in den CPU-Cache passt, bevor Sie mit dem nächsten Block fortfahren.
  • Für rein numerische Aufgaben ist Fortran so gut wie unschlagbar, und C ++ erfordert viel Überredung, um eine ähnliche Geschwindigkeit zu erreichen. Dies ist möglich, und es gibt einige Bibliotheken, die dies demonstrieren (normalerweise unter Verwendung von Ausdrucksvorlagen), aber es ist nicht trivial und es passiert nicht einfach so .
jalf
quelle
Vielen Dank, ich habe den korrekten Code gemäß Justicles Vorschlag hinzugefügt und keine große Verbesserung festgestellt. Ich mag die blockweise Idee. Aus Neugier, ohne die Cache-Größe der CPU zu kennen, wie würde man den richtigen optimalen Code finden?
DeusAduro
2
Das tust du nicht. Um optimalen Code zu erhalten, müssen Sie die Cache-Größe der CPU kennen. Der Nachteil dabei ist natürlich, dass Sie Ihren Code effektiv hartcodieren, um die beste Leistung auf einer CPU-Familie zu erzielen .
Jalf
2
Zumindest die innere Schleife vermeidet hier Schrittlasten. Es sieht so aus, als ob dies für eine Matrix geschrieben wurde, die bereits transponiert wurde. Deshalb ist es "nur" eine Größenordnung langsamer als BLAS! Aber ja, es schlägt immer noch zu, weil es keine Cache-Blockierung gibt. Sind Sie sicher, dass Fortran viel helfen würde? Ich denke, alles, was Sie hier gewinnen würden, ist, dass restrict(kein Aliasing) die Standardeinstellung ist, anders als in C / C ++. (Und leider hat ISO C ++ kein restrictSchlüsselwort, so dass Sie es __restrict__auf Compilern verwenden müssen, die es als Erweiterung bereitstellen).
Peter Cordes
11

Ich weiß nichts über die BLAS-Implementierung, aber es gibt effizientere Alogorithmen für die Matrixmultiplikation, die besser als die Komplexität von O (n3) sind. Ein bekannter ist der Strassen-Algorithmus

Softveda
quelle
8
Der Strassen-Algorithmus wird in der Numerik aus zwei Gründen nicht verwendet: 1) Er ist nicht stabil. 2) Sie sparen einige Berechnungen, aber das ist mit dem Preis verbunden, den Sie für Cache-Hierarchien verwenden können. In der Praxis verlieren Sie sogar die Leistung.
Michael Lehn
4
Für die praktische Implementierung des Strassen-Algorithmus, der eng auf dem Quellcode der BLAS-Bibliothek basiert, gibt es eine kürzlich erschienene Veröffentlichung: " Strassen Algorithm Reloaded " in SC16, die selbst bei der Problemgröße 1000x1000 eine höhere Leistung als BLAS erzielt.
Jianyu Huang
4

Die meisten Argumente für die zweite Frage - Assembler, Aufteilung in Blöcke usw. (aber nicht weniger als N ^ 3-Algorithmen, sie sind wirklich überentwickelt) - spielen eine Rolle. Die geringe Geschwindigkeit Ihres Algorithmus wird jedoch 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 neu anordnen, dass so viel wie möglich in einer Zeile im Cache ausgeführt wird, wodurch die Cache-Aktualisierungen drastisch reduziert werden (BTW-Aufteilung in kleine Blöcke hat einen analogen Effekt, am besten, wenn die Schleifen über den Blöcken ähnlich angeordnet sind). Es folgt eine Modellimplementierung für quadratische Matrizen. Auf meinem Computer betrug der Zeitaufwand etwa 1:10 im Vergleich zur Standardimplementierung (wie bei Ihnen). Mit anderen Worten: Programmieren Sie niemals eine Matrixmultiplikation entlang der "

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

Noch eine Bemerkung: Diese Implementierung ist auf meinem Computer noch besser als das Ersetzen aller durch die BLAS-Routine cblas_dgemm (versuchen Sie es auf Ihrem Computer!). Viel schneller (1: 4) wird jedoch dgemm_ der Fortran-Bibliothek direkt aufgerufen. Ich denke, diese Routine ist in der Tat nicht Fortran, sondern Assembler-Code (ich weiß nicht, was sich in der Bibliothek befindet, ich habe keine Quellen). Völlig unklar ist für mich, warum cblas_dgemm nicht so schnell ist, da es meines Wissens nur ein Wrapper für dgemm_ ist.

Wolfgang Jansen
quelle
3

Dies ist eine realistische Beschleunigung. Ein Beispiel dafür, was mit SIMD Assembler über C ++ - Code gemacht werden kann, finden Sie in einigen Beispielfunktionen der iPhone-Matrix - diese waren über 8x schneller als die C-Version und nicht einmal "optimiert" - es gibt noch keine Rohrverkleidung ist unnötige Stapeloperationen.

Außerdem ist Ihr Code nicht " korrekt einschränken " - woher weiß der Compiler, dass er beim Ändern von C A und B nicht ändert?

Justicle
quelle
Sicher, wenn Sie die Funktion wie mmult (A ..., A ..., A) aufgerufen haben; Sie würden sicherlich nicht das erwartete Ergebnis erhalten. Auch hier habe ich nicht versucht, BLAS zu schlagen / neu zu implementieren, sondern nur gesehen, wie schnell es wirklich ist. Daher wurde die Fehlerprüfung nicht berücksichtigt, sondern nur die grundlegende Funktionalität.
DeusAduro
3
Es tut mir leid, um klar zu sein, ich sage, wenn Sie Ihre Zeiger "einschränken", erhalten Sie viel schnelleren Code. Dies liegt daran, dass der Compiler jedes Mal, wenn Sie C ändern, A und B nicht neu laden muss, was die innere Schleife erheblich beschleunigt. Wenn Sie mir nicht glauben, überprüfen Sie die Demontage.
Justicle
@DeusAduro: Dies ist keine Fehlerprüfung. Möglicherweise kann der Compiler die Zugriffe auf das B [] -Array in der inneren Schleife nicht optimieren, da er möglicherweise nicht herausfinden kann, dass die Zeiger A und C niemals den Alias ​​B haben Array. Wenn es ein Aliasing gäbe, könnte sich der Wert im B-Array ändern, während die innere Schleife ausgeführt wird. Wenn Sie den Zugriff auf den B [] -Wert aus der inneren Schleife herausheben und in eine lokale Variable einfügen, kann der Compiler möglicherweise kontinuierliche Zugriffe auf B [] vermeiden.
Michael Burr
1
Hmmm, also habe ich zuerst versucht, das Schlüsselwort '__restrict' in VS 2008 zu verwenden, das auf A, B und C angewendet wurde. Dies zeigte keine Änderung im Ergebnis. Durch Verschieben des Zugriffs auf B von der innersten Schleife zur äußeren Schleife wurde die Zeit jedoch um ~ 10% verbessert.
DeusAduro
1
Entschuldigung, ich bin mir bei VC nicht sicher, aber mit GCC müssen Sie aktivieren -fstrict-aliasing. Es gibt auch eine bessere Erklärung für "einschränken" hier: cellperformance.beyond3d.com/articles/2006/05/…
Justicle
2

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

Der größte Teil der Beschleunigung ergibt sich aus der Verwendung von Schleifenoptimierungstechniken für diese Dreifachschleifenfunktion bei der MM-Multiplikation. Es werden zwei Hauptschleifenoptimierungstechniken verwendet; Abrollen und Blockieren. In Bezug auf das Abrollen entrollen wir die beiden äußersten Schleifen und blockieren sie für die Wiederverwendung von Daten im Cache. Das Abrollen der äußeren Schleife hilft, den Datenzugriff zeitlich zu optimieren, indem die Anzahl der Speicherreferenzen auf dieselben Daten zu unterschiedlichen Zeiten während des gesamten Vorgangs verringert wird. Das Blockieren des Schleifenindex bei einer bestimmten Nummer hilft beim Speichern der Daten im Cache. Sie können wählen, ob Sie für den L2-Cache oder den L3-Cache optimieren möchten.

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

Pari Rajaram
quelle
-24

Aus vielen Gründen.

Erstens sind Fortran-Compiler stark optimiert, und die Sprache ermöglicht es ihnen, als solche zu sein. C und C ++ sind in Bezug auf die Array-Behandlung sehr locker (z. B. bei 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 optimierter, und der Compiler hat eine bessere Kontrolle darüber, was passiert, sodass 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 achten Sie darauf, 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 Ursache für Ineffizienz.

Drittens hängt es von der von Ihnen verwendeten blas-Implementierung ab. Einige Implementierungen werden möglicherweise in Assembler geschrieben und für den von Ihnen verwendeten Prozessor optimiert. Die Netlib-Version ist in fortran 77 geschrieben.

Außerdem führen Sie viele Vorgänge aus, von denen die meisten wiederholt und redundant sind. Alle diese Multiplikationen, um den Index zu erhalten, sind für die Leistung nachteilig. Ich weiß nicht wirklich, wie das in BLAS gemacht wird, aber es gibt viele Tricks, um teure Operationen zu verhindern.

Sie können Ihren Code beispielsweise 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 sicher, Sie werden etwas sparen.

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

Stefano Borini
quelle
36
Diese Antwort ist völlig falsch, sorry. BLAS-Implementierungen sind nicht in fortran geschrieben. Der leistungskritische Code wird in Assembly geschrieben, und die gängigsten werden heutzutage in C darüber geschrieben. Außerdem gibt BLAS die Zeilen- / Spaltenreihenfolge als Teil der Schnittstelle an, und Implementierungen können jede Kombination verarbeiten.
Andrew Tomazos
10
Ja, diese Antwort ist völlig falsch. Leider ist es voll von gesundem Menschenverstand, zB war die Behauptung BLAS wegen Fortran schneller. 20 (!) Positive Bewertungen zu haben ist eine schlechte Sache. Jetzt verbreitet sich dieser Unsinn aufgrund der Popularität von Stackoverflow sogar noch weiter!
Michael Lehn
12
Ich denke, Sie verwechseln die nicht optimierte Referenzimplementierung mit Produktionsimplementierungen. Die Referenzimplementierung dient nur zur Angabe der Schnittstelle und des Verhaltens der Bibliothek und wurde aus historischen Gründen in Fortran geschrieben. Es ist nicht für die Produktion bestimmt. In der Produktion verwenden die Mitarbeiter optimierte Implementierungen, die das gleiche Verhalten wie die Referenzimplementierung aufweisen. Ich habe die Interna von ATLAS (das Octave - Linux "MATLAB" unterstützt) studiert, von denen ich bestätigen kann, dass sie intern in C / ASM geschrieben sind. Die kommerziellen Implementierungen sind mit ziemlicher Sicherheit auch.
Andrew Tomazos
5
@KyleKanos: Ja, hier ist die Quelle von ATLAS: sourceforge.net/projects/math-atlas/files/Stable/3.10.1 Soweit ich weiß, handelt es sich um die am häufigsten verwendete Open-Source-Implementierung für tragbare BLAS. Es ist in C / ASM geschrieben. Hochleistungs-CPU-Hersteller wie Intel bieten auch BLAS-Implementierungen an, die speziell für ihre Chips optimiert wurden. Ich garantiere, dass Teile der Intel-Bibliothek auf niedriger Ebene in (duuh) x86-Assembly geschrieben sind, und ich bin mir ziemlich sicher, dass die Teile auf mittlerer Ebene in C oder C ++ geschrieben werden.
Andrew Tomazos
9
@ KyleKanos: Du bist verwirrt. Netlib BLAS ist die Referenzimplementierung. Die Referenzimplementierung ist viel langsamer als optimierte Implementierungen (siehe Leistungsvergleich ). Wenn jemand sagt, dass er Netlib BLAS in einem Cluster verwendet, bedeutet dies nicht, dass er tatsächlich die Netlib-Referenzimplementierung verwendet. Das wäre einfach albern. Es bedeutet nur, dass sie eine Bibliothek mit derselben Schnittstelle wie die netlib blas verwenden.
Andrew Tomazos