C++Guns – RoboBlog blogging the bot

Im dritten Teil geht es um die Parallelisierung mit OpenMP. Auf dem Arbeitsrechner werden 1000 zufällig gefüllte Strings in 71.6 sec. Der Algorithmus ist trivial über die äusserste Schleife zu parallelisieren, da ein String unabhänig der anderen Berechnet werden kann.

Um eine schnelle Übersicht über die systime zu bekommen, wird das Programm mittels dem Hilfsprogramm 'time' gestartet. Die systime beträgt 0m0.004s dieses Wert müssen wir nach der Parallelisierung optimalerweise wieder erreichen.

Ein einfaches

#pragma omp parallel for
    for(int i=0; i < strings.size(); ++i) {

und schon brummen alle 8 Threads des AMD FX 8350 los. Die Zeit sagt 12.7 sec mit 0m0.008s systime. Das entspricht einem Speedup von 5.6 70%. Ein normaler Wert für den ersten Versuch. Um ein besseres Gefühl zu bekommen wie der Algorithmus skaliert, wird nochmal mit nur 4 Threads getestet. Es ist stark CPU abhängig ob ob auch alle Kerne einer CPU mit voller Leistung brummen können.

$ OMP_NUM_THREADS=4  time ./LevenshteinTerror

Mit 4 Threads ergeben sich 18.4 sec und ein Speedup von 3.89 97%. Besser geht es eigentlich garnicht.

Damit können wir gleich zu einer dickeren CPU wechseln: Intel Xeon E5-2680 mit 32 Kerne. Hier die Ergebnisse:

Threads Time [sec] Speedup Eff [%]
1       67.7          -      -
2       33.9        2      100 
4       17.0        3.98   99.5
8        9.22       7.34   91.8
16       4.94      13.7    85.6
32       4.54      14.9    46.6

Bis zu der Hälfte der Threads skaliert der Algorithmus gut, danach gibt es keinen Weiteren Speedup. Es werden in diesem Beispiel auch keine Unmengen an Daten verarbeitet, die CPU ist der Flaschenhals, nicht der RAM Bus.

Zum Abschluss noch eine Messung mit 16 und 32 Threads mit der 10-fachen Menge an Strings.

Threads Time [sec] Speedup Eff [%]
1       6709          -      -
16      475         14.1    88.1
32      447         15.0    46.9

Zwei Stunden v.s. 8 Minuten.

In der Regel wird die min() Funktion mittels Verzweigung implementiert. Ohne zusätzliche Optimierung führt dies zu einem Sprung im Code und damit zu Branchs Mispredicts. Dies ist gerade beim HPC Code ziemlich schlecht, da im Fall einer falscher Sprungvorhersage die CPU Pipeline wieder zurückgesetzt werden muss.

Es gibt aber Alternativen: Conditional Move CPU Instruktionen. Je nach Bedingung wird ein Wert in ein anderes Register kopiert, oder eben nicht. Da steckt kein Sprung im Code dahinter, das wird alles auf CPU Ebene realisiert. So wird dieser C++ Code zu folgendem Assember ersetzt:

int func(int x, int y) {
    int a = y;
    if(x < y) {
        a = x;
    }
    return a;
}

func(int, int):
  cmp esi, edi        x < y ?
  mov eax, edi        a = y
  cmovle eax, esi     Wenn x<y dann a = x
  ret                 return a

Es ist aber garnicht nötig so komplizierten C++ Code, oder gar inline Assember zu schreiben, der Compiler erkennt das Muster automatisch. So führt die Benutzung von std::min(x,y) zu exakt dem selben Assember Code.

Das ganze Funktioniert bei meinen Tests mit Optimierung O1 und O2, aber nicht mehr mit O3. Da bei dieser Optimierung wohl zuviel Schleifen Optimierungs Magie zum Einsatz kommt. Die verantwortlichen Flags müssten -fif-conversion -fif-conversion2 sein, welche bei jeder Optimierungsstufe aktiv sein.

-fif-conversion

Attempt to transform conditional jumps into branch-less equivalents. This includes use of conditional moves, min, max, set flags and abs instructions, and some tricks doable by standard arithmetics. The use of conditional execution on chips where it is available is controlled by -fif-conversion2.

Enabled at levels -O, -O2, -O3, -Os.
-fif-conversion2

Use conditional execution (where available) to transform conditional jumps into branch-less equivalents.

Enabled at levels -O, -O2, -O3, -Os.

Selbstverständlich funktioniert das auch für Floatingpoint Typen, da reich auch eine einzelne Assembler Anweisung: std::min(x,y)

Von der Compiler Optimierung O3 habe ich mir keine weiteren Geschwindigkeitsgewinn erhofft, aber dass das Programm doppelt so langsam! wurde, da stand mir dann doch der Mund offen ;)
Das ist btw. auch das erste mal in meinem Programmiererleben, dass O3 die Laufzeit wirklich deutlich verschlechtert. Bisher hatte ich nur davon gehört. Ein schneller Blick in den Profiler zeigt auch warum. Gab es bei O2 und der innersten Schleife ein paar 100k Branch Mispredicts, waren es bei O3 auf einmal ein paar Millionen! Die Stelle im Code war auch im ersten Teil schon etwas verwirrend. Anscheinend gibt es eine Code Optimierung, welche die Verzweigung in der min() Funktion entfernt, ganz ohne Bithacks.
Nun war es schon spät am Abend, aber die Sache musste unbedingt untersucht werden! Also habe ich dem vom Compiler erzeugten Assembler Code studiert und konnte die Sache auf ein Minimalbeispiel runter brechen.

