W jaki sposób BLAS osiąga tak ekstremalną wydajność?

Z ciekawości zdecydowałem się porównać moją własną funkcję mnożenia macierzy z implementacją BLAS-a ... Wynik był, powiem, najmniej zaskoczony:

Implementacja niestandardowa, 10 prób mnożenia macierzy 1000x1000:

Took: 15.76542 seconds.

Implementacja BLAS, 10 prób mnożenia macierzy 1000x1000:

Took: 1.32432 seconds.

To jest używanie liczb zmiennoprzecinkowych o pojedynczej precyzji.

Moje wdrożenie:

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
    if ( ADim2!=BDim1 )
        throw std::runtime_error("Error sizes off");

    memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
    int cc2,cc1,cr1;
    for ( cc2=0 ; cc2<BDim2 ; ++cc2 )
        for ( cc1=0 ; cc1<ADim2 ; ++cc1 )
            for ( cr1=0 ; cr1<ADim1 ; ++cr1 )
                C[cc2*ADim2+cr1] += A[cc1*ADim1+cr1]*B[cc2*BDim1+cc1];
}

Mam dwa pytania:

1. Biorąc pod uwagę, że mnożenie macierzy-macierzy mówi: nxm * mxn wymaga n * n * m mnożenia, czyli w przypadku powyżej 1000 ^ 3 lub 1e9 operacji. W jaki sposób na moim procesorze 2,6 GHz BLAS może wykonać 10 * 1e9 operacji w 1,32 sekundy? Nawet jeśli mnożenie było pojedynczą operacją i nic innego nie zostało zrobione, powinno to zająć ~ 4 sekundy.
2. Dlaczego moja implementacja jest o wiele wolniejsza?

Dobrym punktem wyjścia jest wspaniała książka „ The Science of Programming Matrix Computations” autorstwa Roberta A. van de Geijna i Enrique S. Quintana-Ortí. Zapewniają bezpłatną wersję do pobrania.

BLAS jest podzielony na trzy poziomy:

• Poziom 1 definiuje zestaw funkcji algebry liniowej, które działają tylko na wektorach. Funkcje te czerpią korzyści z wektoryzacji (np. Przy użyciu SSE).

• Funkcje poziomu 2 to operacje macierzowo-wektorowe, np. Jakiś iloczyn macierzowo-wektorowy. Funkcje te można by zaimplementować w kategoriach funkcji poziomu 1. Możesz jednak zwiększyć wydajność tej funkcji, jeśli możesz zapewnić dedykowaną implementację, która wykorzystuje architekturę wieloprocesorową ze współdzieloną pamięcią.

• Funkcje poziomu 3 to operacje takie jak iloczyn macierzy-macierzy. Ponownie możesz zaimplementować je w kategoriach funkcji Level2. Ale funkcje Level3 wykonują operacje O (N ^ 3) na danych O (N ^ 2). Jeśli więc Twoja platforma ma hierarchię pamięci podręcznej, możesz zwiększyć wydajność, jeśli zapewnisz dedykowaną implementację, która jest zoptymalizowana pod kątem pamięci podręcznej / przyjazna dla pamięci podręcznej . Jest to ładnie opisane w książce. Głównym wzmocnieniem funkcji Level3 jest optymalizacja pamięci podręcznej. To przyspieszenie znacznie przewyższa drugie wzmocnienie z równoległości i innych optymalizacji sprzętu.

Nawiasem mówiąc, większość (lub nawet wszystkie) wysokowydajnych implementacji BLAS NIE jest zaimplementowanych w Fortranie. ATLAS jest zaimplementowany w C. GotoBLAS / OpenBLAS jest zaimplementowany w C, a jego części krytyczne dla wydajności w Assemblerze. W Fortranie zaimplementowano tylko referencyjną implementację BLAS. Jednak wszystkie te implementacje BLAS zapewniają interfejs Fortran w taki sposób, że można go połączyć z LAPACK (LAPACK zyskuje całą swoją wydajność z BLAS).

Zoptymalizowane kompilatory odgrywają pod tym względem niewielką rolę (a dla GotoBLAS / OpenBLAS kompilator nie ma żadnego znaczenia).

IMHO żadna implementacja BLAS nie wykorzystuje algorytmów, takich jak algorytm Coppersmith-Winograd lub algorytm Strassena. Nie jestem do końca pewien, dlaczego tak jest, ale zgaduję:

