C++Guns – RoboBlog blogging the bot

Part1: print std::array with std:integer_sequence
Part 2: convert tuple to parameter pack
Part 3: print std::array with std::apply and fold
Part 4: fold over std::tuple und erzeugten Assembler Code
Part 5: fold over std::tuple of std::vector of Types ...
Part 6: apply generic lambda to tuple
Part 7: Play with std::tuple and std::apply

Wird an der Zeit, dass ich das auch mal versuche:

#include <array>
#include <iostream>

template<size_t N, std::size_t... Ints>
void print_impl(const std::array<double,N>& arr, std::index_sequence<Ints...>) {
    ((std::cout << std::get<Ints>(arr)),...);
}

template<size_t N>
void print(const std::array<double,N>& arr) {
    print_impl(arr, std::make_index_sequence<N>{});
}

void func() {
    std::array<double,4> arr{1,2,3,4};
    print(arr);
}

Das selbe Beispiel nochmal mit Pascal Code. Für ein halbwegs vernünftiges Ergebnis musste ich Optimierung O2 auswählen und von Hand inline einschalten. Dead Code elimination funktioniert nicht. Die wesentliche Subroutine hat fünf Subtraktionen, zwei Multiplikationen, elf MOVs, ein LEA und der Stack wird genutzt. Auch zum übergeben der Funktionsparameter. Zur Erinnerung: Wir brauchen fünf Subtraktionen, zwei Multiplikationen und vier explizite Kopier-Befehle (MOV) in C++.
Besser wird es nicht.

unit output;
interface
implementation

type
  Point2D = record
    x: double;
    y: double;
  end;

type 
  Line2D = record
    pt1: Point2D;
    pt2: Point2D;
  end;

Operator -(p1: Point2D; p2: Point2D) res: Point2D inline ;  
begin  
  res.x := p1.x-p2.x;  
  res.y := p1.y-p2.y;  
end; 

  function func(line1: Line2D; line2: Line2D) : double;
var   
    denominator: double; 
    a: Point2D;
    b: Point2D;
begin
    a := line1.pt2 - line1.pt1;
    b := line2.pt1 - line2.pt2;
    denominator := a.y * b.x - a.x * b.y;
    exit(denominator)
end;
  
end.

func(line2d,line2d):
  pushq %rbp
  movq %rsp,%rbp
  leaq -48(%rsp),%rsp
  movsd 32(%rbp),%xmm0
  subsd 16(%rbp),%xmm0
  movsd %xmm0,-16(%rbp)
  movsd 40(%rbp),%xmm0
  subsd 24(%rbp),%xmm0
  movsd %xmm0,-8(%rbp)
  movsd 48(%rbp),%xmm0
  subsd 64(%rbp),%xmm0
  movsd %xmm0,-32(%rbp)
  movsd 56(%rbp),%xmm0
  subsd 72(%rbp),%xmm0
  movsd %xmm0,-24(%rbp)
  movsd -8(%rbp),%xmm0
  mulsd -32(%rbp),%xmm0
  movsd -16(%rbp),%xmm1
  mulsd -24(%rbp),%xmm1
  subsd %xmm1,%xmm0
  leave
  ret

Schaut euch mal folgenden Anfang einer Line2D intersect Funktion an. Erstaunlich wie sauber und perfekt der Assember Code aus meinem C++ Code generiert wird. Da ist absolut kein Overhead zu erkennen. Keine Funktionsaufrufe, kein unnötiges Laden und zwischenspeichern von Daten. Keine temporäre Objekte auf dem Stack. Überhaupt keine Stack Nutzung. Es gibt im Code fünf Subtraktionen und zwei Multiplikationen. Das resultiert im Assember Code mit fünf Subtraktionen, zwei Multiplikationen und vier explizite Kopier-Befehle für doubles. Besser geht es doch gar nicht mehr. Und das alles ist nur mit Optimerungslevel O1 compiliert. Nicht O2, nicht fanzy O3, nein, nur O1. GNU GCC alle gängigen Versionen. Clang kann es natürlich nicht mit O1, nur mit O2, aber dann sieht der ASM Code nicht mehr so schön symmetrisch aus ;) Bei Intel funktioniert es auch mit O1, die ASM Code Anordnung ist etwas anders. Microsoft... nein.
Und die wichtigsten Zeilen vom C++ Code sind doch auch ausdrucksstark.

#include <array>

struct Point2D : public std::array<double,2> {
    inline const auto& x() const { return at(0); }
    inline const auto& y() const { return at(1); }
};

inline const Point2D operator-(const Point2D& p1, const Point2D& p2) {
    return Point2D{p1.x()-p2.x(), p1.y()-p2.y()};
}

struct Line2D : public std::array<Point2D, 2> {
    inline const Point2D& p1() const { return at(0); }
    inline const Point2D& p2() const { return at(1); }
};

auto func(const Line2D& line1, const Line2D& line2) {
    Point2D a = line1.p2() - line1.p1();
    Point2D b = line2.p1() - line2.p2();
    const auto denominator = a.y() * b.x() - a.x() * b.y();
    return denominator;
}

func(Line2D const& line1, Line2D const& line2):
  movsd (%rsi),   %xmm0 # line2 x1
  subsd 16(%rsi), %xmm0 # line2 x2
  movsd 24(%rdi), %xmm1 # line1 y2
  subsd 8(%rdi),  %xmm1 # line1 y1
  mulsd %xmm1,    %xmm0
  movsd 8(%rsi),  %xmm1 # line2 y1
  subsd 24(%rsi), %xmm1 # line2 y2
  movsd 16(%rdi), %xmm2 # line1 x2
  subsd (%rdi),   %xmm2 # line1 x1
  mulsd %xmm2,    %xmm1
  subsd %xmm1,    %xmm0
  ret 