Betrachtet werden muss nur die innerste Schleife, da nur hier die min() Funktion vorkommt.

#include <string>
#include <vector>

void func(size_t N) {
    std::vector<int> v0(N+1), v1(N+1);

    for(size_t i=0; i < N; ++i) {   
        v1[i+1] = std::min(v0[i], v0[i+1]);
    }
}

Es wird paarweise das Minimum gesucht. Wird der Eingangs Array mit Zufallszahlen gefüllt, so muss die spekulative Ausführung der min() Funktion, also Branch Mispredict 50% betragen. Der Compiler erkennt dies irgendwie, und produziert folgenden Assembler Code ohne spekulative min() Funktion.

  mov ecx, DWORD PTR [rbx]             ecx = v0[0]
  xor edx, edx                         i = 0
.L10:                                  for()
  add rdx, 1                           ++i
  mov esi, DWORD PTR [rbx+rdx*4]       esi = v0[i+1]
  cmp esi, ecx                         v0[i+1] <= ecx ?
  cmovle ecx, esi                      Wenn true: ecx = esi  Hier passiert die Magie.
  cmp rdx, rbp                         i < N
  mov DWORD PTR [rdi+rdx*4], ecx       v1[i+1] = ecx
  mov ecx, esi                         ecx = v0[i+1]
  jb .L10

Der Code ist nicht weiter schwer. Ich habe eine grobe Übersetzung in die rechte Spalte geschrieben. In der ersten Zeile wird das erste Element von v0 in das Register ecx kopiert. Und in der zweiten Zeile das Register edx, also die Laufvariable i, auf 0 gesetzt. Danach beginnt in Zeile 3 der Schleifenkörper.
Variable i wird um eins erhöht und Zeile fünf erfolgt der lesende Zugriff auf v0[i+1] und der Wert wird in das Register esi geladen.
In den Register esi und ecx stehen jetzt die beiden Werte die verglichen werden sollen. Dabei steht in ecx immer der Wert der Vorherigen Schleifeniteration. Damit wird nur ein RAM Zugiff pro Iteration benötigt.
Verglichen werden die Werte in Zeile 6 und das Ergebnis des Vergleichs wird in Zeile 7 ausgewertet. Nur wenn der neue Wert an Stelle i+1 kleiner oder gleich als der alte ist, wird mit dem neuen Wert der alte überschrieben.
Das Ergebnis wird in Zeile 9 nach v1[i+1] kopiert und in Zeile 10 wird ecx für die nächste Iteration vorbereitet. Zeile 8 und 11 entspricht noch dem Schleifenkopf mit dem Rücksprung zum Label L10.

Die Magie passiert bei der Assembler Instruktion cmovle. Hier wird ein Kopiervorgang nur dann ausgeführt, wenn das Ergebnis des Vergleichs eine Zeile weiter oben ergeben hat, dass die Zahl kleiner ist. Es gibt also keine Verzweigung im Programmablauf und keine spekulative Ausführung. Daher auch kein Branch Mispredict. Der Code kann mit voller CPU Power laufen. Besser geht es nur noch mit SSE und prefetch.

Dies erklärt also warum die min() Funktion keine Branch Mispredicts erzeugt. Aber warum läuft der O3 Optimierte Code nicht schneller. Auch die Assembler Code hab ich untersucht und der Grund ist unter anderem tree vectorize und andere Schleifen Optimierungen die der Compiler vornimmt. Aber die Schleife ist so super einfach, dass jegliche Optimierung nur eine Verschlechterung verursacht. Auch Loop unroll, welches per default nicht mit O3 aktiviert wird, wird die Laufzeit nicht verbessern. Da der Compiler nicht wissen kann, wie oft die innerste Schleife laufen wird. Daher müssen alle Möglichkeiten betrachtet werden und dies führt zu ein paar mehr Jumps im Code. Und schon führt die verlangsamte spekulative Code Ausführung wieder zu.

