Dlaczego skalarny parametr sqrt (x) SSE jest wolniejszy niż rsqrt (x) * x?

106

Profilowałem część naszej podstawowej matematyki na Intel Core Duo i patrząc na różne podejścia do pierwiastka kwadratowego zauważyłem coś dziwnego: używając operacji skalarnych SSE, szybciej jest wziąć odwrotność pierwiastka kwadratowego i pomnożyć go aby uzyskać sqrt, niż użyć natywnego kodu operacji sqrt!

Testuję to z pętlą coś takiego:

inline float TestSqrtFunction( float in );

void TestFunc()
{
  #define ARRAYSIZE 4096
  #define NUMITERS 16386
  float flIn[ ARRAYSIZE ]; // filled with random numbers ( 0 .. 2^22 )
  float flOut [ ARRAYSIZE ]; // filled with 0 to force fetch into L1 cache

  cyclecounter.Start();
  for ( int i = 0 ; i < NUMITERS ; ++i )
    for ( int j = 0 ; j < ARRAYSIZE ; ++j )
    {
       flOut[j] = TestSqrtFunction( flIn[j] );
       // unrolling this loop makes no difference -- I tested it.
    }
  cyclecounter.Stop();
  printf( "%d loops over %d floats took %.3f milliseconds",
          NUMITERS, ARRAYSIZE, cyclecounter.Milliseconds() );
}

Próbowałem tego z kilkoma różnymi ciałami dla TestSqrtFunction i mam kilka czasów, które naprawdę drapią mnie po głowie. Zdecydowanie najgorsze było użycie natywnej funkcji sqrt () i pozwolenie „inteligentnemu” kompilatorowi na „optymalizację”. Przy 24ns / float, przy użyciu FPU x87 było to żałośnie złe:

inline float TestSqrtFunction( float in )
{  return sqrt(in); }

Następną rzeczą, jaką próbowałem, było użycie funkcji wewnętrznej, aby zmusić kompilator do użycia skalarnego kodu operacji sqrt SSE:

inline void SSESqrt( float * restrict pOut, float * restrict pIn )
{
   _mm_store_ss( pOut, _mm_sqrt_ss( _mm_load_ss( pIn ) ) );
   // compiles to movss, sqrtss, movss
}

To było lepsze, przy 11,9ns / float. Wypróbowałem również zwariowaną technikę przybliżenia Newtona-Raphsona Carmacka , która działała nawet lepiej niż sprzęt, przy 4,3 ns / float, chociaż z błędem 1 na 2 10 (co jest zbyt duże dla moich celów).

Doozy miał miejsce, gdy próbowałem opcją SSE dla odwrotności pierwiastka kwadratowego, a następnie użyłem mnożenia, aby uzyskać pierwiastek kwadratowy (x * 1 / √x = √x). Mimo to trwa dwie operacje zależne było najszybsze rozwiązanie zdecydowanie na 1.24ns / pływaka i dokładnością do 2 -14 :

inline void SSESqrt_Recip_Times_X( float * restrict pOut, float * restrict pIn )
{
   __m128 in = _mm_load_ss( pIn );
   _mm_store_ss( pOut, _mm_mul_ss( in, _mm_rsqrt_ss( in ) ) );
   // compiles to movss, movaps, rsqrtss, mulss, movss
}

Moje pytanie brzmi: co daje ? Dlaczego kod operacji pierwiastka kwadratowego wbudowany w sprzęt SSE jest wolniejszy niż jego synteza z dwóch innych operacji matematycznych?

Jestem pewien, że to tak naprawdę koszt samej operacji, bo zweryfikowałem:

  • Wszystkie dane mieszczą się w pamięci podręcznej, a dostęp jest sekwencyjny
  • funkcje są wbudowane
  • rozwijanie pętli nie robi różnicy
  • flagi kompilatora są ustawione na pełną optymalizację (a montaż jest dobry, sprawdziłem)

( edytuj : stephentyrone poprawnie wskazuje, że operacje na długich ciągach liczb powinny wykorzystywać wektoryzację operacji spakowanych w SIMD, na przykład rsqrtps- ale struktura danych tablicy jest tutaj tylko do celów testowych: to, co naprawdę próbuję zmierzyć, to wydajność skalarna do użycia w kodzie których nie można wektoryzować).