• Może nie jest możliwe zapewnienie implementacji zoptymalizowanej pod kątem pamięci podręcznej tych algorytmów (tj. Straciłbyś więcej niż wygrał)
• Te algorytmy nie są stabilne numerycznie. Ponieważ BLAS jest obliczeniowym jądrem LAPACK, nie można tego zrobić.

Edycja / aktualizacja:

Nowym i przełomowym dokumentem na ten temat są dokumenty BLIS . Są wyjątkowo dobrze napisane. Na wykładzie "Podstawy oprogramowania do obliczeń o wysokiej wydajności" zaimplementowałem produkt macierzowo-macierzowy po ich artykule. Właściwie zaimplementowałem kilka wariantów produktu macierzowo-macierzowego. Najprostsze warianty są w całości napisane w czystym C i mają mniej niż 450 linii kodu. Wszystkie inne warianty tylko optymalizują pętle

    for (l=0; l<MR*NR; ++l) {
        AB[l] = 0;
    }
    for (l=0; l<kc; ++l) {
        for (j=0; j<NR; ++j) {
            for (i=0; i<MR; ++i) {
                AB[i+j*MR] += A[i]*B[j];
            }
        }
        A += MR;
        B += NR;
    }

Ogólna wydajność produktu macierz-macierz zależy tylko od tych pętli. Spędza się tu około 99,9% czasu. W pozostałych wariantach użyłem elementów wewnętrznych i kodu asemblera, aby poprawić wydajność. Możesz zobaczyć samouczek przedstawiający wszystkie warianty tutaj:

ulmBLAS: Samouczek dotyczący GEMM (produkt Matrix-Matrix)

Wraz z dokumentami BLIS dość łatwo można zrozumieć, w jaki sposób biblioteki takie jak Intel MKL mogą uzyskać taką wydajność. I dlaczego nie ma znaczenia, czy używasz pamięci głównej w postaci wierszy czy kolumn!

Końcowe testy porównawcze są tutaj (nazwaliśmy nasz projekt ulmBLAS):

Benchmarki dla ulmBLAS, BLIS, MKL, openBLAS i Eigen

Kolejna edycja / aktualizacja:

Napisałem również kilka poradników na temat tego, jak BLAS jest używany do rozwiązywania problemów z algebry liniowej numerycznej, takich jak rozwiązywanie układu równań liniowych:

Wysokowydajna faktoryzacja LU

(Ta faktoryzacja LU jest na przykład używana przez Matlab do rozwiązywania układu równań liniowych).

Mam nadzieję, że znajdę czas na rozszerzenie tego samouczka, aby opisać i zademonstrować, jak zrealizować wysoce skalowalną równoległą implementację faktoryzacji LU, jak w PLASMA .

OK, gotowe: kodowanie równoległej faktoryzacji LU zoptymalizowanej pod kątem pamięci podręcznej

PS: Zrobiłem też kilka eksperymentów nad poprawą wydajności uBLAS. W rzeczywistości jest całkiem proste, aby zwiększyć (tak, grać słowami :)) wydajność uBLAS:

Eksperymenty na uBLAS .

Tutaj podobny projekt z BLAZE :

Eksperymenty na BLAZE .

Przede wszystkim BLAS to tylko interfejs z około 50 funkcjami. Istnieje wiele konkurencyjnych implementacji interfejsu.

Najpierw wspomnę o rzeczach, które są w dużej mierze niezwiązane:

• Fortran vs C, nie ma znaczenia
• Zaawansowane algorytmy macierzowe, takie jak Strassen, nie używają ich, ponieważ nie pomagają w praktyce

Większość implementacji dzieli każdą operację na małe macierze lub operacje wektorowe w mniej lub bardziej oczywisty sposób. Na przykład duże mnożenie macierzy 1000x1000 może zostać podzielone na sekwencję mnożenia macierzy 50x50.

Te operacje o małych wymiarach o stałym rozmiarze (zwane jądrem) są zakodowane na stałe w kodzie asemblera specyficznym dla procesora przy użyciu kilku funkcji procesora docelowego:

• Instrukcje w stylu SIMD
• Równoległość poziomu instrukcji
• Świadomość pamięci podręcznej

Co więcej, jądra te mogą być wykonywane równolegle względem siebie przy użyciu wielu wątków (rdzeni procesora), w typowym wzorcu projektowym zmniejszania map.

Spójrz na ATLAS, który jest najczęściej używaną implementacją BLASa typu open source. Ma wiele różnych konkurujących ze sobą jąder, a podczas procesu budowania biblioteki ATLAS uruchamia konkurencję między nimi (niektóre są nawet sparametryzowane, więc to samo jądro może mieć różne ustawienia). Próbuje różnych konfiguracji, a następnie wybiera najlepszą dla konkretnego systemu docelowego.