Was haben wir gelernt? Optimierung O2 ist gut. O3 bringt in der Summe genauso viel wie sie wieder zerstört. Bei HPC macht es keinen Sinn immer mit O3 zu kompilieren. Vielmehr muss der Algorithmus und die Arbeitsweise der Ziel CPU bis ins Detail verstanden werden. Danach kann die passende Compiler Optimierung ausgewählt werden.

Getestet mit GNU GCC 7.3

Funktionen können Funktionen als Rückgabe Type/Wert haben, solange sie die selbe Signatur vorweisen. Also z.B. alles Binary Funktionen mit Argumenten von Type double.

#include <iostream>
#include <vector>
#include <functional>

double add(const double lhs, const double rhs) {
    return lhs+rhs;
}

double mul(const double lhs, const double rhs) {
    return lhs*rhs;
}

// does not work
double PI() {
    return 3.14;
}

enum class Operator { add, mul, PI };

auto parse(Operator op) {
    switch(op) {
    case Operator::add: {
        return add;
    }
    case Operator::mul: {
        return mul;
    }
    case Operator::PI: {
        // error: inconsistent deduction for auto return type: ‘double (*)(double, double)’ and then ‘double (*)()’
        // return PI;
    }
    }
}

int main() {
    auto l = parse(Operator::add);
    auto result = l(1,2);
    std::cout << result << "\n";
}

Wen hat es nicht schon einmal gestört, dass std::sort() so transparent arbeitet? Manchmal muss man wissen welche Elemente beim Sortieren wohin verschoben werden. Sozusagen eine Permutation.
Manchmal ist es möglich, den zu sortierenden Daten eine ID mitzugeben. Aber oft nicht. Zum Beispiel arbeite ich mit 2D Punkten. Diese liegen natürlich schön dicht gepackt im std::vector, und die ID ist implizit durch die Position im vector gegeben. Durch das Sortieren würde sich auch die ID ändern, was aber nicht geht, da sie mit andern Daten noch verknüpft ist. Und die ID explicit abspeichern ist auch nicht gut, da dann die x/y Koordinaten vom Punkt im Speicher nicht mehr optimal ausgerichtet sind. Das Speicherlayout ändert sich zu meinen Ungunsten.

Zumindest in meinem Anwendungsfall brauche ich die sortierten Daten auch nur kurzzeitig, und die original Daten sollen unverändert bleiben. Es würde teilweise auch vollkommen langen, wenn die sort Funktion nur ein vector mit neuen IDs in der sortierten Reihenfolge zurück gibt.

Mein erster Ansatz ist, eine sort() Funktion zu schreiben, welche als zusätzliches Argument einen Visitor akzeptiert. Diese Funktion ruft die eigentlich std::sort() Funktion auf, so dass alles vom Standard noch eingehalten wird und sich auch nichts an der asymptotischen Laufzeit ändert.
Um eine ID zu erhalten, werte ich die Adresse der Elemente im std::vector aus. Damit die original Daten (erst einmal) in ihrer ursprünglichen Reihenfolge bleiben, wird ein zweiter std::vector mit Referenz Wrappern erstellt. Dies kostet pro Element den Speicherplatz einer Referenz, also 8byte. Dies ist, denke ich, nie ein Problem.

Nach der eigentlichen Sortierung mit std::sort ist so feststellbar, wie die Elemente im vector getauscht wurden. Damit wird das Visitor Objekt aufgerufen, so dass der Anweder den Permutationsvektor erstellen kann. Und die eigentlichen Daten werden durch swap sortiert.

Hier der Code der ersten Implementierung:

template<typename T, typename Allocator, typename Visitor>
inline void sort(std::vector<T,Allocator>& data,  Visitor&& vis) {
    std::vector<std::reference_wrapper<const T>> sorted;
    sorted.reserve(data.size());
    for(const auto& x : data) {
        sorted.push_back(std::cref(x));
    }

    std::sort(sorted.begin(), sorted.end());
    // ptrdiff_t is the signed equivalent to size_t
    for(size_t i=0; i < sorted.size(); ++i) {
        // get the address of an item in data, not sorted.
        const size_t idx1 = std::addressof(sorted[i].get()) - data.data();
        if(i < idx1) {
            std::swap(data[i], data[idx1]);
            vis(i, idx1);
        }
    }
}

//////

    struct Visitor {
        void operator()(size_t i, size_t j) {
            std::cout << "Swap idx " << i << " with idx " << j << "\n";
        }
    };

    std::vector<int> v {3,2,1,0};
    cpl::sort(v, Visitor());
    QVERIFY(std::is_sorted(v.begin(), v.end()));

Swap idx 0 with idx 3
Swap idx 1 with idx 2