Crashworks
źródło
13
x / sqrt (x) = sqrt (x). Albo inaczej: x ^ 1 * x ^ (- 1/2) = x ^ (1 - 1/2) = x ^ (1/2) = sqrt (x)
Crashworks,
6
oczywiście inline float SSESqrt( float restrict fIn ) { float fOut; _mm_store_ss( &fOut, _mm_sqrt_ss( _mm_load_ss( &fIn ) ) ); return fOut; }. Jest to jednak zły pomysł, ponieważ może łatwo wywołać przeciągnięcie magazynu typu load-hit-store, jeśli procesor zapisuje liczby zmiennoprzecinkowe na stosie, a następnie odczytuje je natychmiast - żonglowanie z rejestru wektorowego do rejestru zmiennoprzecinkowego w szczególności dla wartości zwracanej to zła wiadomość. Poza tym podstawowe instrukcje maszynowe, które reprezentują elementy wewnętrzne SSE, i tak przyjmują operandy adresu.
Crashworks
4
To, ile LHS ma znaczenie, zależy od konkretnego genu i kroku danego x86: z mojego doświadczenia wynika, że ​​na czymkolwiek do i7 przenoszenie danych między zestawami rejestrów (np. FPU do SSE do eax) jest bardzo złe, podczas gdy podróż w obie strony między xmm0 a stosem a z powrotem nie jest ze względu na przekazywanie do sklepu przez firmę Intel. Możesz to zrobić samemu, aby to sprawdzić. Generalnie najłatwiejszym sposobem zobaczenia potencjalnego LHS jest przyjrzenie się emitowanemu złożeniu i zobaczenie, gdzie dane są żonglowane między zestawami rejestrów; Twój kompilator może zrobić mądrą rzecz lub może nie. Co do normalizacji wektorów, moje wyniki spisałem
Crashworks,
2
W przypadku PowerPC tak: IBM ma symulator procesora, który może przewidywać LHS i wiele innych bąbli potoków za pomocą analizy statycznej. Niektóre PPC mają również licznik sprzętowy dla LHS, który można sondować. Jest to trudniejsze dla x86; dobre narzędzia do profilowania są rzadsze (VTune jest obecnie nieco zepsuty), a uporządkowane potoki są mniej deterministyczne. Możesz spróbować zmierzyć to empirycznie, mierząc instrukcje na cykl, co można precyzyjnie wykonać za pomocą liczników wydajności sprzętu. Rejestry „wycofanych instrukcji” i „całkowitych cykli” można odczytać np. Za pomocą PAPI lub PerfSuite ( bit.ly/an6cMt ).
Crashworks
2
Możesz także po prostu napisać kilka permutacji funkcji i ustawić czas, aby sprawdzić, czy któraś z nich cierpi szczególnie na przeciągnięcia. Intel nie publikuje wielu szczegółów na temat sposobu działania ich rurociągów (że w ogóle LHS jest czymś w rodzaju brudnej tajemnicy), więc wiele z tego, czego nauczyłem się, to spojrzenie na scenariusz, który powoduje zatrzymanie się na innych łukach (np. ), a następnie konstruowanie kontrolowanego eksperymentu, aby sprawdzić, czy x86 również go ma.
Crashworks

Odpowiedzi:

216

sqrtssdaje poprawnie zaokrąglony wynik. rsqrtsspodaje przybliżenie odwrotności, z dokładnością do około 11 bitów.

sqrtssgeneruje znacznie dokładniejsze wyniki, gdy wymagana jest dokładność. rsqrtssistnieje w przypadkach, gdy wystarczy przybliżenie, ale wymagana jest prędkość. Jeśli przeczytasz dokumentację Intela, znajdziesz również sekwencję instrukcji (odwrotne przybliżenie pierwiastka kwadratowego, po którym następuje pojedynczy krok Newtona-Raphsona), która zapewnia prawie pełną precyzję (~ 23 bity dokładności, jeśli dobrze pamiętam) i nadal jest nieco szybciej niż sqrtss.

edycja: Jeśli szybkość ma kluczowe znaczenie i naprawdę wywołujesz to w pętli dla wielu wartości, powinieneś używać wektoryzowanych wersji tych instrukcji rsqrtpslub sqrtpsobie przetwarzają cztery zmiennoprzecinkowe na instrukcję.

