Jaki jest koszt obliczeniowy

26

Jednym z głównych problemów, z którymi musimy się zmierzyć w symulacjach molekularnych, jest obliczanie sił zależnych od odległości. Jeśli możemy ograniczyć funkcje siły i odległości do uzyskania równych mocy odległości separacji , możemy po prostu obliczyć kwadrat odległości r 2 = r r i nie musimy się martwić o r . Jeśli jednak istnieją moce nieparzyste, musimy poradzić sobie z r = rr2=rrr .r=r2

Moje pytanie brzmi: jak drogie jest przetwarzanie komputerowe zaimplementowane w bibliotekach popularnych języków (C / C ++, Fortran, Python) itp.? Czy naprawdę trzeba poprawić wydajność ręcznie, dostosowując kod do konkretnych architektur?x

eeismail
źródło

Odpowiedzi:

39

Jako rozszerzenie odpowiedź moyner za , on-Chip sqrtjest zwykle rsqrt, czyli odwrotność pierwiastka kwadratowego, który oblicza 1 / . Więc jeśli w swoim kodzie będziesz używał tylko1/r(jeśli robisz dynamikę molekularną, jesteś), możesz obliczyćbezpośrednio i zapisać sobie podział. Powodem, dla któregooblicza się zamiast tegojest to, że jego iteracja Newtona nie ma podziałów, tylko uzupełnienia i mnożenia.a1/a1/rr = rsqrt(r2)rsqrtsqrt

Na marginesie, podziały są również obliczane iteracyjnie i są prawie tak samo wolne jak rsqrtw sprzęcie. Jeśli szukasz wydajności, lepiej jest usunąć zbędne podziały.

Niektóre bardziej nowoczesne architektury, takie jak architektury POWER IBM, nie zapewniają rsqrtper se, ale oszacowanie z dokładnością do kilku bitów, np . FRSQRTE . Gdy użytkownik dzwoni rsqrt, generuje to oszacowanie, a następnie jedną lub dwie (tyle, ile potrzeba) iteracje algorytmu Newtona lub Goldschmidta przy użyciu regularnych mnożeń i dodatków. Zaletą tego podejścia jest to, że etapy iteracji mogą być potokowane i przeplatane z innymi instrukcjami bez blokowania FPU (bardzo ładny przegląd tej koncepcji, aczkolwiek na starszych architekturach, patrz praca doktorska Rolfa Strebela ).

W przypadku potencjałów interakcji sqrtoperacji można całkowicie uniknąć, stosując interpolant wielomianowy funkcji potencjalnej, ale moja własna praca (zaimplementowana w mdcore) w tym obszarze pokazuje, że przynajmniej w architekturach typu x86 sqrtinstrukcja jest wystarczająco szybka.

Aktualizacja

Ponieważ wydaje się, że ta odpowiedź przyciąga sporo uwagi, chciałbym również odnieść się do drugiej części twojego pytania, tj. Czy naprawdę warto spróbować ulepszyć / wyeliminować podstawowe operacje, takie jak sqrt?

W kontekście symulacji dynamiki molekularnej lub jakiejkolwiek symulacji opartej na cząsteczkach z interakcjami ograniczonymi przez granicę odcięcia można wiele zyskać dzięki lepszym algorytmom znajdowania sąsiadów. Jeśli używasz list komórek lub czegoś podobnego, aby znaleźć sąsiadów lub utworzyć listę Verlet , będziesz obliczał dużą liczbę fałszywych odległości parami. W naiwnym przypadku tylko 16% skontrolowanych par cząstek faktycznie znajduje się w odległości odcięcia od siebie. Chociaż dla takich par nie jest obliczana interakcja, dostęp do danych cząstek i obliczenie fałszywej odległości parami wiąże się z dużymi kosztami.

Moja własna praca w tym obszarze ( tutaj , tutaj i tutaj ), a także innych (np. Tutaj ), pokazuje, jak można uniknąć tych fałszywych obliczeń. Te algorytmy wyszukiwania sąsiadów przewyższają nawet listy Verletów, jak opisano tutaj .

