C++Guns – RoboBlog

18.09.2021

Vektorisieren leicht gemacht

Filed under: Allgemein — Tags: — Thomas @ 09:09

Als Beispiel sollen ein paar reduzierte Zeilen Code aus der Datei computeFluxes.h [1] aus dem Vplna-OP2 [2] Projekt dienen, welche auf unterschiedliche Arten vektorisiert werden sollen. Um so nah wie möglich am original Code zu bleiben, werden die Pointer Argumente der Funktionen, welche Arrays darstellen, wo es möglich ist nicht durch einen passerenden Typen ersetzt. Zu erst werden die GCC Vector Extensions benutzt, danach der Vektorisierer von GCC und zum Schluss OpenMP mit der SIMD Direktiven. Mit OpenMP ist es nicht nur möglich mehrere Tasks auf viele CPU Kerne laufen zu lassen, es lassen sich auch einzelne Schleifen vektorisieren. Untersucht wird der erzeugte ASM Code mit godbolt [3]

Compiliert wird mit GCC 11.2 -O2. Die stärke Optimierung -O3 schaltet zwar den Vektorisierer ein, aber erzeugt auch ASM Code der nicht mehr so leicht nachzuvollziehen ist. Der Vektorisierer kann manuell mit -ftree-vectorize [4] aktiviert werden. Daher wird der Quelltest mit dem Optionen -O2 -ftree-vectorize untersucht.

Zusammenfassung

Die GCC Vector Extensions bieten für dieses Beispiel die beste Art zu Vektorisieren, da der C++ Code am besten verständlich ist und ein hohen Grad an Vektorisierung ermöglicht. Es müssen keine Schleifen von Hand eingebaut oder mit OpenMP spezialisiert werden. Es ist Möglich den Code "einfach so hin zu schreiben". Es ist noch zu erwähnen, dass die "GCC Vector Extensions" Variante nur vier Zeilen C++ Code benötigt!

Der direkte Vergleich der Beispiele zeigt es deutlich. Verglichen wird die Anzahl relevanter C++ Code Zeilen (LOC), die totale Anzahl ASM Instruktionen, Anzahl SIMD Instruktionen die Vektorisiert werden könnten, Anzahl SIMD Instruktionen die vectorisiert wurden in Prozent. Sowie die Anzahl Sprung Anweisungen (Schleifen). Da bei jedem Sprung die Branch Prediction falsch liegen kann sind Sprünge als Performance Kritisch anzusehen.

Variante LOC Anzahl Instruktionen Anzahl SIMD Instruktionen Vektorisiert % Sprung Anweisung
GCC Vector Extensions 4 28 25 84.0 Nein
-O2 -ftree-vectorize Original Code 12 35 34 0.0 Nein
-O2 -ftree-vectorize Händisch Schleifen eingefügt 13 30 27 85.2 Nein
OpenMP 15 37 26 11.5 Ja

Es kann niemals eine 100%ige Vektorisierung stattfinden, da es wie indiesem Beispiel immer Stellen gibt, bei denen einzelne Skalare in unterschiedlichen Register stehen zusammen gezählt werden müssen. Der Unterschied zwischen 84% und 85.2% ist nicht relevant, da dieses Beispiel nur sehr wenige Instruktionen hat. Ein Prozent würde weniger als eine Instruktion entsprechen.

GCC Vector Extensions

GCC bietet mit Vector Extension die Möglichkeit SIMD Register direkt zu nutzen. Einfache Operationen wie Addition, Multiplikation sind darauf definiert. Damit lässt sich eine einfache Vector Klasse erstellen die es ermöglicht, lesbaren Code zu schreiben. Zusammen mit -ftree-vectorize ergibt sich die bestmögliche Abdeckung an Vektorizationen.

Das Scalarprodukt ist in Zeilen 5-16 makiert.

GCC -O2 -ftree-vectorize

float computeFluxes(const Vec4D& cellLeft, const float  *alphaleft,
                    const Vec2D& leftcellCenters, const Vec2D& edgeCenters,
                    const std::array<Vec2D,4>& leftGradient
) {
    Vec4D leftCellValues = cellLeft;
    
    Vec2D dxyl = edgeCenters - leftcellCenters;

    leftCellValues += alphaleft[0] * dot_product(dxyl, leftGradient);
 
    return leftCellValues[0] + leftCellValues[1] + leftCellValues[2] + leftCellValues[3]; 
}


