Voxelzeug mit Visual C++ 2017 Update 2 x64. Umständlich zu erklären, aber die Wirkung von Mikrooptimierungen grenzt hier an Magie.
Ich habe 8×8×8
float-Voxel eines Signed Distance Fields in einem Block (Standardmaße eines Leaf Nodes bei OpenVDB). Ich möchte die 2×2×2 Voxel bei Position XYZ haben, um sie kubisch interpolieren zu können. (Ich nehme niemals die Voxel ganz am Rand, dafür habe ich bereits vorgesorgt.)
OpenVDB schreibt die ersten acht Werte entlang Z in ein Array. Dann rutscht es in Y weiter und holt sich wieder acht Werte entlang der Z-Achse. Das klingt unintuitiv (serialisiert man nicht normalerweise in X-Y-Z statt Z-Y-X?), erlaubt aber, Schleifen in
for(x) for(y) for(z) zu schreiben, was wohl einfacher lesbar sein soll.
Okay. Hier mein erster Entwurf (schon doppelt so schnell wie OpenVDB):
Code: Alles auswählen
auto firstIndex = 64 * x + 8 * y + z;
auto toDistances = leaf.buffer().mData;
float result[2 * 2 * 2];
result[0] = toDistances[firstIndex ]; // x y z
result[1] = toDistances[firstIndex + 1]; // x y z+1
result[2] = toDistances[firstIndex + 8 ]; // x y+1 z
result[3] = toDistances[firstIndex + 8 + 1]; // x y+1 z+1
result[4] = toDistances[firstIndex + 64 ]; // x+1 y z
result[5] = toDistances[firstIndex + 64 + 1]; // x+1 y z+1
result[6] = toDistances[firstIndex + 64 + 8 ]; // x+1 y+1 z
result[7] = toDistances[firstIndex + 64 + 8 + 1]; // x+1 y+1 z+1
Ihr erkennt den Trick: Wenn man zum Index 8 addiert, entspricht das einem Schritt entlang Y, weil der Block in Z acht Voxel groß ist. Addiert man 64, geht man einen Schritt entlang X.
Laufzeit im Testfall: 14.76 s.
Schaut man in den Assembler-Code, wird einem mulmig:
Code: Alles auswählen
; result[2] = toDistances[firstIndex + 8];
lea eax,[rdx+8] ; firstIndex + 8
movsxd rcx,eax ; Index kopieren; Grund ist CPU-Komplexität
movss xmm1,dword ptr [r8+rcx*4] ; mit sizeof(float) multiplizieren; zur Adresse des Arrays addieren = Adresse des floats; von dort laden
Drei Befehle, um die Adresse eines Array-Elements auszurechnen. Pfff. Das können wir selber besser: Die Offsets sind konstant, und da können wir
sizeof(float) von Hand aufmultiplizieren! Statt acht Elemente weiterzugehen, gehen wir also 32 Bytes weiter. Statt 64 Elemente, 256 B.
Code: Alles auswählen
template <typename T> T * byteOffset(T * pointer, size_t offsetInBytes) { // eine der meistgenutzten Funktionen bei mir
return (T *)(((char *)pointer) + offsetInBytes);
}
auto firstOffset = 256 * x + 32 * y + 4 * z;
auto toDistances = l0.buffer().mData;
float result[2 * 2 * 2];
result[0] = *byteOffset(toDistances, firstOffset ); // x y z
result[1] = *byteOffset(toDistances, firstOffset + 4); // x y z+1
result[2] = *byteOffset(toDistances, firstOffset + 32 ); // x y+1 z
result[3] = *byteOffset(toDistances, firstOffset + 32 + 4); // x y+1 z+1
result[4] = *byteOffset(toDistances, firstOffset + 256 ); // x+1 y z
result[5] = *byteOffset(toDistances, firstOffset + 256 + 4); // x+1 y z+1
result[6] = *byteOffset(toDistances, firstOffset + 256 + 32 ); // x+1 y+1 z
result[7] = *byteOffset(toDistances, firstOffset + 256 + 32 + 4); // x+1 y+1 z+1
Die Befehle haben sich so geringfügig geändert, dass man den Unterschied kaum erkennt:
Code: Alles auswählen
lea eax,[r9+20h]
movsxd rcx,eax
movss xmm1,dword ptr [rcx+rdx] ; keine Multiplikation mehr!
Trotzdem …
Laufzeit: 13.89 s (6 % schneller).
… da haben wir dem Decoder wohl eine µOp abgenommen, in die er
„lade r8+rcx*4“ normalerweise zerlegt hätte. Nett.
Aber da geht noch was: Warum eigentlich zwei Mal addieren? Erst Offset zum ersten Index, dann das Ergebnis zur Array-Adresse? Also:
Code: Alles auswählen
auto toDistances = byteOffset(l0.buffer().mData, 256 * x + 32 * y + 4 * z);
float result[2 * 2 * 2];
result[0] = *toDistances; // x y z
result[1] = *byteOffset(toDistances, + 4); // x y z+1
result[2] = *byteOffset(toDistances, + 32 ); // x y+1 z
result[3] = *byteOffset(toDistances, + 32 + 4); // x y+1 z+1
result[4] = *byteOffset(toDistances, 256 ); // x+1 y z
result[5] = *byteOffset(toDistances, 256 + 4); // x+1 y z+1
result[6] = *byteOffset(toDistances, 256 + 32 ); // x+1 y+1 z
result[7] = *byteOffset(toDistances, 256 + 32 + 4); // x+1 y+1 z+1
Hier war ich skeptisch, denn eigentlich vertiefen wir so die Abhängigkeitskette. Vorher hätte jede Adressberechnung unabhängig von der ersten Addition ausgeführt werden können, aber so … ich poste einfach mal das Disassembly der kompletten Funktion vorher:
Code: Alles auswählen
sub rsp,38h
movaps xmmword ptr [rsp+20h],xmm6
movaps xmmword ptr [rsp+10h],xmm7
movaps xmm7,xmm0
shufps xmm7,xmm0,55h
movd eax,xmm1
pshufd xmm2,xmm1,55h
movaps xmmword ptr [rsp],xmm8
movaps xmm8,xmm0
movhlps xmm6,xmm8
movd edx,xmm2
punpckhdq xmm1,xmm1
lea r8d,[rdx+rax*8]
mov rdx,qword ptr [rcx]
movd eax,xmm1
lea r9d,[rax+r8*8]
shl r9d,2
movsxd rax,r9d
movss xmm0,dword ptr [rax+rdx]
lea eax,[r9+4]
movsxd rcx,eax
lea eax,[r9+20h]
movss xmm4,dword ptr [rcx+rdx]
subss xmm4,xmm0
movsxd rcx,eax
lea eax,[r9+100h]
movss xmm1,dword ptr [rcx+rdx]
movsxd rcx,eax
lea eax,[r9+104h]
mulss xmm4,xmm6
addss xmm4,xmm0
movss xmm0,dword ptr [rcx+rdx]
movsxd rcx,eax
lea eax,[r9+120h]
movss xmm5,dword ptr [rcx+rdx]
subss xmm5,xmm0
movsxd rcx,eax
lea eax,[r9+24h]
movss xmm2,dword ptr [rcx+rdx]
movsxd rcx,eax
lea eax,[r9+124h]
mulss xmm5,xmm6
movss xmm3,dword ptr [rcx+rdx]
addss xmm5,xmm0
movsxd rcx,eax
subss xmm3,xmm1
movss xmm0,dword ptr [rcx+rdx]
subss xmm0,xmm2
mulss xmm3,xmm6
mulss xmm0,xmm6
addss xmm3,xmm1
movaps xmm6,xmmword ptr [rsp+20h]
addss xmm0,xmm2
subss xmm3,xmm4
subss xmm0,xmm5
mulss xmm3,xmm7
mulss xmm0,xmm7
addss xmm3,xmm4
movaps xmm7,xmmword ptr [rsp+10h]
addss xmm0,xmm5
subss xmm0,xmm3
mulss xmm0,xmm8
movaps xmm8,xmmword ptr [rsp]
addss xmm0,xmm3
add rsp,38h
ret
… und nacher:
Code: Alles auswählen
sub rsp,38h
movaps xmmword ptr [rsp+20h],xmm6
movaps xmmword ptr [rsp+10h],xmm7
movaps xmm7,xmm0
shufps xmm7,xmm0,55h
movd eax,xmm1
movaps xmmword ptr [rsp],xmm8
movaps xmm8,xmm0
movhlps xmm6,xmm8
pshufd xmm2,xmm1,55h
movd edx,xmm2
punpckhdq xmm1,xmm1
lea r8d,[rdx+rax*8]
movd eax,xmm1
lea eax,[rax+r8*8]
shl eax,2
cdqe
add rax,qword ptr [rcx]
movss xmm0,dword ptr [rax+100h]
movss xmm2,dword ptr [rax+120h]
movss xmm5,dword ptr [rax+104h]
movss xmm3,dword ptr [rax+4]
subss xmm5,xmm0
subss xmm3,dword ptr [rax]
movss xmm4,dword ptr [rax+24h]
subss xmm4,dword ptr [rax+20h]
mulss xmm5,xmm6
mulss xmm3,xmm6
addss xmm5,xmm0
mulss xmm4,xmm6
addss xmm3,dword ptr [rax]
movss xmm0,dword ptr [rax+124h]
addss xmm4,dword ptr [rax+20h]
subss xmm0,xmm2
subss xmm4,xmm3
mulss xmm0,xmm6
movaps xmm6,xmmword ptr [rsp+20h]
addss xmm0,xmm2
mulss xmm4,xmm7
addss xmm4,xmm3
subss xmm0,xmm5
mulss xmm0,xmm7
movaps xmm7,xmmword ptr [rsp+10h]
addss xmm0,xmm5
subss xmm0,xmm4
mulss xmm0,xmm8
movaps xmm8,xmmword ptr [rsp]
addss xmm0,xmm4
add rsp,38h
ret
Da muss man kein Assembler-Champion sein, um den Unterschied zu sehen. Laufzeit:
13.46 s (3 % schneller).
Beachtet, dass das ein realer Anwendungsfall ist, in dem drumherum sehr viel gerechnet wird; kein synthetischer Benchmark, der einfach nur die Werte lädt. Ich habe die ganze Nacht Cache-Lokalität von Bäumen optimiert, und es hat
nichts gebracht. Jetzt ist mein Programm zehn Prozent schneller – durch so einen Kleinkram.