PDA

Archiv verlassen und diese Seite im Standarddesign anzeigen : Skalarprodukt/Vektor-Op. in SSE3 - Diskrepanzen


Marscel
2016-12-27, 03:12:34
Ich hab hier ein kleines privates Spaßprojekt um advanced Zeug auszuprobieren. Ein Teil davon ist es, Daten aus Bitmap-RGB Kanälen in Bitmap-JPEG-YCC (https://de.wikipedia.org/wiki/YCbCr-Farbmodell) umzuschreiben, um damit später noch weiterzuarbeiten.

Davon hab ich gleich drei Implementierungen. C, Rust und Ruby. Die ersten beiden produzieren bitgenau dieselben Daten, Ruby produziert nur eine Handvoll abweichender Bytes. Meine Herausforderung war nun, die C-Implementierung nochmal mit x86-Assembler umzuschreiben, und mangels AVX-CPU mit dem, was SSE3 noch hergibt.

Das Bitmap, das dabei gerade rauskommt, zeigt Zeichen von Überlaufoperationen. Jetzt dachte ich mir aus Debugging-Gründen, ich gucke mir erstmal einfach nur den Y-Kanal an (siehe Code) und kann an den Blauwerten schon ein wenig mehr sehen. Dabei kommt etwas raus (Anhang), das mir sagt, dass ich nicht grundsätzlich großen Mist produziert habe, aber man erkennt, dass das Blauspektrum bei der Assemblerversion etwas mehr Dynamik zeigt. Nun gut, lass ich mir das Histogramm der beiden Bilder anzeigen ... und bin gerade etwas ratlos, wo das Problem liegt. Das ist definitiv kein Rundungs-Modus-Problem mehr, sondern wohl strukturiert numerisch kaputt. Anhand der Werte bin ich aber unschlüssig.

Schaut mal bitte über die naive C-Version:

void transform_data(uint8_t *data, const size_t size) {
size_t t = size / 3;

for(size_t i = 0; i < t; i++) {
size_t s = i * 3;
int32_t b = data[s];
int32_t g = data[s+1];
int32_t r = data[s+2];

data[s] = 0.0f + ( 0.299f * r + 0.587f * g + 0.114f * b);
data[s+1] = 0.0f; // 128.0f + (-0.1687f * r + -0.3313f * g + 0.5f * b);
data[s+2] = 0.0f; // 128.0f + ( 0.5f * r + -0.4187f * g + -0.0813f * b);
}
}



Und dann die mit Assembler ausgestattete:


static float jc[3][4] = {
{ 0.299f, 0.587f, 0.114f, 0.0f },
{ 0.0f, 0.0f, 0.0f, 0.0f },
{ 0.0f, 0.0f, 0.0f, 0.0f },
/*{ -0.1687f, -0.3313f, 0.5f, 0.0f },
{ 0.5f, -0.4187f, -0.0813f, 0.0f }, */
};

static float jo[3] = { 0.0f, 128.0f, 128.0f };

#define ALIGNED_16B __attribute__((aligned(16)))

void transform_data(uint8_t *data, const size_t size) {
size_t t = size / 3;

for(size_t i = 0; i < t; i++) {
size_t s = i * 3;
__m128 vx ALIGNED_16B = { data[s], data[s+1], data[s+2], 0.0f };

for(size_t j = 0; j < 3; j++) {
__m128 jox ALIGNED_16B = { 0, 0, 0, jo[j] };
__m128 jcx ALIGNED_16B = *((__m128*)&jc[j]);
uint8_t r = 0;

asm(
"movaps %1, %%xmm0\n"
"movaps %2, %%xmm1\n"
"mulps %%xmm1, %%xmm0\n"
"haddps %%xmm0, %%xmm0\n"
"haddps %%xmm0, %%xmm0\n"
"addss %3, %%xmm0\n"
"cvtss2si %%xmm0, %%eax\n"
"movb %%al, %0\n"
: "=r"(r)
: "g"(jcx), "g"(vx), "g"(jox)
: "xmm0", "xmm1", "rax"
);

data[s+j] = r;
}
}
}



Seht jemand irgendwas, das ich falsch mache, in welcher Version auch immer? Die Assembler-Instruktionen sind viel nach Handbuch ausgesucht. Compiler ist GCC 6.2 auf Ubuntu 64bit.

Das Assembly der portablen C-Version (ohne -O Parameter) zeigt einen Haufen aneinandergeketteter SSE-Skalarwert-Operationen, nichts grundsätzlich Weltfremdes.

Könnte es wohl auch nochmal mit anderem Input versuchen, aber so aktuell find ich das zu interessant. ;)

Foobar2001
2016-12-27, 04:45:01
Ich wuerde dir ans Herz lesen Intrinsics fuer SSE zu verwenden und nicht den Inline-Assembler. Dann kannst du einfach Breakpoints setzen und dir die Zwischenwerte anschauen.

Ausserdem solltest du fuer maximale Performance SoA benutzen, also statt zu versuchen an einem Pixel rumzurechnen vier auf einmal laden und bearbeiten. Das sollte in dem Fall einfach sein.

Marscel
2016-12-27, 08:51:41
Intrinsics hatte ich schon, hatte aber, ich glaube, an der Intrinsic für cvtss2si ein wenig erfolglos probiert und dann erstmal sein gelassen. Ich werd es dann aber nochmal anschauen.

SoA würde doch aber im Wesentlichen nur noch den ungenutzten Höchstindexwert mitnehmen, oder? Dort, wo ich gerade 0*0 rechne. Wenn ich das andere Problem gelöst habe, schau ich mal weiter :)

Marscel
2016-12-27, 21:09:51
So, typischer Fall von Bäumen und Wäldern vor Augen: Die R- und B-Kanäle waren in der einen Version vertauscht.

Jetzt krieg tatsächlich nur noch Rundungsabweichungen, aber davon einen großen Haufen, die such ich schon noch.

Foobar2001
2016-12-28, 08:00:28
Intrinsics hatte ich schon, hatte aber, ich glaube, an der Intrinsic für cvtss2si ein wenig erfolglos probiert und dann erstmal sein gelassen. Ich werd es dann aber nochmal anschauen.

SoA würde doch aber im Wesentlichen nur noch den ungenutzten Höchstindexwert mitnehmen, oder? Dort, wo ich gerade 0*0 rechne. Wenn ich das andere Problem gelöst habe, schau ich mal weiter :)
Nein. Dein Datenlayout muss sich ändern. Getrennte arrays für rot, grün und blau. Dann lädst du die Daten für die Kanäle für jeweils n Pixel in SIMD-Register. Du rechnest also nicht mit RGB-Werten in einem Register, sondern mit vier R-Werten, vier G-Werten etc.

Was du da schreibst ist kein sinnvoler SIMD-Code, sorry. Die Hälfte davon sind Scalar-Ops und du nutzt mit den Vektor-Ops nicht mal die Einheiten aus. Und es wird auch nie mit AVX funktionieren.

gravitino
2017-01-14, 02:20:40
Herr Melax hat schönen Code zu SOA to AOS shuffling:

http://www.melax.com/normvec.pdf

Marscel
2017-01-14, 16:43:48
Danke, hab leider nicht allzu viel brauchbares seitdem geschafft.

Mein Ansatz für den nächsten Versuch - weil -O2 naiver Ansatz deutlich besser performt als mein händischer Code, in dem ich mittlerweile das for ausgerollt habe - wäre einmal komplett ohne Shuffle und horizontales Adding auszukommen.

Möglicherweise auch nicht die allerbeste Idee, aber der Spaß- und Lernfaktor steht im Vordergrund.