Übrigends bekommt man den selben Assember Code auch mit float statt double. Und für int bleibt die Struktur auch die selbe, es werden nur die normalen Register benutzt. Nur bei long double kommen ein paar Kopierbefehle dazu, da die normalen floating point Register der CPU genutzt werden, statt SSE. Da SSE aber 128bit floating point Zahlen verarbeiten kann, muss ich wohl die CPU beim Compiler angeben.

// Nachtrag
Einen setz ich noch drauf: SSE mit GNU gcc Vector Instructions Extension. Point2D mit double als Type lässt sich wunderbar per SSE verarbeiten. So können zwei Subtraktionen/Multiplikationen parallel ausgeführt werden. Damit lässt sich extrem kompakte (und super effizienter) ASM Code generieren.

#include <smmintrin.h>

using v2sd = double __attribute__ ((vector_size (16)));
struct  Point2D : public std::tuple<v2sd> {
....
}

auto func(const Line2D& line1, const Line2D& line2) {
    const Point2D a = line1.p2() - line1.p1();
    const Point2D b = (line2.p1() - line2.p2());
    const auto denominator = a.y()*b.x() - a.x()*b.y();
    return denominator;
}

Nur noch drei Substraktionen. Leider kommen so unpack Anweisungen dazu.

func(Line2D const&, Line2D const&):
  movapd 16(%rdi), %xmm1
  subpd (%rdi), %xmm1
  movapd (%rsi), %xmm2
  subpd 16(%rsi), %xmm2
  movsd %xmm2, %xmm0
  movapd %xmm1, %xmm4
  unpckhpd %xmm4, %xmm4
  mulsd %xmm4, %xmm0
  unpckhpd %xmm2, %xmm2
  mulsd %xmm2, %xmm1
  subsd %xmm1, %xmm0

Hier mal ein Versuch ein hübscheren tokenizer zu bauen:

/// \todo TH use string_view
std::list<std::string> tokenize(const std::string& str) {
    std::list<std::string> tokens;
    std::string token;

    auto add = [&tokens](std::string& token) {
        if(not token.empty()) {
            tokens.emplace_back(std::move(token));
        }
    };

    for(const char c : str) {
        if(c == ' ') {
            add(token);
        } else if(c == '(' or c == ')') {
            add (token);
            tokens.emplace_back(std::string{c});
        } else {
            token.push_back(c);
        }
    }

    add(token);

    return tokens;
}

Ist auf alle Fälle Millionen mal besser als Dieses hier:
https://gist.github.com/ofan/721464

std::list<std::string> tokenize(const std::string & str)
{
    std::list<std::string> tokens;
    const char * s = str.c_str();
    while (*s) {
        while (*s == ' ')
            ++s;
        if (*s == '(' || *s == ')')
            tokens.push_back(*s++ == '(' ? "(" : ")");
        else {
            const char * t = s;
            while (*t && *t != ' ' && *t != '(' && *t != ')')
                ++t;
            tokens.push_back(std::string(s, t));
            s = t;
        }
    }
    return tokens;
}

hf

Also das Video 15 Sorting Algorithms in 6 Minutes ist ja mal sau genial. Am besten hört sich der Radix Sort an ;) Sehr interessant diese vielen Sortieralgorithmen. Muss man sich unbedingt mal anschaun.
Aber der Code für dieses Projekt ist natürlich wieder .... ;)

In dem Video werden die Anzahl der Vergleiche und Array Zugriffe gezählt. Ich frage mich ob so etwas nicht einfacher zu realisieren ist. Da std::sort eine Vergleichsfunktion als Argument annimmt, sollte es sehr einfach sein, die Anzahl der Vergleiche zu zählen. Die anderen beiden Argumente sind Iteratoren. Also müsste man den indirection operator überladen. Das sollte machbar sein.

Hier mein Versuch. Ständig die Iterator Logik bei so Probleme neu zu schreiben ist nicht schön, aber zumindest kann man mit dieser Technik JEDEN Algorithmus untersuchen der Iteratoren akzeptiert. Besser als das gekrüppel was der andere da macht.

#include <vector>
#include <algorithm>
#include <iostream>

struct CounterIteratorAdapter : public std::vector<int>::iterator {
public:
    using Base = std::vector<int>::iterator;

    int& operator*() {
        ++counter;
        return Base::operator *();
    }

    auto operator+(int rhs) const {
        return CounterIteratorAdapter{Base::operator +(rhs)};
    }

    auto operator-(int rhs) const {
        return CounterIteratorAdapter{Base::operator -(rhs)};
    }

    static int counter;
};

int CounterIteratorAdapter::counter = 0;

int main() {
    std::vector<int> data(100);
    for(int i=0; i < data.size(); ++i) {
        data[i] = data.size()-i;
    }

    int comparisons = 0;
    auto compare = [&comparisons](const int lhs, const int rhs) { ++comparisons; return lhs < rhs;};

    std::sort(CounterIteratorAdapter{data.begin()}, CounterIteratorAdapter{data.end()}, compare);

    std::cout << comparisons << " comparisons "<< CounterIteratorAdapter::counter << " array accesses\n";
}

Random:
830 comparisons 2223 array accesses
Reserve:
526 comparisons 1355 array accesses
Sorted:
629 comparisons 1474 array accesses

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