Tego rodzaju pytania powracają i należy na nie odpowiedzieć jaśniej niż „MATLAB używa wysoce zoptymalizowanych bibliotek” lub „MATLAB używa MKL” raz na przepełnieniu stosu.
Historia:
Mnożenie macierzy (wraz z wektorem macierzy, mnożeniem wektora i wieloma rozkładami macierzy) jest (są) najważniejszymi problemami w algebrze liniowej. Inżynierowie rozwiązują te problemy z komputerami od pierwszych dni.
Nie jestem ekspertem od historii, ale najwyraźniej wtedy wszyscy przepisali jego wersję FORTRAN za pomocą prostych pętli. Potem pojawiła się pewna standaryzacja, wraz z identyfikacją „jąder” (podstawowych procedur), które większość problemów algebry liniowej potrzebowała do rozwiązania. Te podstawowe operacje zostały następnie znormalizowane w specyfikacji o nazwie: podprogramy podstawowej algebry liniowej (BLAS). Inżynierowie mogliby wówczas nazwać te standardowe, dobrze przetestowane procedury BLAS w swoim kodzie, co znacznie ułatwi ich pracę.
BLAS:
BLAS ewoluował z poziomu 1 (pierwsza wersja, która definiowała operacje na wektorach skalarnych i wektorach) do poziomu 2 (operacje na macierzach wektorowych) do poziomu 3 (operacje na macierzach) i dostarczał coraz więcej „jąder”, więc bardziej znormalizowanych i więcej podstawowych operacji algebry liniowej. Oryginalne implementacje FORTRAN 77 są nadal dostępne na stronie internetowej Netlib .
W kierunku lepszej wydajności:
Tak więc z biegiem lat (zwłaszcza między wydaniami BLAS poziomu 1 i poziomu 2: wczesne lata 80.) zmienił się sprzęt wraz z nadejściem operacji wektorowych i hierarchii pamięci podręcznej. Te zmiany umożliwiły znaczne zwiększenie wydajności podprogramów BLAS. Następnie pojawili się różni dostawcy, którzy wdrażali procedury BLAS, które były coraz bardziej wydajne.
Nie znam wszystkich historycznych wdrożeń (wtedy nie urodziłem się ani nie byłem dzieckiem), ale na początku 2000 roku pojawiły się dwa najbardziej znaczące: Intel MKL i GotoBLAS. Twój Matlab korzysta z Intel MKL, który jest bardzo dobrym, zoptymalizowanym BLAS, a to wyjaśnia doskonałą wydajność, którą widzisz.
Szczegóły techniczne dotyczące mnożenia macierzy:
Dlaczego więc Matlab (MKL) jest tak szybki w dgemm
(dwukrotne ogólne mnożenie macierzy-macierzy)? Mówiąc prosto: ponieważ wykorzystuje wektoryzację i dobre buforowanie danych. W bardziej złożonych terminach: patrz artykuł dostarczony przez Jonathana Moore'a.
Zasadniczo, gdy wykonujesz mnożenie w dostarczonym kodzie C ++, wcale nie jesteś przyjazny dla pamięci podręcznej. Ponieważ podejrzewam, że stworzyłeś tablicę wskaźników do tablic wierszowych, twoje dostępy w wewnętrznej pętli do k-tej kolumny „matice2”: matice2[m][k]
są bardzo wolne. Rzeczywiście, kiedy uzyskujesz dostęp matice2[0][k]
, musisz uzyskać k-ty element tablicy 0 macierzy. Następnie w następnej iteracji musisz uzyskać dostęp matice2[1][k]
, który jest k-tym elementem innej tablicy (tablica 1). Następnie w następnej iteracji uzyskujesz dostęp do kolejnej tablicy i tak dalej ... Ponieważ cała matryca matice2
nie mieści się w najwyższych pamięciach podręcznych (ma 8*1024*1024
duże bajty), program musi pobrać żądany element z pamięci głównej, tracąc dużo czas.
Jeśli właśnie przetransponowałeś macierz, aby dostęp był w ciągłych adresach pamięci, twój kod działałby już znacznie szybciej, ponieważ teraz kompilator może ładować całe wiersze w pamięci podręcznej w tym samym czasie. Po prostu wypróbuj zmodyfikowaną wersję:
timer.start();
float temp = 0;
//transpose matice2
for (int p = 0; p < rozmer; p++)
{
for (int q = 0; q < rozmer; q++)
{
tempmat[p][q] = matice2[q][p];
}
}
for(int j = 0; j < rozmer; j++)
{
for (int k = 0; k < rozmer; k++)
{
temp = 0;
for (int m = 0; m < rozmer; m++)
{
temp = temp + matice1[j][m] * tempmat[k][m];
}
matice3[j][k] = temp;
}
}
timer.stop();
Możesz więc zobaczyć, jak po prostu pamięć podręczna znacznie zwiększyła wydajność kodu. Teraz prawdziwe dgemm
implementacje wykorzystują to na bardzo szerokim poziomie: wykonują mnożenie na blokach macierzy określonych przez rozmiar TLB (bufor lookaside tłumaczenia, krótka historia: co można skutecznie buforować), aby przesyłać strumieniowo do procesora dokładnie ilość danych, które może przetworzyć. Innym aspektem jest wektoryzacja, używają one wektoryzowanych instrukcji procesora w celu uzyskania optymalnej przepustowości instrukcji, czego tak naprawdę nie można zrobić z wieloplatformowego kodu C ++.
Wreszcie, ludzie twierdzący, że jest to spowodowane algorytmem Strassena lub Coppersmitha-Winograda, są w błędzie, oba te algorytmy nie są możliwe do wdrożenia w praktyce ze względu na wspomniane wyżej względy sprzętowe.