Stephen Canon
źródło
3
Krok n / r zapewnia 22-bitową dokładność (podwaja ją); 23 bity oznaczałyby dokładnie pełną dokładność.
Jasper Bekkers
7
@Jasper Bekkers: Nie, nie byłoby. Po pierwsze, zmiennoprzecinkowa ma 24 bity precyzji. Po drugie, sqrtssjest poprawnie zaokrąglony , co wymaga ~ 50 bitów przed zaokrągleniem i nie można go osiągnąć za pomocą prostej iteracji N / R z pojedynczą precyzją.
Stephen Canon,
1
To jest zdecydowanie powód. Aby rozszerzyć ten wynik: projekt Embree firmy Intel ( software.intel.com/en-us/articles/… ) wykorzystuje w matematyce wektoryzację. Możesz pobrać źródło pod tym linkiem i zobaczyć, jak robią swoje wektory 3/4 D. Ich normalizacja wektorowa wykorzystuje rsqrt, po której następuje iteracja Newtona-Raphsona, która jest wtedy bardzo dokładna i nadal szybsza niż 1 / ssqrt!
Brandon Pelfrey
7
Małe zastrzeżenie: x rsqrt (x) daje NaN, jeśli x jest równe zero lub nieskończoność. 0 * rsqrt (0) = 0 * INF = NaN. INF rsqrt (INF) = INF * 0 = NaN. Z tego powodu CUDA na procesorach graficznych NVIDIA oblicza przybliżone pierwiastki kwadratowe o pojedynczej precyzji jako przepis (rsqrt (x)), przy czym sprzęt zapewnia zarówno szybkie przybliżenie odwrotności, jak i odwrotności pierwiastka kwadratowego. Oczywiście możliwe są również jawne kontrole obsługujące dwa specjalne przypadki (ale byłyby wolniejsze na GPU).
njuffa
@BrandonPelfrey W którym pliku znalazłeś krok Newtona Rhapsona?
fredoverflow
7

Dotyczy to również podziału. MULSS (a, RCPSS (b)) jest znacznie szybszy niż DIVSS (a, b). W rzeczywistości jest nadal szybszy, nawet jeśli zwiększysz jego precyzję za pomocą iteracji Newtona-Raphsona.

Intel i AMD zalecają tę technikę w swoich podręcznikach optymalizacji. W aplikacjach, które nie wymagają zgodności ze standardem IEEE-754, jedynym powodem używania div / sqrt jest czytelność kodu.

Sprzeczka
źródło
1
Broadwell i później mają lepszą wydajność dzielenia FP, więc kompilatory takie jak clang nie używają odwrotności + Newton dla wartości skalarnych na najnowszych procesorach, ponieważ zwykle nie jest to szybsze. W większości pętli divnie jest to jedyna operacja, więc łączna przepustowość UOP jest często wąskim gardłem, nawet jeśli występuje divpslub divss. Zobacz dzielenie zmiennoprzecinkowe a mnożenie zmiennoprzecinkowe , gdzie moja odpowiedź zawiera sekcję wyjaśniającą, dlaczego rcppsnie jest już wygrana w przepustowości. (Lub wygrana z opóźnieniem) i liczby dotyczące dzielenia przepustowości / opóźnienia.
Peter Cordes,
Jeśli twoje wymagania dotyczące dokładności są tak niskie, że możesz pominąć iterację Newtona, to tak a * rcpss(b)może być szybsze, ale nadal jest to więcej niż ups a/b!
Peter Cordes,
5

Zamiast udzielać odpowiedzi, która w rzeczywistości może być niepoprawna (nie zamierzam też sprawdzać ani dyskutować o pamięci podręcznej i innych rzeczach, powiedzmy, że są identyczne) spróbuję wskazać źródło, które może odpowiedzieć na twoje pytanie.
Różnica może polegać na sposobie obliczania sqrt i rsqrt. Więcej informacji można znaleźć tutaj http://www.intel.com/products/processor/manuals/ . Proponuję zacząć od przeczytania o funkcjach procesora, których używasz, jest trochę informacji, szczególnie o rsqrt (procesor używa wewnętrznej tabeli przeglądowej z ogromnym przybliżeniem, co znacznie ułatwia uzyskanie wyniku). Może się wydawać, że rsqrt jest o wiele szybszy niż sqrt, że 1 dodatkowa operacja mul (co nie jest zbyt kosztowna) może nie zmienić sytuacji tutaj.