(Wskazówka: dlatego, jeśli używasz ATLAS, lepiej jest zbudować i dostroić bibliotekę ręcznie dla konkretnej maszyny, niż używając wstępnie zbudowanej).

Po pierwsze, istnieją bardziej wydajne algorytmy mnożenia macierzy niż ten, którego używasz.

Po drugie, twój procesor może wykonać więcej niż jedną instrukcję naraz.

Twój procesor wykonuje 3-4 instrukcje na cykl, a jeśli używane są jednostki SIMD, każda instrukcja przetwarza 4 liczby zmiennoprzecinkowe lub 2 podwójne. (oczywiście ta liczba również nie jest dokładna, ponieważ procesor zwykle może przetwarzać tylko jedną instrukcję SIMD na cykl)

Po trzecie, twój kod jest daleki od optymalnego:

• Używasz surowych wskaźników, co oznacza, że ​​kompilator musi założyć, że mogą one alias. Istnieją słowa kluczowe lub flagi specyficzne dla kompilatora, które można określić, aby poinformować kompilator, że nie mają aliasów. Alternatywnie powinieneś użyć innych typów niż surowe wskaźniki, które rozwiązują problem.
• Niszczysz pamięć podręczną, wykonując naiwne przeglądanie każdego wiersza / kolumny macierzy wejściowych. Możesz użyć blokowania, aby wykonać jak najwięcej pracy na mniejszym bloku macierzy, który mieści się w pamięci podręcznej procesora, przed przejściem do następnego bloku.
• W przypadku zadań czysto numerycznych Fortran jest prawie nie do pokonania, a C ++ wymaga wiele wysiłku, aby uzyskać podobną prędkość. Można to zrobić i jest kilka bibliotek demonstrujących to (zwykle przy użyciu szablonów wyrażeń), ale nie jest to trywialne i nie dzieje się tak po prostu .

Nie wiem konkretnie o implementacji BLAS-a, ale istnieją bardziej wydajne alogorytmy dla mnożenia macierzy, które mają lepszą złożoność niż O (n3). Dobrze znanym jest Algorytm Strassena

Większość argumentów na drugie pytanie - asembler, dzielenie na bloki itp. (Ale nie mniej niż algorytmy N ^ 3, są one naprawdę nadmiernie rozbudowane) - odgrywa rolę. Ale niska prędkość twojego algorytmu jest spowodowana głównie rozmiarem macierzy i niefortunnym ułożeniem trzech zagnieżdżonych pętli. Twoje macierze są tak duże, że nie mieszczą się od razu w pamięci podręcznej. Możesz przestawić pętle w taki sposób, aby jak najwięcej było zrobionych w wierszu w pamięci podręcznej, w ten sposób radykalnie zmniejszając odświeżanie pamięci podręcznej (przy okazji dzielenie na małe bloki ma efekt analogowy, najlepiej, jeśli pętle nad blokami są ułożone podobnie). Następuje modelowa implementacja macierzy kwadratowych. Na moim komputerze jego zużycie czasu wyniosło około 1:10 w porównaniu ze standardową implementacją (taką jak Twoja). Innymi słowy: nigdy nie programuj mnożenia macierzy wzdłuż "

    void vector(int m, double ** a, double ** b, double ** c) {
      int i, j, k;
      for (i=0; i<m; i++) {
        double * ci = c[i];
        for (k=0; k<m; k++) ci[k] = 0.;
        for (j=0; j<m; j++) {
          double aij = a[i][j];
          double * bj = b[j];
          for (k=0; k<m; k++)  ci[k] += aij*bj[k];
        }
      }
    }

Jeszcze jedna uwaga: ta implementacja jest nawet lepsza na moim komputerze niż zastąpienie wszystkich przez procedurę BLAS cblas_dgemm (wypróbuj ją na swoim komputerze!). Ale znacznie szybsze (1: 4) jest bezpośrednie wywołanie dgemm_ z biblioteki Fortran. Myślę, że ta procedura nie jest w rzeczywistości Fortranem, ale kodem assemblera (nie wiem, co jest w bibliotece, nie mam źródeł). Zupełnie niejasne jest dla mnie, dlaczego cblas_dgemm nie jest tak szybki, skoro według mojej wiedzy jest tylko opakowaniem dla dgemm_.