Chciałbym podkreślić, że chociaż mogą istnieć pewne ulepszenia, które można uzyskać dzięki lepszej znajomości / wykorzystaniu podstawowej architektury sprzętowej, istnieją również potencjalnie większe korzyści z ponownego przemyślenia algorytmów wyższego poziomu.

Pedro
źródło
6
SSE rsqrtpsi AVX vrsqrtpssą również szacunkami, mają poprawne pierwsze 11 do 12 bitów i powinieneś udoskonalić iteracją Newtona lub dwoma, jeśli chcesz większej dokładności. Są to instrukcje 5/1 i 7/1 (opóźnienie / odwrotna przepustowość) na Sandy Bridge (patrz dokumentacja Intela lub tabele instrukcji Agner Fog, które są porównywalne do mnożenia. W przeciwieństwie do tego, pełna dokładność (v)sqrtps(lub podwójna precyzja (v)sqrtpd) zajmuje 10-43 / 10–43 (szczegóły w tabelach instrukcji)
Jed Brown
@JedBrown: Dziękujemy za zwrócenie na to uwagi! Zapomniałem, że SSE i jego rozszerzenia również to zapewniają.
Pedro
16

Pierwiastek kwadratowy jest zaimplementowany w sprzęcie na większości procesorów, co oznacza, że ​​istnieją specjalne instrukcje montażu, a wydajność powinna być porównywalna w większości języków, ponieważ bardzo trudno jest zepsuć implementację. Prawdopodobnie nigdy nie będziesz w stanie pokonać instrukcji FSQRT, ponieważ została ona zaprojektowana przez jakiegoś inteligentnego projektanta sprzętu.

Sposób implementacji w sprzęcie może się różnić, ale prawdopodobnie jest to pewnego rodzaju iteracja punktu stałego, na przykład metoda Newtona-Raphsona, która wykonuje określoną liczbę iteracji do momentu obliczenia wymaganej liczby cyfr. Metody iteracyjne w sprzęcie są ogólnie znacznie wolniejsze niż inne operacje, ponieważ kilka cykli musi zostać zakończonych, zanim wynik będzie gotowy.

Istnieją również pewne instrukcje Streaming SIMD , które mogą być stosowane w tych rejestrach xmm dla szybkich obliczeń wektorowych znaleźć tutaj . Rejestry te są dość małe, ale jeśli masz znaną liczbę współrzędnych (powiedzmy trójwymiarowy kartezjański układ współrzędnych), mogą one być nieco szybsze.

Jeśli Twój język jest wystarczająco niski, zawsze możesz rzutować czcionką na niższą dokładność lub użyć mniejszej dokładności dla swoich współrzędnych. Pojedyncza precyzja jest często więcej niż wystarczająca i z tego, co pamiętam, będzie szybciej przy obliczaniu pierwiastków kwadratowych, ponieważ iteracje można zakończyć wcześniej.

Analiza porównawcza różnych języków powinna być łatwa: po prostu napisz do pliku długą serię liczb losowych, załaduj go w różnych językach, a następnie zmierz pierwiastki kwadratowe.

Moyner
źródło
0

Mogą istnieć ulepszenia wydajności, ale najpierw należy profilować, aby wiedzieć, że obliczenie odwrotności sqrt jest szyjką butelki (a nie, powiedzmy, ładowaniem pozycji i oszczędzaniem sił).

Projekt GROMACS MD wyrósł z pomysłu wykorzystania szczegółów formatu zmiennoprzecinkowego IEEE w celu zaszczepienia schematu iteracji Newtona-Raphsona w celu obliczenia akceptowalnego przybliżenia do odwrotności pierwiastka kwadratowego (patrz Dodatek B.3 http: / /www.gromacs.org/Documentation/Manual ), ale nie ma używanych procesorów HPC, w których GROMACS nadal korzysta z tego pomysłu.

Mabraham
źródło