computeFluxes(Vec4D const&, float const*, Vec2D const&, Vec2D const&, std::array<Vec2D, 4ul> const&):
        movq    (%rdx), %xmm1
        movq    (%rcx), %xmm0
        subps   %xmm1, %xmm0
        movups  (%r8), %xmm2
        movups  16(%r8), %xmm4
        movaps  %xmm2, %xmm1
        shufps  $221, %xmm4, %xmm2
        shufps  $136, %xmm4, %xmm1
        movaps  %xmm0, %xmm3
        shufps  $0xe5, %xmm0, %xmm0
        shufps  $0, %xmm0, %xmm0
        shufps  $0, %xmm3, %xmm3
        mulps   %xmm0, %xmm2
        mulps   %xmm3, %xmm1
        addps   %xmm2, %xmm1
        movss   (%rsi), %xmm0
        shufps  $0, %xmm0, %xmm0
        mulps   %xmm0, %xmm1
        addps   (%rdi), %xmm1
        movaps  %xmm1, %xmm0
        movaps  %xmm1, %xmm2
        shufps  $85, %xmm1, %xmm0
        addss   %xmm1, %xmm0
        unpckhps        %xmm1, %xmm2
        shufps  $255, %xmm1, %xmm1
        addss   %xmm2, %xmm0
        addss   %xmm1, %xmm0
        ret

Vektorisieren mit -O2 -ftree-vectorize Original Code

Das erste Listening zeigt auf der linken Seite die zu untersuchenden Codezeilen. Auf der rechten Seite die erzeugte ASM Ausabe. Die verwendeten ASM Befehle movss, addss, subss, mulss bedeuten: Kopieren, Addieren, Subtrahieren, Multiplizieren. Jeweils mit dem Kürzel "ss" welches man sich als "scalar single (precision)" merken kann. Die vektorisierte Variante wäre dementsprechend "ps" für "packed single (precision)".

Es fällt auf, dass keine "ps" Befehle vorkommen. Auch mit eingeschalteten Vektorisierer mittels -O3 oder -ftree-vectorize werden keine "ps" ASM Befehle erzeugt, da schlicht keine Schleifen im Code vorhanden sind, die vektorisiert werden können.

Die beispielhaft markierten Zeilen 11 und 12 im C++ Code entsprechend den ASM Befehlen 2-5.

GCC -O2

float computeFluxes(const float *cellLeft, const float *alphaleft,
                    const float *leftcellCenters, const float *edgeCenters,
                    const float *leftGradient
) {
    float leftCellValues[4];
    leftCellValues[0] = cellLeft[0];
    leftCellValues[1] = cellLeft[1];
    leftCellValues[2] = cellLeft[2];
    leftCellValues[3] = cellLeft[3];

    float dxl = (edgeCenters[0] - leftcellCenters[0]);
    float dyl = (edgeCenters[1] - leftcellCenters[1]);

    leftCellValues[0] += alphaleft[0] * ((dxl * leftGradient[0])+(dyl * leftGradient[1]));
    leftCellValues[1] += alphaleft[0] * ((dxl * leftGradient[2])+(dyl * leftGradient[3]));
    leftCellValues[2] += alphaleft[0] * ((dxl * leftGradient[4])+(dyl * leftGradient[5]));
    leftCellValues[3] += alphaleft[0] * ((dxl * leftGradient[6])+(dyl * leftGradient[7]));

    return leftCellValues[0] + leftCellValues[1] + leftCellValues[2] + leftCellValues[3]; 
}
computeFluxes(float const*, float const*, float const*, float const*, float const*):
        movss   (%rcx), %xmm1
        movss   4(%rcx), %xmm2
        subss   (%rdx), %xmm1
        subss   4(%rdx), %xmm2
        movss   (%r8), %xmm0
        movss   4(%r8), %xmm3
        movss   12(%r8), %xmm5
        movss   (%rsi), %xmm4
        mulss   %xmm2, %xmm3
        mulss   %xmm1, %xmm0
        mulss   %xmm2, %xmm5
        addss   %xmm3, %xmm0
        movss   8(%r8), %xmm3
        mulss   %xmm1, %xmm3
        mulss   %xmm4, %xmm0
        addss   (%rdi), %xmm0
        addss   %xmm5, %xmm3
        movss   20(%r8), %xmm5
        mulss   %xmm2, %xmm5
        mulss   %xmm4, %xmm3
        addss   4(%rdi), %xmm3
        mulss   28(%r8), %xmm2
        addss   %xmm3, %xmm0
        movss   16(%r8), %xmm3
        mulss   %xmm1, %xmm3
        mulss   24(%r8), %xmm1
        addss   %xmm5, %xmm3
        addss   %xmm2, %xmm1
        mulss   %xmm4, %xmm3
        addss   8(%rdi), %xmm3
        mulss   %xmm4, %xmm1
        addss   12(%rdi), %xmm1
        addss   %xmm3, %xmm0
        addss   %xmm1, %xmm0
        ret

