Bezpiecznie i szybko FFT

Question

Zainspirowany przez przekonujący wykład Herb Sutter'a Not your father's C++, Postanowiłem spojrzeć jeszcze raz na najnowszą wersję C++ przy użyciu Visual Studio 2010 Microsoftu. Szczególnie interesowało mnie twierdzenie Herba, że ​​C++ jest "bezpieczny i szybki", ponieważ napisz dużo kodu krytycznego pod względem wydajności.

Jako benchmark postanowiłem spróbować napisać ten sam prosty algorytm FFT w różnych językach.

wpadłem na następujący kod C++ 11, który wykorzystuje wbudowany w complex rodzaju i vector kolekcji:

#include <complex> 
#include <vector> 

using namespace std; 

// Must provide type or MSVC++ barfs with "ambiguous call to overloaded function" 
double pi = 4 * atan(1.0); 

void fft(int sign, vector<complex<double>> &zs) { 
    unsigned int j=0; 
    // Warning about signed vs unsigned comparison 
    for(unsigned int i=0; i<zs.size()-1; ++i) { 
     if (i < j) { 
      auto t = zs.at(i); 
      zs.at(i) = zs.at(j); 
      zs.at(j) = t; 
     } 
     int m=zs.size()/2; 
     j^=m; 
     while ((j & m) == 0) { m/=2; j^=m; } 
    } 
    for(unsigned int j=1; j<zs.size(); j*=2) 
     for(unsigned int m=0; m<j; ++m) { 
      auto t = pi * sign * m/j; 
      auto w = complex<double>(cos(t), sin(t)); 
      for(unsigned int i = m; i<zs.size(); i+=2*j) { 
       complex<double> zi = zs.at(i), t = w * zs.at(i + j); 
       zs.at(i) = zi + t; 
       zs.at(i + j) = zi - t; 
      } 
     } 
} 

Należy pamiętać, że funkcja ta działa tylko dla n -elementowe wektorami gdzie n jest integralną moc dwóch. Każdy, kto szuka szybkiego kodu FFT, który działa dla każdego n powinien wyglądać na FFTW.

Jak rozumiem, tradycyjna składnia xs[i] z C dla indeksowania vector nie sprawdza granic, a zatem nie jest bezpieczna w pamięci i może być źródłem błędów pamięci, takich jak niedeterministyczna korupcja i naruszenia praw dostępu do pamięci . Więc zamiast tego użyłem xs.at(i).

Teraz chcę, aby ten kod był "bezpieczny i szybki", ale nie jestem ekspertem C++ 11, więc chciałbym poprosić o ulepszenia tego kodu, które uczyniłyby go bardziej idiotycznym lub wydajnym?

Answer 1

Myślę, że jesteś zbyt "bezpieczny" w używaniu funkcji at(). W większości przypadków użyty indeks jest trywialnie dający się sklasyfikować jako ograniczony przez rozmiar kontenera w pętli for.

np.

for(unsigned int i=0; i<zs.size()-1; ++i) { 
    ... 
    auto t = zs.at(i); 

Jedynymi, że zostawię na dzień (s) są (i + j) s. Nie jest od razu oczywiste, czy zawsze będą ograniczeni (chociaż gdybym był naprawdę niepewny, to prawdopodobnie sprawdziłbym ręcznie - ale nie znam się dobrze na FFT, żeby mieć opinię w tym przypadku).

Istnieją także środki Obliczenia powtarza się dla każdej iteracji:

int m=zs.size()/2; 
pi * sign 
2*j 

i The zs.at (i + j) są obliczane dwa razy.

Jest możliwe, że optymalizator może je złapać - ale jeśli traktujesz to jako krytyczne dla wydajności i masz testujących go timerów, wyciągnę je z pętli (lub, w przypadku z.at) (i + j), po prostu weź referencję na temat pierwszego użycia) i zobacz, czy to wpływa na timer.

Mówiąc o drugim odgadnięciu optymalizatora: Jestem pewien, że wywołania funkcji .size() zostaną zainicjowane jako przynajmniej bezpośrednie połączenie z wewnętrzną zmienną składową - ale zważywszy ile razy to nazwiesz Eksperymentował także z wprowadzaniem lokalnych zmiennych dla zs.size() i zs.size() - 1 z góry. Bardziej prawdopodobne jest, że będą również umieszczane w rejestrach.

Nie wiem, ile różnicy (jeśli w ogóle) ma to wszystko w twoim całkowitym czasie wykonywania - niektóre z nich mogą już zostać złapane przez optymalizatora, a różnice mogą być niewielkie w porównaniu do obliczeń - ale wart strzał.

Jeśli chodzi o bycie idiomatycznym, to moim jedynym komentarzem jest to, że size() zwraca std :: size_t (zwykle jest to typedef dla unsigned int - ale bardziej idiomatyczne jest używanie tego typu). Jeśli chcesz używać auto, ale unikaj ostrzeżenia, możesz spróbować dodać sufiks ul do 0 - nie jestem pewien, czy powiedziałbym, że jest to idiomatyczne. Przypuszczam, że już nie jesteś idiomatyczny, nie używając tutaj iteratorów, ale widzę, dlaczego nie możesz tego zrobić (łatwo).

Aktualizacja

dałem wszystkie moje sugestie spróbować i wszyscy mieli wymierny wzrost wydajności - z wyjątkiem I + J i 2 * j precalcs - faktycznie spowodował nieznaczne spowolnienie! Zakładam, że albo uniemożliwiły optymalizację kompilatora, albo uniemożliwiły mu korzystanie z rejestrów dla niektórych rzeczy.

Ogólnie mam> 10% perf. ulepszenie dzięki tym sugestiom. Podejrzewam, że można by było uzyskać więcej, gdyby drugi blok pętli został nieco refaktoryzowany, aby uniknąć skoków - i po wykonaniu tej czynności zestaw instrukcji SSE2 może dać znaczący impuls (spróbowałem go tak, jak jest i zauważyłem lekkie spowolnienie).

Myślę, że refaktoryzacja, wraz z użyciem czegoś takiego jak MKL dla wywołań cos i grzechów powinna dawać większe i mniej kruche, ulepszenia. Żadna z tych rzeczy nie byłaby zależna od języka (wiem, że pierwotnie porównywano ją do implementacji F #).

Aktualizacja 2

Zapomniałem wspomnieć, że wstępne obliczenia zs.size() zrobił zrobić różnicę.

Update 3

zapomniał również powiedzieć (do przypomnienia, przez @xeo w komentarzu do OP), że blok po j czeku i < można sprowadzić do std :: wymiany. Jest to bardziej idiomatyczne i co najmniej tak samo wydajne - w najgorszym przypadku powinno być zgodne z tym samym kodem, jaki został napisany. Rzeczywiście, kiedy to zrobiłem, nie widziałem żadnej zmiany w występie. W innych przypadkach może to prowadzić do zwiększenia wydajności, jeśli konstruktorzy ruchu są dostępni.