Edycja: Kilka faktów, o których warto wspomnieć:
1. Kiedyś robiłem mikro optymalizacje dla mojej biblioteki graficznej i użyłem rsqrt do obliczania długości wektorów. (zamiast sqrt pomnożyłem sumę kwadratu przez rsqrt, co jest dokładnie tym, co zrobiłeś w swoich testach) i wypadło lepiej.
2. Obliczenie rsqrt przy użyciu prostej tabeli przeglądowej może być łatwiejsze, ponieważ dla rsqrt, gdy x dochodzi do nieskończoności, 1 / sqrt (x) idzie do 0, więc dla małych x wartości funkcji się nie zmieniają (dużo), podczas gdy dla sqrt - dąży do nieskończoności, więc to taki prosty przypadek;).

Ponadto wyjaśnienie: nie jestem pewien, gdzie znalazłem to w książkach, do których linkowałem, ale jestem prawie pewien, że czytałem, że rsqrt używa jakiejś tabeli odnośników i powinno być używane tylko wtedy, gdy wynik nie musi być dokładne, chociaż - ja też mogę się mylić, tak jak to było jakiś czas temu :).

Marcin Deptuła
źródło
4

Newton-Raphson zbiega się do zera f(x)przy użyciu przyrostów równych -f/f' gdzie f'jest pochodną.

Ponieważ x=sqrt(y)możesz spróbować rozwiązać f(x) = 0za xpomocą f(x) = x^2 - y;

Wtedy przyrost jest następujący: dx = -f/f' = 1/2 (x - y/x) = 1/2 (x^2 - y) / x który ma powolny podział.

Możesz wypróbować inne funkcje (np. f(x) = 1/y - 1/x^2), Ale będą one równie skomplikowane.

Spójrzmy 1/sqrt(y)teraz. Możesz spróbować f(x) = x^2 - 1/y, ale będzie to równie skomplikowane: dx = 2xy / (y*x^2 - 1)na przykład. Jednym z nieoczywistych alternatywnych opcji f(x)jest:f(x) = y - 1/x^2

Następnie: dx = -f/f' = (y - 1/x^2) / (2/x^3) = 1/2 * x * (1 - y * x^2)

Ach! To nie jest trywialne wyrażenie, ale masz w nim tylko mnożenia, bez dzielenia. => Szybciej!

I: pełny krok aktualizacji new_x = x + dxbrzmi:

x *= 3/2 - y/2 * x * x co też jest łatwe.

skal
źródło
2

Istnieje wiele innych odpowiedzi na to pytanie sprzed kilku lat. Oto, co jest słuszne w konsensusie:

  • Instrukcje rsqrt * obliczają przybliżenie do odwrotności pierwiastka kwadratowego, dobre do około 11-12 bitów.
  • Jest zaimplementowany z tablicą przeglądową (tj. ROM) indeksowaną przez mantysę. (W rzeczywistości jest to skompresowana tablica przeglądowa, podobna do starych tablic matematycznych, wykorzystująca korekty mniej znaczących bitów, aby zaoszczędzić na tranzystorach).
  • Powodem, dla którego jest dostępny, jest to, że jest to wstępne oszacowanie używane przez FPU dla "prawdziwego" algorytmu pierwiastka kwadratowego.
  • Istnieje również przybliżona wzajemna instrukcja, rcp. Obie te instrukcje są wskazówką, jak FPU implementuje pierwiastek kwadratowy i dzielenie.

Oto dlaczego konsensus się nie zgadzał:

  • Jednostki FPU ery SSE nie używają Newtona-Raphsona do obliczania pierwiastków kwadratowych. To świetna metoda w oprogramowaniu, ale byłoby błędem wdrażanie jej w ten sposób w sprzęcie.

Algorytm NR do obliczania odwrotności pierwiastka kwadratowego ma ten krok aktualizacji, jak zauważyli inni:

x' = 0.5 * x * (3 - n*x*x);

To dużo mnożenia zależnego od danych i jedno odejmowanie.

Poniżej znajduje się algorytm, którego faktycznie używają nowoczesne jednostki FPU.

Biorąc pod uwagę b[0] = n, załóżmy, że możemy znaleźć szereg liczb Y[i], które b[n] = b[0] * Y[0]^2 * Y[1]^2 * ... * Y[n]^2zbliżają się do 1. Następnie rozważ:

