Obsah
1. Instrukční sady SIMD a automatické vektorizace prováděné překladačem GCC (2)
2. První příklad: součet všech prvků v krátkém jednorozměrném poli se čtyřmi prvky
3. Výsledek překladu se zákazem i povolením automatické vektorizace
4. Automatická vektorizace po specifikaci přepínače -ffast-math
6. Druhý příklad: součet všech prvků v krátkém jednorozměrném poli s osmi prvky
7. Výsledné kódy vygenerované překladačem: „skalární“ i „vektorová“ varianta
8. Třetí příklad: součet pole se 100 prvky
9. Výsledné kódy vygenerované překladačem: „skalární“ i „vektorová“ varianta
10. Čtvrtý příklad: vyhledání největšího prvku v poli
11. Nevektorizovaná a vektorizovaná podoba výsledného kódu
12. Vliv velikosti pole na výsledný vygenerovaný kód
13. Jeden z nejdůležitějších algoritmů současnosti
14. Pátý příklad: skalární součin dvou krátkých polí se čtyřmi prvky
15. Nevektorizovaná a vektorizovaná podoba výsledného kódu
16. Další možné optimalizace výpočtů skalárního součinu
17. Vysvětlení postupu výpočtu
18. Repositář s demonstračními příklady
19. Seznam všech předchozích částí tohoto seriálu a článků o SIMD instrukcích
1. Instrukční sady SIMD a automatické vektorizace prováděné překladačem GCC (2)
V úvodním článku o technologii automatické vektorizace prováděné při překladu zdrojových kódů do strojového kódu (nebo do assembleru) jsme si ukázali základní techniky podporované překladačem GCC v případě, že se nějaká operace provádí se všemi prvky pole (resp. v terminologii SSE vektoru). Jednalo se například o operace vynulování všech prvků pole, přičtení hodnoty ke všem prvkům pole nebo součet odpovídajících si prvků dvou polí.
Ovšem technologie automatické vektorizace dokáže v některých případech korektně zpracovat i smyčky, ve kterých se provádí „redukce“ – součet všech prvků v poli, nalezení nejmenšího nebo naopak největšího prvku v poli, ale i výpočet skalárního součinu atd. A právě na způsob překladu takových zdrojových kódů se zaměříme dnes. Pro jednoduchost se stále budeme omezovat na instrukce SSE.
2. První příklad: součet všech prvků v krátkém jednorozměrném poli se čtyřmi prvky
V dnešním prvním demonstračním příkladu si otestujeme, jakým způsobem je možné vektorizovat zdrojový kód, ve kterém se sice zpracovávají všechny prvky nějakého pole (resp. vektoru), ovšem výsledkem tohoto zpracování bude jediná skalární hodnota vzniklá zkombinováním hodnot všech prvků. Nejjednodušší takový příklad počítá sumu všech prvků pole a realizovat je ho možné velmi snadno. Začneme tím nejmenším možným případem, který je ještě možné vektorizovat – součtem všech (čtyř) prvků čtyřprvkového pole:
float array_sum_4(float *a) {
#define SIZE 4
int i;
float result = 0.0;
for (i=0; i<SIZE; i++) {
result += a[i];
}
return result;
}
Překladač tedy musí umět rozeznat, že se provádí nějaká forma algoritmu typu reduce (někdy se setkáme s označením reduction operator, které je přesnější).
3. Výsledek překladu se zákazem i povolením automatické vektorizace
Podívejme se nyní, jak bude vypadat výsledek překladu céčkovského kódu do assembleru (resp. později do strojového kódu). První kód byl překladačem vytvořen ve chvíli, kdy byly zakázány veškeré optimalizace (neuvedli jsme totiž přepínač -O2). Povšimněte si, že je realizována počítaná programová smyčka, ve které se volají „skalární“ instrukce MOVSS a ADDSS. Dokonce ani nedošlo k pokusu o rozbalení této smyčky (což je dobře, když jsme optimalizace zakázali):
array_sum:
push rbp
mov rbp, rsp
mov QWORD PTR [rbp-24], rdi
pxor xmm0, xmm0
movss DWORD PTR [rbp-8], xmm0
mov DWORD PTR [rbp-4], 0
jmp .L2
.L3:
mov eax, DWORD PTR [rbp-4]
cdqe
lea rdx, [0+rax*4]
mov rax, QWORD PTR [rbp-24]
add rax, rdx
movss xmm0, DWORD PTR [rax]
movss xmm1, DWORD PTR [rbp-8]
addss xmm0, xmm1
movss DWORD PTR [rbp-8], xmm0
add DWORD PTR [rbp-4], 1
.L2:
cmp DWORD PTR [rbp-4], 3
jle .L3
movss xmm0, DWORD PTR [rbp-8]
pop rbp
ret
Po povolení optimalizací, ovšem bez povolení automatické vektorizace, je výsledek mnohem lepší. Namísto programové smyčky se pouze provede součet všech prvků pole. Zajímavé je, že se stále začíná s nulovou hodnotou uloženou do registru XMM0, i když by teoreticky bylo možné do tohoto registru uložit přímo první prvek pole (vyžadujeme ovšem chování podle normy IEEEE 754):
array_sum:
pxor xmm0, xmm0
addss xmm0, DWORD PTR [rdi]
addss xmm0, DWORD PTR [rdi+4]
addss xmm0, DWORD PTR [rdi+8]
addss xmm0, DWORD PTR [rdi+12]
ret
Zajímavé je, že i když automatickou vektorizaci povolíme, výsledný kód se nijak nezmění a bude stále používat „skalární“ instrukce:
array_sum:
pxor xmm0, xmm0
addss xmm0, DWORD PTR [rdi]
addss xmm0, DWORD PTR [rdi+4]
addss xmm0, DWORD PTR [rdi+8]
addss xmm0, DWORD PTR [rdi+12]
ret
4. Automatická vektorizace po specifikaci přepínače -ffast-math
Jak ovšem dosáhneme toho, aby se náš zdrojový kód přeložil do „vektorizované“ podoby? V tomto případě budeme muset použít přepínač -ffast-math, který překladači umožní vložit do výsledného kódu instrukce, jejichž chování nebude přesně odpovídat normě IEEE754 (pořadí operací, práce s denormalizovanými hodnotami atd.), takže může nastat situace, kdy ten stejný výpočet (na úrovni původního zdrojového kódu) povede k nepatrně odlišným výsledkům. To nám ovšem mnohdy nevadí.
Výsledná vektorizovaná podoba vypadá takto:
array_sum:
movups xmm0, XMMWORD PTR [rdi]
movaps xmm1, xmm0
movhlps xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
addps xmm0, xmm1
ret
5. Vysvětlení postupu výpočtu
Co vlastně výše uvedená sekvence instrukcí znamená? Zkusme si ji popsat postupně, a to na vstupním vektoru obsahujícím hodnoty [1000, 100, 10, 1] (schválně nepoužívám hodnoty 1, 2, 3, 4, protože jejich mezisoučty jsou nejasné). Nejdříve se čtveřice hodnot [1000, 100, 10, 1] přesune do registru XMM0 a posléze i do registru XMM1. Obsah obou registrů tedy bude vypadat následovně:
XMM0: 00001000 00000100 00000010 00000001 XMM1: 00001000 00000100 00000010 00000001
Posléze se instrukcí MOVHLPS přenesou nejvyšší dva prvky z registru XMM0 do spodních dvou prvků registru XMM1, takže obsah registrů bude vypadat takto:
XMM0: 00001000 00000100 00000010 00000001 XMM1: 00001000 00000100 00001000 00000100
Po vektorovém součtu instrukcí ADDPS bude v registru XMM1 uložen vektor s prvky [2000, 200, 1010, 101], tj. spodní dva prvky již budou obsahovat dva mezisoučty (označme 101 jako první mezisoučet a 1010 jako mezisoučet druhý):
XMM0: 00001000 00000100 00000010 00000001 XMM1: 00002000 00000200 00001010 00000101
Instrukcemi MOVAPS XMM0, XMM1 a SHUFPS xmm0, xmm1, 01010101 se do registru XMM0 uloží prvky s indexem jedna z registrů XMM0 a XMM1, tj. konkrétně druhý mezisoučet 1010. V registru XMM1 je v horních prvcích uloženo smetí a v nejnižším prvku pak stále máme první mezisoučet 101:
XMM0: 00001010 00001010 00001010 00001010 XMM1: 00002000 00000200 00001010 00000101
Poslední součet instrukcí ADDPS nám již dá v nejnižším prvku vektoru XMM0 kýžený výsledek:
XMM0: 00003010 00001210 00002020 00001111
Vše si můžeme otestovat tímto kódem (přeložitelným v Linuxu):
[bits 32]
%include "linux_macros.asm"
;-----------------------------------------------------------------------------
section .data
hex_message:
times 8 db '?'
db ' '
hex_message_length equ $ - hex_message
align 16
sse_val_1 dd 0x1, 0x10, 0x100, 0x1000
;-----------------------------------------------------------------------------
section .bss
sse_tmp resb 16
;-----------------------------------------------------------------------------
section .text
global _start ; tento symbol ma byt dostupny i linkeru
_start:
mov ebx, sse_val_1
movaps xmm0, [ebx] ; nacteni prvni hodnoty do registru XMM0
movaps xmm1, xmm0 ; xmm1 := xmm0
print_sse_reg_as_hex xmm0 ; tisk hodnoty registru XMM0
print_sse_reg_as_hex xmm1 ; tisk hodnoty registru XMM1
movhlps xmm1, xmm0 ; přesun nejvyšších dvou prvků vektoru
print_sse_reg_as_hex xmm0 ; tisk hodnoty registru XMM0
print_sse_reg_as_hex xmm1 ; tisk hodnoty registru XMM1
addps xmm1, xmm0 ; vektorový součet
print_sse_reg_as_hex xmm0 ; tisk hodnoty registru XMM0
print_sse_reg_as_hex xmm1 ; tisk hodnoty registru XMM1
movaps xmm0, xmm1
shufps xmm0, xmm1, 85 ; proložení prvků vektorů
print_sse_reg_as_hex xmm0 ; tisk hodnoty registru XMM0
print_sse_reg_as_hex xmm1 ; tisk hodnoty registru XMM1
addps xmm0, xmm1 ; vektorový součet
print_sse_reg_as_hex xmm0 ; tisk hodnoty registru XMM0
exit ; ukonceni procesu
%include "hex2string.asm"
6. Druhý příklad: součet všech prvků v krátkém jednorozměrném poli s osmi prvky
Zkusme si nyní celé pole, jehož prvky se sčítají, zvětšit ze čtyř prvků na prvků osm. Tím donutíme překladač, aby součet provedl minimálně dvakrát, protože do XMM registru je možné uložit maximálně čtyři prvky typu single/float. Příklad je upraven pouze nepatrně – pozměnili jsme konstantu SIZE:
float array_sum_8(float *a) {
#define SIZE 8
int i;
float result = 0.0;
for (i=0; i<SIZE; i++) {
result += a[i];
}
return result;
}
7. Výsledné kódy vygenerované překladačem: „skalární“ i „vektorová“ varianta
Opět se podívejme na výsledné sekvence instrukcí vygenerovaných překladačem. V případě, že jsou optimalizace zakázány, je výsledkem (podle očekávání) stále stejná neoptimalizovaná počítaná programová smyčka:
array_sum:
push rbp
mov rbp, rsp
mov QWORD PTR [rbp-24], rdi
pxor xmm0, xmm0
movss DWORD PTR [rbp-8], xmm0
mov DWORD PTR [rbp-4], 0
jmp .L2
.L3:
mov eax, DWORD PTR [rbp-4]
cdqe
lea rdx, [0+rax*4]
mov rax, QWORD PTR [rbp-24]
add rax, rdx
movss xmm0, DWORD PTR [rax]
movss xmm1, DWORD PTR [rbp-8]
addss xmm0, xmm1
movss DWORD PTR [rbp-8], xmm0
add DWORD PTR [rbp-4], 1
.L2:
cmp DWORD PTR [rbp-4], 7
jle .L3
movss xmm0, DWORD PTR [rbp-8]
pop rbp
ret
Při zapnutí optimalizací, ovšem bez použití autovektorizace, nyní získáme počítanou smyčku; sice krátkou a efektivní, ale jedná se o smyčku:
array_sum:
lea rax, [rdi+32]
pxor xmm0, xmm0
.L2:
addss xmm0, DWORD PTR [rdi]
add rdi, 8
addss xmm0, DWORD PTR [rdi-4]
cmp rdi, rax
jne .L2
ret
Ovšem při povolení automatické vektorizace a současně i výpočtů neodpovídajících přesně normě IEEE 754 („rychlá matematika“) je výsledek mnohem lepší. Dostaneme vlastně podobný kód, jako u pole se čtyřmi prvky, nyní ovšem hned na začátku proběhne součet prvku číslo 0 s prvkem číslo 4, dále prvku číslo 1 s prvkem číslo 5 atd. – tedy první mezisoučet, jehož výsledkem jsou čtyři částečné součty. A následnou sekvenci instrukcí již známe:
array_sum:
movups xmm1, XMMWORD PTR [rdi+16]
movups xmm0, XMMWORD PTR [rdi]
addps xmm0, xmm1
movaps xmm1, xmm0
movhlps xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
addps xmm0, xmm1
ret
8. Třetí příklad: součet pole se 100 prvky
Naposledy se ještě podívejme na problematiku součtu prvků polí, tentokrát pro rozsáhlé pole se sto prvky. Takový počet by již překladač měl donutit k vygenerování programové smyčky (i když si dovedu představit situaci, kdy bude vhodnější celý kód rozbalit). Naivní a nijak neoptimalizovaná funkce pro součet prvků pole o délce 100 bude vypadat takto:
float array_sum_100(float *a) {
#define SIZE 100
int i;
float result = 0.0;
for (i=0; i<SIZE; i++) {
result += a[i];
}
return result;
}
9. Výsledné kódy vygenerované překladačem: „skalární“ i „vektorová“ varianta
Skalární, ovšem optimalizovaná varianta výpočtu, bude vypadat poněkud zajímavě, protože z ní můžeme vyčíst, že se v jedné iteraci smyčky provádí přičtení hodnoty dvou prvků. Smyčka je tedy částečně rozbalena. Taktéž je užitečné, že se hodnota registru RDI (adresa prvku) zvyšuje před provedením součtu, což sice poněkud komplikuje adresování, ale po součtu již bude možné novou hodnotu porovnat s mezní hodnotou. Tím se vynechá (překryje) čekání na výsledek přičtení osmičky k RDI:
array_sum100:
lea rax, [rdi+400]
pxor xmm0, xmm0
.L2:
addss xmm0, DWORD PTR [rdi]
add rdi, 8
addss xmm0, DWORD PTR [rdi-4]
cmp rdi, rax
jne .L2
ret
Vektorová varianta je vlastně variací na známé téma, tentokrát je ovšem akumulátor (výsledek celého součtu) uložen v registru XMM0 a vektorové součty se provádí s registry XMM1 a XMM2:
array_sum100:
lea rax, [rdi+400]
pxor xmm0, xmm0
.L2:
movups xmm2, XMMWORD PTR [rdi]
add rdi, 16
addps xmm0, xmm2
cmp rax, rdi
jne .L2
movaps xmm1, xmm0
movhlps xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
addps xmm0, xmm1
ret
10. Čtvrtý příklad: vyhledání největšího prvku v poli
Další ukázkou operace typu reduce je vyhledání největšího prvku v poli. To vlastně znamená opakované testování hodnoty prvků vůči předchozí maximální hodnotě, takže bude zajímavé zjistit, zda i tento kód je možné vektorizovat, resp. jestli to dokáže překladač GCC. Začneme opět tím nejjednodušším možným případem: polem se čtyřmi prvky:
#include <float.h>
float find_max_4(float *a) {
#define SIZE 4
int i;
float max = FLT_MIN_EXP;
for (i=0; i<SIZE; i++) {
if (a[i] > max) {
max = a[i];
}
}
return max;
}
11. Nevektorizovaná a vektorizovaná podoba výsledného kódu
Podobně jako u předchozích demonstračních příkladů se i nyní podíváme na výsledný kód vygenerovaný překladačem GCC. Nejprve je ukázána skalární varianta založená na programové smyčce, ve které se volá „skalární“ instrukce MAXSS, s níž jsme se již setkali.
find_max_4:
movss xmm0, DWORD PTR .LC0[rip]
lea rax, [rdi+16]
.L2:
movss xmm1, DWORD PTR [rdi]
add rdi, 4
maxss xmm1, xmm0
movaps xmm0, xmm1
cmp rdi, rax
jne .L2
ret
.LC0:
.long -1023803392
Optimalizovaná skalární varianta explicitně zavolá instrukci MAXSS postupně pro všechny prvky pole. Jedná se tedy o rozbalenou variantu programové smyčky, což je pochopitelně rychlejší řešení:
find_max_4:
movss xmm0, DWORD PTR [rdi]
movss xmm1, DWORD PTR [rdi+8]
maxss xmm0, DWORD PTR [rdi+4]
maxss xmm1, DWORD PTR .LC0[rip]
maxss xmm0, DWORD PTR [rdi+12]
maxss xmm0, xmm1
ret
.LC0:
.long -1023803392
A konečně se podívejme na vektorizovanou variantu. Povšimněte si, jak se tato varianta podobá výpočtu součtu prvků čtyřprvkového pole – ovšem pochopitelně se nyní namísto instrukce ADDPS volá instrukce MAXPS:
find_max_4:
movups xmm0, XMMWORD PTR [rdi]
movaps xmm1, xmm0
movhlps xmm1, xmm0
maxps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
maxps xmm0, xmm1
maxss xmm0, DWORD PTR .LC0[rip]
ret
.LC0:
.long -1023803392
Pro zajímavost si oba vektorizované výpočty porovnejme po jednotlivých instrukcích, aby ještě více vynikla jejich podobnost:
součet všech prvků pole
movups xmm0, XMMWORD PTR [rdi] movups xmm0, XMMWORD PTR [rdi]
movaps xmm1, xmm0 movaps xmm1, xmm0
movhlps xmm1, xmm0 movhlps xmm1, xmm0
addps xmm1, xmm0 maxps xmm1, xmm0
movaps xmm0, xmm1 movaps xmm0, xmm1
shufps xmm0, xmm1, 85 shufps xmm0, xmm1, 85
addps xmm0, xmm1 maxps xmm0, xmm1
ret maxss xmm0, DWORD PTR .LC0[rip]
ret
12. Vliv velikosti pole na výsledný vygenerovaný kód
Samozřejmě si otestujeme, jaký vliv bude mít velikost pole (vektoru) na vygenerovaný strojový kód nebo assembler. Pro pole o počtu prvků osm dostaneme následující výsledky.
Skalární varianta, ve které se už používá programová smyčka a v každé iteraci se zpracují dva prvky pole:
find_max_8:
movss xmm0, DWORD PTR .LC0[rip]
lea rax, [rdi+32]
.L2:
maxss xmm0, DWORD PTR [rdi]
add rdi, 8
maxss xmm0, DWORD PTR [rdi-4]
cmp rdi, rax
jne .L2
ret
.LC0:
.long -1023803392
Vektorová varianta, ve které překladač smyčku rozbalil, ovšem stále potřebuje čtyři instrukce MAXPS:
find_max_8:
movss xmm0, DWORD PTR .LC1[rip]
movups xmm1, XMMWORD PTR [rdi]
shufps xmm0, xmm0, 0
maxps xmm1, xmm0
movups xmm0, XMMWORD PTR [rdi+16]
maxps xmm0, xmm1
movaps xmm1, xmm0
movhlps xmm1, xmm0
maxps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
maxps xmm0, xmm1
ret
.LC1:
.long -1023803392
Následuje výsledek překladu s automatickou vektorizací pro pole, které obsahuje 100 prvků. To již pochopitelně vede k vygenerování programové smyčky namísto jejího rozbalení. Nicméně stále se v jedné iteraci nalezne maximální hodnota čtyř nových prvků a porovná se s průběžně „akumulovanou“ maximální hodnotou:
find_max_100:
movss xmm0, DWORD PTR .LC1[rip]
lea rax, [rdi+400]
shufps xmm0, xmm0, 0
.L2:
movups xmm2, XMMWORD PTR [rdi]
add rdi, 16
maxps xmm0, xmm2
cmp rax, rdi
jne .L2
movaps xmm1, xmm0
movhlps xmm1, xmm0
maxps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
maxps xmm0, xmm1
ret
.LC1:
.long -1023803392
13. Jeden z nejdůležitějších algoritmů současnosti
S rozvojem neuronových sítí a umělé inteligence (ale nejenom zde) se začal masivně využívat známý algoritmus pro výpočet skalárního součinu (dot product, scalar product). Tento algoritmus se používá v oblasti velkých jazykových modelů (LLM) pro zjišťování podobnosti dlouhých vektorů s numerickými hodnotami (vector similarity). Kromě klasického skalárního součinu se v této oblasti používá i tzv. cosinus similarity, což je varianta skalárního součinu, v níž nezáleží na délce vektorů, ale pouze na jejich vzájemné orientaci (výpočet je tedy doplněn o normalizaci vektorů). A toto porovnávání vektorů se v LLM provádí neustále a většinou je optimalizováno a výpočty běží na GPU.
To však není zdaleka vše. Pokud se zaměříme na oblast klasických neuronových sítí (NN – neural networks), zjistíme, že se tyto sítě skládají z takzvaných perceptronů, což je vlastně značně zjednodušený model neuronu. A na vstup perceptronu se přivádí nějaké množství numerických vstupů. Každý z těchto vstupů je váhován, tj. vynásoben určitou konstantou a výsledky tohoto váhování jsou nakonec sečteny. Když se ovšem nad touto operací zamyslíme, zjistíme, že se vlastně nejedná o nic jiného, než o aplikaci výpočtu skalárního součinu. První z vektorů, který do tohoto součinu vstupuje jako operand, jsou vstupy do neuronu, druhým vektorem je pak vektor vah, které si neuron zapamatoval. A samotný trénink neuronové sítě vlastně není nic jiného, než rekonfigurace těchto vah – vektorů:
Obrázek 1: Idealizovaný model neuronu.
Obrázek 2: Idealizovaný model neuronu s biasem.
14. Pátý příklad: skalární součin dvou krátkých polí se čtyřmi prvky
V dnešním pátém ukázkovém příkladu je realizován algoritmus pro výpočet skalárního součinu dvou čtyřprvkových vektorů. Funkce dot_product4 je přitom napsána zcela naivním způsobem, bez jakékoli snahy o její optimalizaci na úrovni céčkového kódu:
float dot_product_4(float *a, float *b) {
#define SIZE 4
int i;
float result = 0.0;
for (i=0; i<SIZE; i++) {
result += a[i] * b[i];
}
return result;
}
15. Nevektorizovaná a vektorizovaná podoba výsledného kódu
Opět se, podobně jako u předchozích příkladů, podíváme na výsledek překladu výše uvedeného výpočtu do strojového kódu. Nejprve je uvedena nevektorizovaná varianta. Ta je v praxi dosti nezajímavá – moc pomalá, protože provádí čtyři součiny a čtyři součty:
dot_product_4:
xor eax, eax
pxor xmm1, xmm1
.L2:
movss xmm0, DWORD PTR [rdi+rax]
mulss xmm0, DWORD PTR [rsi+rax]
add rax, 4
addss xmm1, xmm0
cmp rax, 16
jne .L2
movaps xmm0, xmm1
ret
Vektorizovaná varianta vypadá následovně. Jedná se o variantu kompatibilní s normou IEEE 754:
dot_product_4:
movups xmm0, XMMWORD PTR [rsi]
movups xmm1, XMMWORD PTR [rdi]
pxor xmm2, xmm2
mulps xmm1, xmm0
movaps xmm0, xmm1
addss xmm0, xmm2
movaps xmm2, xmm1
shufps xmm2, xmm1, 85
addss xmm0, xmm2
movaps xmm2, xmm1
unpckhps xmm2, xmm1
shufps xmm1, xmm1, 255
addss xmm0, xmm2
addss xmm0, xmm1
ret
16. Další možné optimalizace výpočtů skalárního součinu
Nejkratší kód (bez smyčky) dostaneme v případě povolení autovektorizace a současně i povolením režimu „rychlé matematiky“ (fast-math), takže výsledky obecně nemusí přesně odpovídat normě IEEE 754. Tato varianta je založena na součinu prvků vektorů, součtu prvků vektorů a následně ještě jednoho součtu, který vrátí výslednou skalární hodnotu:
dot_product_4:
movups xmm1, XMMWORD PTR [rsi]
movups xmm0, XMMWORD PTR [rdi]
mulps xmm0, xmm1
movaps xmm1, xmm0
movhlps xmm1, xmm0
addps xmm1, xmm0
movaps xmm0, xmm1
shufps xmm0, xmm1, 85
addps xmm0, xmm1
ret
17. Vysvětlení postupu výpočtu
Sekvenci instrukcí z předchozí kapitoly si opět popíšeme podrobněji. Přitom provedeme výpočet skalárního součinu pro dvojici vektorů:
a = (1, 2, 3, 4) b = (5, 6, 7, 1)
Ručně by se výpočet provedl takto:
a · b = a1·b1 + a2·b2 + a3·b3 + a4·b4 = 1·5 + 2·6 + 3·7 + 4·1 = 5 + 12 + 21 + 4 = 42
Nejdříve jsou instrukcemi MOVUPS načteny oba vektory do pracovních registrů XMM1 a XMM0 (u každého registru je nejdříve vypsána jeho hexadecimální hodnota a poté odpovídající FP hodnoty):
XMM0: 40A00000 40C00000 40E00000 3F800000 5.0 6.0 7.0 1.0 XMM1: 3F800000 40000000 40400000 40800000 1.0 2.0 3.0 4.0
Instrukcí MULPS je proveden součin odpovídajících si prvků vektorů a výsledek je uložen do registru XMM0:
XMM0: 40A00000 41400000 41A80000 40800000 5.0 12.0 21.0 4.0 XMM1: 3F800000 40000000 40400000 40800000 1.0 2.0 3.0 4.0
XMM1 nyní neobsahuje žádné použitelné hodnoty, takže ho využijeme pro další výpočty. Nejdříve kopie XMM0 do XMM1:
XMM0: 40A00000 41400000 41A80000 40800000 5.0 12.0 21.0 4.0 XMM1: 40A00000 41400000 41A80000 40800000 5.0 12.0 21.0 4.0
A následně kopie pouze vyšších dvou prvků do prvků nižších v cílovém registru (provedeno instrukcí MOVHLPS. Nyní máme první dva mezivýsledky součinů v registru XMM1 jak na horních dvou, tak i na dolních dvou pozicích:
XMM0: 40A00000 41400000 41A80000 40800000 5.0 12.0 21.0 4.0 XMM1: 40A00000 41400000 40A00000 41400000 5.0 12.0 5.0 12.0
Oba vektory sečteme, takže získáme první dva mezisoučty:
XMM0: 40A00000 41400000 41A80000 40800000 5.0 12.0 21.0 4.0 XMM1: 41200000 41C00000 41D00000 41800000 10.0 24.0 26.0 16.0
Následuje známý trik spočívající v kombinovaném použití MOVAPS+SHUFPS s konstantou 85=0b01010101:
XMM0: 41D00000 41D00000 41D00000 41D00000 26.0 26.0 26.0 26.0 XMM1: 41200000 41C00000 41D00000 41800000 10.0 24.0 26.0 16.0
Nakonec nám postačuje součet (dokonce by mohlo jít o skalární součet, protože nás zajímá jen nejnižší prvek):
XMM0: 42100000 42480000 42500000 42280000 36.0 50.0 52.0 42.0 XMM1: 41200000 41C00000 41D00000 41800000 10.0 24.0 26.0 16.0
V registru XMM0 je v nejnižším prvku skutečně výsledek 42.0!
18. Repositář s demonstračními příklady
Demonstrační příklady naprogramované v jazyku, které jsou určené pro překlad s využitím assembleru gcc, byly uloženy do Git repositáře, který je dostupný na adrese https://github.com/tisnik/8bit-fame. Kromě zdrojových kódů příkladů jsou do repositáře přidány i výsledky překladu do assembleru v syntaxi kompatibilní s Intelem. Jednotlivé demonstrační příklady si můžete v případě potřeby stáhnout i jednotlivě bez nutnosti klonovat celý (dnes již poměrně rozsáhlý) repositář:
19. Seznam všech předchozích částí tohoto seriálu a článků o SIMD instrukcích
Podporou SIMD instrukcí na úrovni intrinsic jsme se už na Rootu zabývali, stejně jako samotnými SIMD instrukcemi na úrovni assembleru. Pro úplnost jsou v této příloze uvedeny odkazy na příslušné články:
- Užitečné rozšíření GCC: podpora SIMD (vektorových) instrukcí
https://www.root.cz/clanky/uzitecne-rozsireni-gcc-podpora-simd-vektorovych-instrukci/ - Užitečné rozšíření GCC – podpora SIMD (vektorových) instrukcí: nedostatky technologie
https://www.root.cz/clanky/uzitecne-rozsireni-gcc-podpora-simd-vektorovych-instrukci-nedostatky-technologie/ - Podpora SIMD (vektorových) instrukcí na RISCových procesorech
https://www.root.cz/clanky/podpora-simd-vektorovych-instrukci-na-riscovych-procesorech/ - Podpora SIMD operací v GCC s využitím intrinsic pro nízkoúrovňové optimalizace
https://www.root.cz/clanky/podpora-simd-operaci-v-gcc-s-vyuzitim-intrinsic-pro-nizkourovnove-optimalizace/ - Podpora SIMD operací v GCC s využitím intrinsic: technologie SSE
https://www.root.cz/clanky/podpora-simd-operaci-v-gcc-s-vyuzitim-intrinsic-technologie-sse/ - Rozšíření instrukční sady „Advanced Vector Extensions“ na platformě x86–64
https://www.root.cz/clanky/rozsireni-instrukcni-sady-advanced-vector-extensions-na-platforme-x86–64/ - Rozšíření instrukční sady F16C, FMA a AVX-512 na platformě x86–64
https://www.root.cz/clanky/rozsireni-instrukcni-sady-f16c-fma-a-avx-512-na-platforme-x86–64/ - Rozšíření instrukční sady AVX-512 na platformě x86–64 (dokončení)
https://www.root.cz/clanky/rozsireni-instrukcni-sady-avx-512-na-platforme-x86–64-dokonceni/ - SIMD instrukce na platformě 80×86: instrukční sada MMX
https://www.root.cz/clanky/simd-instrukce-na-platforme-80×86-instrukcni-sada-mmx/ - SIMD instrukce na 80×86: dokončení popisu MMX, instrukce 3DNow!
https://www.root.cz/clanky/simd-instrukce-na-80–86-dokonceni-popisu-mmx-instrukce-3dnow/ - SIMD instrukce v rozšíření SSE
https://www.root.cz/clanky/simd-instrukce-v-rozsireni-sse/ - SIMD instrukce v rozšíření SSE (2. část)
https://www.root.cz/clanky/simd-instrukce-v-rozsireni-sse-2-cast/ - Pokročilejší SSE operace: přeskupení, promíchání a rozbalování prvků vektorů
https://www.root.cz/clanky/pokrocilejsi-sse-operace-preskupeni-promichani-a-rozbalovani-prvku-vektoru/ - Od instrukční sady SSE k sadě SSE2
https://www.root.cz/clanky/od-instrukcni-sady-sse-k-sade-sse2/ - Instrukční sady SIMD a automatické vektorizace prováděné překladačem GCC
https://www.root.cz/clanky/instrukcni-sady-simd-a-automaticke-vektorizace-provadene-prekladacem-gcc/
20. Odkazy na Internetu
- Auto-vectorization in GCC
https://gcc.gnu.org/projects/tree-ssa/vectorization.html - GCC documentation: Extensions to the C Language Family
https://gcc.gnu.org/onlinedocs/gcc/C-Extensions.html#C-Extensions - GCC documentation: Using Vector Instructions through Built-in Functions
https://gcc.gnu.org/onlinedocs/gcc/Vector-Extensions.html - SSE (Streaming SIMD Extentions)
http://www.songho.ca/misc/sse/sse.html - Timothy A. Chagnon: SSE and SSE2
http://www.cs.drexel.edu/~tc365/mpi-wht/sse.pdf - Intel corporation: Extending the Worldr's Most Popular Processor Architecture
http://download.intel.com/technology/architecture/new-instructions-paper.pdf - SIMD architectures:
http://arstechnica.com/old/content/2000/03/simd.ars/ - Tour of the Black Holes of Computing!: Floating Point
http://www.cs.hmc.edu/~geoff/classes/hmc.cs105…/slides/class02_floats.ppt - 3Dnow! Technology Manual
AMD Inc., 2000 - Intel MMXTM Technology Overview
Intel corporation, 1996 - MultiMedia eXtensions
http://softpixel.com/~cwright/programming/simd/mmx.phpi - AMD K5 („K5“ / „5k86“)
http://www.pcguide.com/ref/cpu/fam/g5K5-c.html - Sixth Generation Processors
http://www.pcguide.com/ref/cpu/fam/g6.htm - Great Microprocessors of the Past and Present
http://www.cpushack.com/CPU/cpu1.html - Very long instruction word (Wikipedia)
http://en.wikipedia.org/wiki/Very_long_instruction_word - CPU design (Wikipedia)
http://en.wikipedia.org/wiki/CPU_design - Bulldozer (microarchitecture)
https://en.wikipedia.org/wiki/Bulldozer_(microarchitecture) - SIMD Instructions Considered Harmful
https://www.sigarch.org/simd-instructions-considered-harmful/ - GCC Compiler Intrinsics
https://iq.opengenus.org/gcc-compiler-intrinsics/ - Scalable_Vector_Extension_(SVE)
https://en.wikipedia.org/wiki/AArch64#Scalable_Vector_Extension_(SVE) - Improve the Multimedia User Experience
https://www.arm.com/technologies/neon - NEON Technology (stránky ARM)
https://developer.arm.com/technologies/neon - SIMD Assembly Tutorial: ARM NEON – Xiph.org
https://people.xiph.org/~tterribe/daala/neon_tutorial.pdf - Ne10
http://projectne10.github.io/Ne10/ - NEON and Floating-Point architecture
http://infocenter.arm.com/help/index.jsp?topic=/com.arm.doc.den0024a/BABIGHEB.html - An Introduction to ARM NEON
http://peterdn.com/post/an-introduction-to-ARM-NEON.aspx - ARM NEON Intrinsics Reference
http://infocenter.arm.com/help/topic/com.arm.doc.ihi0073a/IHI0073A_arm_neon_intrinsics_ref.pdf - Arm Neon Intrinsics vs hand assembly
https://stackoverflow.com/questions/9828567/arm-neon-intrinsics-vs-hand-assembly - ARM NEON Optimization. An Example
http://hilbert-space.de/?p=22 - AArch64 NEON instruction format
https://developer.arm.com/docs/den0024/latest/7-aarch64-floating-point-and-neon/73-aarch64-neon-instruction-format - ARM SIMD instructions
https://developer.arm.com/documentation/dht0002/a/Introducing-NEON/What-is-SIMD-/ARM-SIMD-instructions - Learn the architecture – Migrate Neon to SVE Version 1.0
https://developer.arm.com/documentation/102131/0100/?lang=en - 1.2.2. Comparison between NEON technology and other SIMD solutions
https://developer.arm.com/documentation/den0018/a/Introduction/Comparison-between-ARM-NEON-technology-and-other-implementations/Comparison-between-NEON-technology-and-other-SIMD-solutions?lang=en - NEON Programmer’s Guide
https://documentation-service.arm.com/static/63299276e68c6809a6b41308 - Brain Floating Point – nový formát uložení čísel pro strojové učení a chytrá čidla
https://www.root.cz/clanky/brain-floating-point-ndash-novy-format-ulozeni-cisel-pro-strojove-uceni-a-chytra-cidla/ - Other Built-in Functions Provided by GCC
https://gcc.gnu.org/onlinedocs/gcc/Other-Builtins.html - GCC: 6.60 Built-in Functions Specific to Particular Target Machines
https://gcc.gnu.org/onlinedocs/gcc/Target-Builtins.html#Target-Builtins - Advanced Vector Extensions
https://en.wikipedia.org/wiki/Advanced_Vector_Extensions - Top 10 Craziest Assembly Language Instructions
https://www.youtube.com/watch?v=Wz_xJPN7lAY - Intel x86: let's take a look at one of the most complex instruction set!
https://www.youtube.com/watch?v=KBLy23B38-c - x64 Assembly Tutorial 58: Intro to AVX
https://www.youtube.com/watch?v=yAvuHd8cBJY - AVX512 (1 of 3): Introduction and Overview
https://www.youtube.com/watch?v=D-mM6X5×nTY - AVX512 (2 of 3): Programming AVX512 in 3 Different Ways
https://www.youtube.com/watch?v=I3efQKLgsjM - AVX512 (3 of 3): Deep Dive into AVX512 Mechanisms
https://www.youtube.com/watch?v=543a1b-cPmU - AVX-512
https://en.wikipedia.org/wiki/AVX-512 - AVX-512
https://iq.opengenus.org/avx512/ - SIMD Algorithms Youtube course
https://denisyaroshevskiy.github.io/presentations/ - Compiler explorer
https://godbolt.org/ - Restricting pointers
https://gcc.gnu.org/onlinedocs/gcc/Restricted-Pointers.html - Does the restrict keyword provide significant benefits in gcc/g++
https://stackoverflow.com/questions/1965487/does-the-restrict-keyword-provide-significant-benefits-in-gcc-g - Demystifying The Restrict Keyword
https://cellperformance.beyond3d.com/articles/2006/05/demystifying-the-restrict-keyword.html - Basics of Vectorization for Fortran Applications
https://inria.hal.science/hal-01688488/document - What does the restrict keyword mean in C++?
https://stackoverflow.com/questions/776283/what-does-the-restrict-keyword-mean-in-c - restrict keyword (Wikipedia)
https://en.wikipedia.org/wiki/Restrict - Reduction operator
https://en.wikipedia.org/wiki/Reduction_operator - The Power of the Dot Product in Artificial Intelligence
https://medium.com/data-science/the-power-of-the-dot-product-in-artificial-intelligence-c002331e1829 - Can any one explain why dot product is used in neural network and what is the intitutive thought of dot product
https://stats.stackexchange.com/questions/291680/can-any-one-explain-why-dot-product-is-used-in-neural-network-and-what-is-the-in
ČLÁNKY DO MAILU