Wektoryzacja Obliczanie iloczynu punktowego za pomocą SSE4

Question

Próbuję poprawić ten kod za pomocą produktu kropki SSE4, ale mam problem ze znalezieniem rozwiązania. Ta funkcja pobiera parametry qi i tj, które zawierają tablice float z 80 komórkami, a następnie obliczają iloczyn kropki. Wartość zwracana to wektor z produktami z czterema kropkami. Więc próbuję to zrobić równolegle obliczając czteropunktowe produkty o dwudziestu wartościach.

Masz pojęcie, jak ulepszyć ten kod?

inline __m128 ScalarProd20Vec(__m128* qi, __m128* tj) 
{ 
    __m128 res=_mm_add_ps(_mm_mul_ps(tj[0],qi[0]),_mm_mul_ps(tj[1],qi[1])); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[2],qi[2]),_mm_mul_ps(tj[3],qi[3]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[4],qi[4]),_mm_mul_ps(tj[5],qi[5]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[6],qi[6]),_mm_mul_ps(tj[7],qi[7]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[8],qi[8]),_mm_mul_ps(tj[9],qi[9]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[10],qi[10]),_mm_mul_ps(tj[11],qi[11]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[12],qi[12]),_mm_mul_ps(tj[13],qi[13]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[14],qi[14]),_mm_mul_ps(tj[15],qi[15]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[16],qi[16]),_mm_mul_ps(tj[17],qi[17]))); 
    res=_mm_add_ps(res,_mm_add_ps(_mm_mul_ps(tj[18],qi[18]),_mm_mul_ps(tj[19],qi[19]))); 
    return res; 
} 

Answer 1

Z setek przykładów SSE, które widziałem na SO, Twój kod jest jednym z niewielu, który jest już w całkiem dobrej formie od samego początku. Nie potrzebujesz instrukcji dot. Produktu SSE4. (Można to zrobić lepiej!)

Jednak jest jedna rzecz, można spróbować: (. Mówię spróbować, bo nie timed go jeszcze)

Obecnie masz łańcuch danych o zależnościach res. Dodawanie wektorowe to 3-4 cykle na większości komputerów.Więc kod zajmie minimum 30 cykli uruchomić ponieważ masz:

(10 additions on critical path) * (3 cycles addps latency) = 30 cycles 

Co można zrobić, to do węzła-Split zmienną res następująco:

__m128 res0 = _mm_add_ps(_mm_mul_ps(tj[ 0],qi[ 0]),_mm_mul_ps(tj[ 1],qi[ 1])); 
__m128 res1 = _mm_add_ps(_mm_mul_ps(tj[ 2],qi[ 2]),_mm_mul_ps(tj[ 3],qi[ 3])); 

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[ 4],qi[ 4]),_mm_mul_ps(tj[ 5],qi[ 5]))); 
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[ 6],qi[ 6]),_mm_mul_ps(tj[ 7],qi[ 7]))); 

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[ 8],qi[ 8]),_mm_mul_ps(tj[ 9],qi[ 9]))); 
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[10],qi[10]),_mm_mul_ps(tj[11],qi[11]))); 

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[12],qi[12]),_mm_mul_ps(tj[13],qi[13]))); 
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[14],qi[14]),_mm_mul_ps(tj[15],qi[15]))); 

res0 = _mm_add_ps(res0,_mm_add_ps(_mm_mul_ps(tj[16],qi[16]),_mm_mul_ps(tj[17],qi[17]))); 
res1 = _mm_add_ps(res1,_mm_add_ps(_mm_mul_ps(tj[18],qi[18]),_mm_mul_ps(tj[19],qi[19]))); 

return _mm_add_ps(res0,res1); 

To prawie odcina swoją krytyczny ścieżka na pół. Zauważ, że ze względu na nie-łączność zmiennoprzecinkową, ta optymalizacja jest niedozwolona dla kompilatorów.

---

Oto alternatywna wersja z 4-stronnymi instrukcjami dzielenia węzłów i AMD FMA4. Jeśli nie możesz użyć dodanych fuzji-mnożenia, możesz je podzielić. To może być lepsze niż pierwsza wersja powyżej.