x[n] = b[0] * Y[0] * Y[1] * ... * Y[n]
y[n] = Y[0] * Y[1] * ... * Y[n]

Jasne x[n]podejście sqrt(n)i y[n]podejście 1/sqrt(n).

Możemy użyć kroku aktualizacji Newtona-Raphsona do odwrotności pierwiastka kwadratowego, aby uzyskać dobry Y[i]:

b[i] = b[i-1] * Y[i-1]^2
Y[i] = 0.5 * (3 - b[i])

Następnie:

x[0] = n Y[0]
x[i] = x[i-1] * Y[i]

i:

y[0] = Y[0]
y[i] = y[i-1] * Y[i]

Następna kluczowa obserwacja jest taka b[i] = x[i-1] * y[i-1]. Więc:

Y[i] = 0.5 * (3 - x[i-1] * y[i-1])
     = 1 + 0.5 * (1 - x[i-1] * y[i-1])

Następnie:

x[i] = x[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = x[i-1] + x[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))
y[i] = y[i-1] * (1 + 0.5 * (1 - x[i-1] * y[i-1]))
     = y[i-1] + y[i-1] * 0.5 * (1 - x[i-1] * y[i-1]))

Oznacza to, że mając początkowe x i y, możemy użyć następującego kroku aktualizacji:

r = 0.5 * (1 - x * y)
x' = x + x * r
y' = y + y * r

Lub, nawet bardziej wyszukane, możemy ustawić h = 0.5 * y. To jest inicjalizacja:

Y = approx_rsqrt(n)
x = Y * n
h = Y * 0.5

A to jest krok aktualizacji:

r = 0.5 - x * h
x' = x + x * r
h' = h + h * r

To jest algorytm Goldschmidta, który ma ogromną zaletę, jeśli implementujesz go w sprzęcie: „pętla wewnętrzna” to trzy wielokrotne dodawanie i nic więcej, a dwa z nich są niezależne i można je potokować.

W 1999 r. Jednostki FPU potrzebowały już potokowego obwodu dodawania / odejmowania i potokowego obwodu wielokrotnego, w przeciwnym razie SSE nie byłoby zbyt „strumieniowe”. W 1999 r. Potrzebny był tylko jeden z każdego obwodu, aby zaimplementować tę wewnętrzną pętlę w sposób w pełni potokowy bez marnowania dużej ilości sprzętu tylko na pierwiastek kwadratowy.

Dziś oczywiście połączyliśmy mnożenie i dodawanie ujawnione programiście. Ponownie, pętla wewnętrzna to trzy potokowe FMA, które są (znowu) ogólnie przydatne, nawet jeśli nie obliczasz pierwiastków kwadratowych.

Pseudonim
źródło
1
Powiązane: Jak działa sqrt () GCC po skompilowaniu? Która metoda rootowania jest używana? Newton-Raphson? ma pewne linki do projektów jednostek wykonawczych div / sqrt. Szybki wektoryzowany rsqrt i odwrotność z SSE / AVX w zależności od precyzji - jedna iteracja Newtona w oprogramowaniu, z lub bez FMA, do użytku z _mm256_rsqrt_psanalizą perf Haswell. Zwykle jest to dobry pomysł tylko wtedy, gdy nie masz innej pracy w pętli i utrudniłoby to przepustowość rozdzielacza. HW sqrt jest pojedynczym uopem, więc można go mieszać z innymi pracami.
Peter Cordes,
-2

Jest to szybsze, ponieważ te instrukcje ignorują tryby zaokrąglania i nie obsługują wyjątków zmiennoprzecinkowych ani zdernormalizowanych liczb. Z tych powodów znacznie łatwiej jest potokować, spekulować i wykonać inne instrukcje fp.

Witek
źródło
Oczywiście źle. FMA zależy od bieżącego trybu zaokrąglania, ale ma przepustowość dwóch na zegar w Haswell i nowszych. Dzięki dwóm w pełni rurociągowym jednostkom FMA Haswell może mieć w locie do 10 FMA jednocześnie. Prawidłowa odpowiedź to rsqrtjest znacznie niższa dokładność, co oznacza znacznie mniej pracy do zrobienia (albo wcale?) Po stołowych odnośnika dostać przypuszczenie wyjściowej.
Peter Cordes,