Wydajność degradacja mnożenie macierzy pojedynczy vs podwójnych tablic dokładność na maszynie wielordzeniowej

Question

UPDATE

Niestety, ze względu na moje niedopatrzenie, miałem starszą wersję MKL (11,1) umieszczonego przed numpy. Nowsza wersja MKL (11.3.1) daje taką samą wydajność w C i gdy jest wywoływana z Pythona.

To, co było ukryte, było nawet wtedy, gdy łączenie kompilowanych bibliotek współdzielonych jawnie z nowszym MKL i wskazywanie zmiennych LD_ * na nich, a następnie w pythonie podczas importowania numpy, w jakiś sposób powodowało, że python wywoływał stare biblioteki MKL. Tylko przez zastąpienie w folderze Pythona lib wszystkich libmkl _ *., Tak więc przy nowszym MKL byłem w stanie dopasować wydajność w pythonie i wywołaniach C.

Informacje o tle/bibliotece.

Mnożenie macierzy zostało wykonane za pomocą sgemm (pojedyncza precyzja) i dgemm (podwójna precyzja) Wywołania bibliotek MKL Intela, za pomocą funkcji numpy.dot. Faktyczne wywołanie funkcji biblioteki można zweryfikować za pomocą np. oprof.

Tutaj używa się procesora 2x18 rdzeniowego E5-2699 v3, a więc łącznie 36 rdzeni fizycznych. KMP_AFFINITY = scatter. Działa na Linuksie.

TL; DR

1) Dlaczego numpy.dot, mimo że dzwoni te same funkcje biblioteczne MKL, dwa razy wolniej w porównaniu do C najlepiej skompilowanego kodu?

2) Dlaczego poprzez numpy.dot wydajność spada wraz ze wzrostem liczby rdzeni, podczas gdy ten sam efekt nie jest przestrzegany w kodzie C (wywoływanie tych samych funkcji bibliotecznych).

Problem

ja zauważył, że robi mnożenie macierzy pojedynczej podwójnej precyzji/pływa w numpy.dot, jak również wywołanie cblas_sgemm/dgemm bezpośrednio od skompilowany C udostępnionej biblioteki w daje zauważalnie gorzej wydajność w porównaniu do wywoływania funkcji MKL cblas_sgemm/dgemm z wewnątrz czystego kodu C.

import numpy as np 
import mkl 
n = 10000 
A = np.random.randn(n,n).astype('float32') 
B = np.random.randn(n,n).astype('float32') 
C = np.zeros((n,n)).astype('float32') 

mkl.set_num_threads(3); %time np.dot(A, B, out=C) 
11.5 seconds 
mkl.set_num_threads(6); %time np.dot(A, B, out=C) 
6 seconds 
mkl.set_num_threads(12); %time np.dot(A, B, out=C) 
3 seconds 
mkl.set_num_threads(18); %time np.dot(A, B, out=C) 
2.4 seconds 
mkl.set_num_threads(24); %time np.dot(A, B, out=C) 
3.6 seconds 
mkl.set_num_threads(30); %time np.dot(A, B, out=C) 
5 seconds 
mkl.set_num_threads(36); %time np.dot(A, B, out=C) 
5.5 seconds 

Robi dokładnie to samo co powyżej, ale z podwójną precyzją A, B i C, można uzyskać: 3 rdzenie 20s, 6 rdzeni: 10s, 12 żył: 5s, 18 żył: 4.3s, 24 rdzenie: 3s, 30 rdzeni: 2,8s, 36 rdzeni: 2,8s.

Uzupełnienie prędkości dla punktów pływających o pojedynczej precyzji wydaje się być związane z chybieniami w pamięci podręcznej. Dla 28 rdzeń, tutaj jest wyjście z perf. Dla pojedynczej precyzji:

perf stat -e task-clock,cycles,instructions,cache-references,cache-misses ./ptestf.py 
631,301,854 cache-misses # 31.478 % of all cache refs 

i podwójną precyzją:

93,087,703 cache-misses # 5.164 % of all cache refs 

C biblioteki współdzielonej, sporządzoną z

/opt/intel/bin/icc -o comp_sgemm_mkl.so -openmp -mkl sgem_lib.c -lm -lirc -O3 -fPIC -shared -std=c99 -vec-report1 -xhost -I/opt/intel/composer/mkl/include 

#include <stdio.h> 
#include <stdlib.h> 
#include "mkl.h" 

void comp_sgemm_mkl(int m, int n, int k, float *A, float *B, float *C); 

void comp_sgemm_mkl(int m, int n, int k, float *A, float *B, float *C) 
{ 
    int i, j; 
    float alpha, beta; 
    alpha = 1.0; beta = 0.0; 

    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 
       m, n, k, alpha, A, k, B, n, beta, C, n); 
} 

Pythona funkcji otoki, nazywając powyższy skompilowaną bibliotekę:

