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. Zu erst wird der Vektorisierer von GCC benutzt, danach OpenMP mit der SIMD Direktive. Und zum Schluss GCC Vector Extensions. 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.
Vektorisieren mit -O3 / -ftree-vectorize
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 -O oder -ftree-vectorize werden keine "ps" ASM Befehle erzeugt, da schlicht keine Schleifen im Code vorhanden sind, die vektorisiert werden können.
Die 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
|
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
|
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 Skalar Produkt ist in den Codeausschnitten makiert.
GCC -O2 -ftree-vectorize
#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;
}
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
|
[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