-O2 -ftree-vectorize Händisch Schleifen eingefügt

Die markierten Zeilen im ersten Listening 1 werden durch eine Schleife ersetzt und mittels -ftree-vectorize vektorisiert. Dabei fiel auf, dass der Vektorisierer erst ab einer Schleifenlänge von 3 aktiv wird. Ob dies eine Compiler Einstellung oder eine allgemeine Beschränkung von SIMD ist, ist nicht bekannt.
Die nun markierten Zeilen 2-4, insbesondere Zeile 4, zeigt nun die "ps" Version der Subtraktion.

Die Zeilen 16-19 zeigen eine weitere Möglichkeit zum Vektorisieren.

GCC -O2 -ftree-vectorize

float computeFluxes(const float *cellLeft, const float *alphaleft,
                    const float *leftcellCenters, const float *edgeCenters,
                    const float *leftGradient
) {
    float leftCellValues[4];
    leftCellValues[0] = cellLeft[0];
    leftCellValues[1] = cellLeft[1];
    leftCellValues[2] = cellLeft[2];
    leftCellValues[3] = cellLeft[3];

    float dxyl[3];
    for(int i=0; i < 3; ++i) {
        dxyl[i] = (edgeCenters[i] - leftcellCenters[i]);
    }
    
    leftCellValues[0] += alphaleft[0] * ((dxyl[0] * leftGradient[0])+(dxyl[1] * leftGradient[1]));
    leftCellValues[1] += alphaleft[0] * ((dxyl[0] * leftGradient[2])+(dxyl[1] * leftGradient[3]));
    leftCellValues[2] += alphaleft[0] * ((dxyl[0] * leftGradient[4])+(dxyl[1] * leftGradient[5]));
    leftCellValues[3] += alphaleft[0] * ((dxyl[0] * leftGradient[6])+(dxyl[1] * leftGradient[7]));
 
    return leftCellValues[0] + leftCellValues[1] + leftCellValues[2] + leftCellValues[3]; 
}
computeFluxes(float const*, float const*, float const*, float const*, float const*):
        movq    (%rdx), %xmm0
        movq    (%rcx), %xmm2
        subps   %xmm0, %xmm2
        movss   4(%r8), %xmm3
        movss   12(%r8), %xmm5
        movss   (%r8), %xmm0
        movss   (%rsi), %xmm4
        movaps  %xmm2, %xmm1
        shufps  $0xe5, %xmm2, %xmm2
        mulss   %xmm1, %xmm0
        mulss   %xmm2, %xmm3
        mulss   %xmm2, %xmm5
        addss   %xmm3, %xmm0
        movss   8(%r8), %xmm3
        mulss   %xmm1, %xmm3
        mulss   %xmm4, %xmm0
        addss   (%rdi), %xmm0
        addss   %xmm5, %xmm3
        movss   20(%r8), %xmm5
        mulss   %xmm2, %xmm5
        mulss   %xmm4, %xmm3
        addss   4(%rdi), %xmm3
        mulss   28(%r8), %xmm2
        addss   %xmm3, %xmm0
        movss   16(%r8), %xmm3
        mulss   %xmm1, %xmm3
        mulss   24(%r8), %xmm1
        addss   %xmm5, %xmm3
        addss   %xmm2, %xmm1
        mulss   %xmm4, %xmm3
        addss   8(%rdi), %xmm3
        mulss   %xmm4, %xmm1
        addss   12(%rdi), %xmm1
        addss   %xmm3, %xmm0
        addss   %xmm1, %xmm0
        ret