def comp_sgemm_mkl(A, B, out=None): 
    lib = CDLL(omplib) 
    lib.cblas_sgemm_mkl.argtypes = [c_int, c_int, c_int, 
           np.ctypeslib.ndpointer(dtype=np.float32, ndim=2), 
           np.ctypeslib.ndpointer(dtype=np.float32, ndim=2), 
           np.ctypeslib.ndpointer(dtype=np.float32, ndim=2)] 
    lib.comp_sgemm_mkl.restype = c_void_p 
    m = A.shape[0] 
    n = B.shape[0] 
    k = B.shape[1] 
    if np.isfortran(A): 
     raise ValueError('Fortran array') 
    if m != n: 
     raise ValueError('Wrong matrix dimensions') 
    if out is None: 
     out = np.empty((m,k), np.float32) 
    lib.comp_sgemm_mkl(m, n, k, A, B, out) 

Jednak jawne wywołania z C skompilowanego binarnie wywołującego MKL cblas_sgemm/cblas_dgemm, z tablicami przydzielonymi przez malloc w C, dają prawie 2x lepszą wydajność w porównaniu do kodu Pythona, tj. Wywołania numpy.dot. Nie obserwuje się również wpływu degradacji wydajności przy rosnącej liczbie rdzeni. Najlepsza wydajność wyniosła 900 ms dla pojedynczej precyzji mnożenia macierzy i została osiągnięta przy użyciu wszystkich 36 rdzeni fizycznych za pośrednictwem mkl_set_num_cores i uruchomienia kodu C za pomocą numactl --interleave = all.

Być może jakieś wymyślne narzędzia lub porady dotyczące profilowania/kontroli/zrozumienia tej sytuacji dalej? Każdy materiał do czytania jest również bardzo doceniany.

UPDATE Za radą @Hristo Iliev, bieganie numactl --interleave = all ./ipython nie zmieniać czasy (w promieniu hałasu), ale poprawia czystym C czasy pracy binarnych.

Answer 1

Podejrzewam, że jest to spowodowane niefortunnym planowaniem wątków. Udało mi się odtworzyć efekt podobny do twojego. Python działał w czasie ~ 2.2 s, podczas gdy wersja C wyświetlała ogromne odmiany z 1,4-2,2 s.

Zastosowanie: KMP_AFFINITY=scatter,granularity=thread Zapewnia to, że 28 wątków jest zawsze uruchomionych na tym samym wątku procesora.

Zmniejsza obie czasy pracy do stabilniejszej ~ 1.24 s dla C i ~ 1.26 s dla pythona.

Jest to system 28-rdzeniowy z podwójnym gniazdem Xeon E5-2680 v3.

Co ciekawe, na bardzo podobnym dwurdzeniowym dwurdzeniowym systemie Haswell, zarówno python, jak i C działają niemal identycznie, nawet bez powinowactwa/przypinania nici.

Dlaczego python wpływa na harmonogram? Zakładam, że wokół niego jest więcej środowiska uruchomieniowego. Najważniejsze jest to, że bez przypinania wyników wyniki będą niedeterministyczne.

Należy również wziąć pod uwagę, że środowisko wykonawcze Intel OpenMP spawns dodatkowy wątek zarządzania, który może zmylić program planujący. Istnieje więcej opcji przypinania, na przykład KMP_AFFINITY=compact - ale z jakiegoś powodu jest to całkowicie pomieszane w moim systemie. Do zmiennej można dodać ,verbose, aby zobaczyć, jak środowisko wykonawcze przypina wątki.

likwid-pin to przydatna alternatywa zapewniająca wygodniejszą kontrolę.

Zwykle pojedyncza precyzja powinna być co najmniej tak duża, jak podwójna precyzja.Podwójna precyzja może być wolniejsza, ponieważ:

• Potrzebujesz większej przepustowości pamięci/pamięci podręcznej dla podwójnej precyzji.
• Można zbudować jednostki ALU, które mają wyższą przepustowość dla pojedynczej precyzji, ale zwykle nie dotyczy to procesorów, ale procesorów graficznych.

Myślę, że gdy pozbędziesz się anomalii wydajności, będzie to odzwierciedlone w twoich liczbach.

Podczas skalowania liczbę wątków dla MKL/* GEMM rozważyć

• Przepustowość pamięci/pamięci podręcznej może stać się wąskim gardłem, ograniczając skalowalność
tryb
• Turbo będzie skutecznie zmniejszyć częstotliwość rdzenia gdy jest zwiększenie wykorzystania. Odnosi się to nawet w przypadku pracy z częstotliwością nominalną: w procesorach Haswell-EP instrukcje AVX będą narzucać niższą "częstotliwość bazową AVX" - ale procesor może przekroczyć to, gdy jest używanych mniej rdzeni/dostępna jest przestrzeń cieplna, a nawet nawet więcej przez krótki czas. Jeśli chcesz uzyskać idealnie neutralne wyniki, będziesz musiał użyć częstotliwości bazowej AVX, która wynosi 1,9 GHz. Jest to udokumentowane here i wyjaśnione w one picture.

Nie sądzę, że istnieje naprawdę prosty sposób pomiaru wpływu złego planowania na aplikację. Możesz to wyrazić przy pomocy perf trace -e sched:sched_switch i jest tam some software, aby to sobie wyobrazić, ale będzie to miało wysoką krzywą uczenia się. I jeszcze raz - w celu równoległej analizy wydajności i tak należy przypiąć wątki.