__m128 res0 = _mm_mul_ps(tj[ 0],qi[ 0]); 
__m128 res1 = _mm_mul_ps(tj[ 1],qi[ 1]); 
__m128 res2 = _mm_mul_ps(tj[ 2],qi[ 2]); 
__m128 res3 = _mm_mul_ps(tj[ 3],qi[ 3]); 

res0 = _mm_macc_ps(tj[ 4],qi[ 4],res0); 
res1 = _mm_macc_ps(tj[ 5],qi[ 5],res1); 
res2 = _mm_macc_ps(tj[ 6],qi[ 6],res2); 
res3 = _mm_macc_ps(tj[ 7],qi[ 7],res3); 

res0 = _mm_macc_ps(tj[ 8],qi[ 8],res0); 
res1 = _mm_macc_ps(tj[ 9],qi[ 9],res1); 
res2 = _mm_macc_ps(tj[10],qi[10],res2); 
res3 = _mm_macc_ps(tj[11],qi[11],res3); 

res0 = _mm_macc_ps(tj[12],qi[12],res0); 
res1 = _mm_macc_ps(tj[13],qi[13],res1); 
res2 = _mm_macc_ps(tj[14],qi[14],res2); 
res3 = _mm_macc_ps(tj[15],qi[15],res3); 

res0 = _mm_macc_ps(tj[16],qi[16],res0); 
res1 = _mm_macc_ps(tj[17],qi[17],res1); 
res2 = _mm_macc_ps(tj[18],qi[18],res2); 
res3 = _mm_macc_ps(tj[19],qi[19],res3); 

res0 = _mm_add_ps(res0,res1); 
res2 = _mm_add_ps(res2,res3); 

return _mm_add_ps(res0,res2); 

Answer 2

Po pierwsze, najważniejszą optymalizacją, jaką możesz zrobić, jest upewnienie się, że Twój kompilator ma włączone wszystkie ustawienia optymalizacji.

---

Kompilatory są dość inteligentne, więc jeśli napisać go jako pętla, jest prawdopodobne, aby rozwinąć go:

__128 res = _mm_setzero(); 
for (int i = 0; i < 10; i++) { 
    res = _mm_add_ps(res, _mm_add_ps(_mm_mul_ps(tj[2*i], qi[2*i]), _mm_mul_ps(tj[2*i+1], qi[2*i+1]))); 
} 
return res; 

(z gcc trzeba przejść -funroll-loops, a następnie będzie ona go rozwinąć zrobić 5 powtórzeń na raz)

można również zdefiniować makro i rozwinąć go ręcznie, jeśli wersja pętli jest wolniejszy, np.

__128 res = _mm_setzero(); 

#define STEP(i) res = _mm_add_ps(res, _mm_add_ps(_mm_mul_ps(tj[2*i], qi[2*i]), _mm_mul_ps(tj[2*i+1], qi[2*i+1]))) 

STEP(0); STEP(1); STEP(2); STEP(3); STEP(4); 
STEP(5); STEP(6); STEP(7); STEP(8); STEP(9); 

#undef STEP 

return res; 

---

Można nawet uruchomić pętlę od 0 do 20 (lub zrobić to samo z makro wersji), czyli:

__128 res = _mm_setzero(); 
for (int i = 0; i < 20; i++) { 
    res = _mm_add_ps(res, _mm_mul_ps(tj[i], qi[i])); 
} 
return res; 

(z gcc i -funroll-loops to rozwinął się zrobić 10 powtórzeń na raz , tj. taki sam jak pętla "dwa w czasie" powyżej).

Answer 3

Twoje dane nie są uporządkowane w pamięci w odpowiednim formacie dla wyspecjalizowanych instrukcji produktu SSE4 dot (dpps). Instrukcje te oczekują wymiary pojedynczego wektora być przylegający ,, tak:

| dim0 | dim1 | dim2 | ... | dim19 | 

natomiast dane wydaje się mieć wektory przeplecione ze sobą:

| v0-dim0 | v1-dim0 | v2-dim0 | v3-dim0 | v0-dim1 | ... 

Aktualna ogólne podejście wydaje się właściwe - możesz poprawić rzeczy, zmieniając kolejność instrukcji tak, aby wyniki mnożenia nie były używane natychmiast po ich wygenerowaniu, ale tak naprawdę kompilator powinien być w stanie sam to sobie wyobrazić.