To realistyczne przyspieszenie. Aby zapoznać się z przykładem tego, co można zrobić za pomocą asemblera SIMD w kodzie C ++, zobacz kilka przykładowych funkcji macierzy iPhone'a - były one ponad 8 razy szybsze niż wersja C i nie są nawet „zoptymalizowane” montażem - nie ma jeszcze wykładania rur i tam to niepotrzebne operacje na stosie.

Również twój kod nie jest „ ogranicz poprawny ” - skąd kompilator wie, że modyfikując C, nie modyfikuje A i B?

W odniesieniu do oryginalnego kodu w mnożeniu MM, odwołanie do pamięci dla większości operacji jest główną przyczyną złej wydajności. Pamięć działa 100-1000 razy wolniej niż pamięć podręczna.

Większość przyspieszenia wynika z zastosowania technik optymalizacji pętli dla tej funkcji potrójnej pętli w mnożeniu MM. Stosowane są dwie główne techniki optymalizacji pętli; rozwijanie i blokowanie. Jeśli chodzi o rozwijanie, rozwijamy dwie zewnętrzne najbardziej pętle i blokujemy je w celu ponownego wykorzystania danych w pamięci podręcznej. Odwijanie pętli zewnętrznej pomaga tymczasowo zoptymalizować dostęp do danych, zmniejszając liczbę odwołań do pamięci do tych samych danych w różnym czasie podczas całej operacji. Zablokowanie indeksu pętli pod określonym numerem pomaga w zachowaniu danych w pamięci podręcznej. Możesz wybrać optymalizację dla pamięci podręcznej L2 lub L3.

https://en.wikipedia.org/wiki/Loop_nest_optimization

Z wielu powodów.

Po pierwsze, kompilatory Fortrana są wysoce zoptymalizowane, a język pozwala im takimi być. C i C ++ są bardzo luźne pod względem obsługi tablic (np. Przypadek wskaźników odnoszących się do tego samego obszaru pamięci). Oznacza to, że kompilator nie może z góry wiedzieć, co robić i jest zmuszony do utworzenia kodu ogólnego. W Fortranie sprawy są bardziej uproszczone, a kompilator ma lepszą kontrolę nad tym, co się dzieje, co pozwala mu na większą optymalizację (np. Przy użyciu rejestrów).

Inną rzeczą jest to, że Fortran przechowuje dane kolumnowo, podczas gdy C przechowuje dane wierszowo. Nie sprawdziłem kodu, ale uważaj na sposób wykonywania produktu. W C musisz skanować mądrze wierszami: w ten sposób skanujesz swoją tablicę wzdłuż ciągłej pamięci, zmniejszając błędy pamięci podręcznej. Brak pamięci podręcznej jest pierwszym źródłem nieefektywności.

Po trzecie, zależy to od używanej implementacji blas. Niektóre implementacje mogą być napisane w asemblerze i zoptymalizowane dla konkretnego używanego procesora. Wersja netlib jest napisana w Fortran 77.

Ponadto wykonujesz wiele operacji, większość z nich jest powtarzana i zbędna. Wszystkie te mnożenia w celu uzyskania indeksu są szkodliwe dla wydajności. Naprawdę nie wiem, jak to się robi w BLAS, ale jest wiele sztuczek, aby zapobiec kosztownym operacjom.

Na przykład możesz w ten sposób przerobić swój kod

template<class ValT>
void mmult(const ValT* A, int ADim1, int ADim2, const ValT* B, int BDim1, int BDim2, ValT* C)
{
if ( ADim2!=BDim1 ) throw std::runtime_error("Error sizes off");

memset((void*)C,0,sizeof(ValT)*ADim1*BDim2);
int cc2,cc1,cr1, a1,a2,a3;
for ( cc2=0 ; cc2<BDim2 ; ++cc2 ) {
    a1 = cc2*ADim2;
    a3 = cc2*BDim1
    for ( cc1=0 ; cc1<ADim2 ; ++cc1 ) {
          a2=cc1*ADim1;
          ValT b = B[a3+cc1];
          for ( cr1=0 ; cr1<ADim1 ; ++cr1 ) {
                    C[a1+cr1] += A[a2+cr1]*b;
           }
     }
  }
} 

Spróbuj, na pewno coś uratujesz.

Jeśli chodzi o twoje pytanie nr 1, powodem jest to, że mnożenie macierzy skaluje się jako O (n ^ 3), jeśli używasz trywialnego algorytmu. Istnieją algorytmy, które skalują się znacznie lepiej .