Die Zeilen 16-19 im 2. Listening wurde durch eine Schleife ersetzt. Diese entsprechen jetzt Zeile 16-18 im 3. Listening. Die Schleife wurde vektorisiert, wie im ASM Output Zeile 5-22 zu sehen ist.
Es wurden einige shufps (shuffel packed single precision) Befehle erzeugt, um die Werte innerhalb eines SIMD Registers zu tauschen. Die erzeugte ASM Ausabe besteht nun hauptsächlich aus "ps" Befehlen und ist hiermit vektorisiert.

GCC -O2 -ftree-vectorize

float computeFluxes(const float *cellLeft, const float *alphaleft,
                    const float *leftcellCenters, const float *edgeCenters,
                    const float *leftGradient
) {
    float leftCellValues[4];
    leftCellValues[0] = cellLeft[0];
    leftCellValues[1] = cellLeft[1];
    leftCellValues[2] = cellLeft[2];
    leftCellValues[3] = cellLeft[3];

    float dxyl[3];
    for(int i=0; i < 3; ++i) {
        dxyl[i] = (edgeCenters[i] - leftcellCenters[i]);
    }

    for(int i=0; i < 4; ++i) {
        leftCellValues[i] += alphaleft[0] * ((dxyl[0] * leftGradient[i*2])+(dxyl[1] * leftGradient[i*2+1]));
    }
 
    return leftCellValues[0] + leftCellValues[1] + leftCellValues[2] + leftCellValues[3]; 
}
computeFluxes(float const*, float const*, float const*, float const*, float const*):
        movq    (%rdx), %xmm1
        movq    (%rcx), %xmm0
        subps   %xmm1, %xmm0
        movups  (%r8), %xmm2
        movups  16(%r8), %xmm4
        movaps  %xmm2, %xmm1
        shufps  $136, %xmm4, %xmm2
        shufps  $221, %xmm4, %xmm1
        movaps  %xmm0, %xmm3
        shufps  $0xe5, %xmm0, %xmm0
        shufps  $0, %xmm0, %xmm0
        mulps   %xmm0, %xmm1
        movaps  %xmm3, %xmm0
        shufps  $0, %xmm0, %xmm0
        mulps   %xmm0, %xmm2
        movss   (%rsi), %xmm0
        shufps  $0, %xmm0, %xmm0
        addps   %xmm2, %xmm1
        mulps   %xmm0, %xmm1
        movups  (%rdi), %xmm0
        addps   %xmm0, %xmm1
        movaps  %xmm1, %xmm0
        movaps  %xmm1, %xmm2
        shufps  $85, %xmm1, %xmm0
        addss   %xmm1, %xmm0
        unpckhps        %xmm1, %xmm2
        shufps  $255, %xmm1, %xmm1
        addss   %xmm2, %xmm0
        addss   %xmm1, %xmm0
        ret

Vektorisieren mit OpenMP

OpenMP besitzt seit Version 4.0 die SIMD Direktive welche ebenfalls "ps" ASM Befehle erzeugen kann. Listenin 4 zeigt die OpenMP Version mit entsprechenden ASM Output.

Die Vektorisierung der ersten Schleife funktioniert nun auch mit einer Länge von 2 statt 3, wie in Zeile 10 des ASM Output zu sehen ist.
Die Vektorisierung der zweiten Schleife schlägt allerdings fehl. Zu erkennen an den "ss" Instruktionen und dem Label .L2 welcher der Rücksprungort der for() Schleife ist.

GCC -O2 -fopenmp

float computeFluxes(const float *cellLeft, const float *alphaleft,
                    const float *leftcellCenters, const float *edgeCenters,
                    const float *leftGradient
) {
    float leftCellValues[4];
    leftCellValues[0] = cellLeft[0];
    leftCellValues[1] = cellLeft[1];
    leftCellValues[2] = cellLeft[2];
    leftCellValues[3] = cellLeft[3];

    float dxyl[2];
#pragma omp simd    
    for(int i=0; i < 2; ++i) {
        dxyl[i] = (edgeCenters[i] - leftcellCenters[i]);
    }

#pragma omp simd        
    for(int i=0; i < 4; ++i) {
        leftCellValues[i] += alphaleft[0] * ((dxyl[0] * leftGradient[i*2])+(dxyl[1] * leftGradient[i*2+1]));
    }
    
    return leftCellValues[0] + leftCellValues[1] + leftCellValues[2] + leftCellValues[3]; 
}
computeFluxes(float const*, float const*, float const*, float const*, float const*):
        movss   4(%rdi), %xmm0
        movq    (%rcx), %xmm1
        movq    %rdi, %rax
        movss   (%r8), %xmm2
        movss   (%rsi), %xmm3
        movss   %xmm0, -20(%rsp)
        movq    (%rdx), %xmm0
        movq    8(%rdi), %rdi
        subps   %xmm0, %xmm1
        movss   4(%r8), %xmm0
        movq    %rdi, -16(%rsp)
        movaps  %xmm1, %xmm4
        shufps  $0xe5, %xmm1, %xmm1
        mulss   %xmm4, %xmm2
        mulss   %xmm1, %xmm0
        addss   %xmm2, %xmm0
        mulss   %xmm3, %xmm0
        addss   (%rax), %xmm0
        movl    $1, %eax
        movss   %xmm0, -24(%rsp)
.L2:
        movss   (%r8,%rax,8), %xmm0
        movss   4(%r8,%rax,8), %xmm2
        mulss   %xmm4, %xmm0
        mulss   %xmm1, %xmm2
        addss   %xmm2, %xmm0
        mulss   %xmm3, %xmm0
        addss   -24(%rsp,%rax,4), %xmm0
        movss   %xmm0, -24(%rsp,%rax,4)
        addq    $1, %rax
        cmpq    $4, %rax
        jne     .L2
        movss   -24(%rsp), %xmm0
        addss   -20(%rsp), %xmm0
        addss   -16(%rsp), %xmm0
        addss   -12(%rsp), %xmm0
        ret

Anhang

Die für Vector Klasse für die GCC Vector Extensions Variante.

#include <array>

struct Vec2D {
    typedef float vec_type __attribute__ ((vector_size (2*sizeof(float))));
    vec_type _data;

    const float& operator[](int i) const {
        return _data[i];
    }   
};

inline Vec2D operator-(const Vec2D& lhs, const Vec2D& rhs) {
    return Vec2D{lhs._data - rhs._data};
}

struct Vec4D {
    typedef float vec_type __attribute__ ((vector_size (4*sizeof(float))));
    vec_type _data;

    float& operator[](int i) {
        return _data[i];
    }
};

inline Vec4D& operator+=(Vec4D& lhs, const Vec4D& rhs) {
    lhs._data += rhs._data;
    return lhs;
}

inline Vec4D operator*(const float lhs, const Vec4D& rhs) {
    return Vec4D{lhs * rhs._data};
}

inline Vec4D dot_product(const Vec2D& v1, const std::array<Vec2D,4>& v2) {
    Vec4D result;
    for(int i=0; i < 4; ++i) {
        result[i] = ((v1[0] * v2[i][0]) + (v1[1] * v2[i][1]));
    }
    return result;
}

[1] https://github.com/reguly/volna/blob/master/sp/computeFluxes.h
[2] https://gmd.copernicus.org/articles/11/4621/2018/
[3] https://godbolt.org/
[4] https://gcc.gnu.org/onlinedocs/gcc-11.2.0/gcc/Optimize-Options.html#index-ftree-vectorize
[5] https://gcc.gnu.org/onlinedocs/gcc/Vector-Extensions.html

No Comments

No comments yet.

RSS feed for comments on this post.

Sorry, the comment form is closed at this time.

Powered by WordPress