prestazioni degradazione della moltiplicazione di matrici di single vs array doppia precisione su una macchina multi-core

Question

UPDATE

Purtroppo, a causa della mia svista, ho avuto una versione precedente di MKL (11.1) collegato contro NumPy. La versione più recente di MKL (11.3.1) offre le stesse prestazioni in C e quando viene chiamata da python.

Ciò che oscurava le cose, era anche se collegava le librerie condivise compilate esplicitamente con il più recente MKL, e indicando le variabili LD_ * a loro, e poi in python facendo import numpy, stava facendo in modo che Python chiamasse le vecchie librerie MKL. Solo rimpiazzando nella cartella lib python tutti libmkl _ *. Quindi con MKL più recente ero in grado di eguagliare le prestazioni nelle chiamate python e C.

Informazioni su sfondo/libreria.

La moltiplicazione della matrice è stata effettuata tramite chiamate della libreria MKL di Intel (single-precision) e dgemm (doppia precisione), tramite la funzione numpy.dot. L'effettiva chiamata delle funzioni della libreria può essere verificata con ad es. oprof.

Utilizzando qui CPU core 2x18 E5-2699 v3, quindi un totale di 36 core fisici. KMP_AFFINITY = scatter. Funzionando su Linux.

TL; DR

1) Perché numpy.dot, anche se si sta chiamando le stesse funzioni di libreria MKL, due volte più lento nella migliore rispetto a C codice compilato?

2) Perché tramite numpy.dot si ottiene una riduzione delle prestazioni con l'aumento del numero di core, mentre lo stesso effetto non si osserva nel codice C (chiamando le stesse funzioni della libreria).

Il problema

ho osservato che facendo la moltiplicazione di matrici di singola doppia precisione/galleggia in numpy.dot, così come chiamare cblas_sgemm/dgemm direttamente da un C libreria condivisa compilato dare notevolmente peggiore prestazioni rispetto a chiamare le stesse funzioni MKL cblas_sgemm/dgemm dall'interno del puro codice 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 

facendo esattamente la stessa di cui sopra, ma con doppia precisione A, B e C, si ottiene: 3 nuclei: 20s, 6 core: 10s, 12 core: 5s, 18 core: 4.3s, 24 core: 3, 30 core: 2,8 s, 36 core: 2,8 s.

Il rabbocco di velocità per punti flottanti a precisione singola sembra essere associato a errori di cache. Per 28 core run, ecco l'output di perf. Per precisione singola:

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

E doppia precisione:

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

C libreria condivisa, compilato con

/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); 
} 

funzione wrapper Python, chiamando il sopra libreria compilata:

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) 

Tuttavia, le chiamate esplicite da un binario C-compilato che chiama MKL's cblas_sgemm/cblas_dgemm, con le matrici allocate tramite malloc in C, offrono prestazioni quasi 2x rispetto al codice python, ovvero la chiamata numpy.dot. Inoltre, l'effetto del degrado delle prestazioni con un numero crescente di core NON viene osservato. La migliore prestazione è stata di 900 ms per la moltiplicazione di matrice a precisione singola ed è stata ottenuta utilizzando tutti i 36 core fisici tramite mkl_set_num_cores ed eseguendo il codice C con numactl --interleave = all.

Forse qualche strumento o consiglio di fantasia per profilare/ispezionare/comprendere ulteriormente questa situazione? Anche qualsiasi materiale di lettura è molto apprezzato.

UPDATE Seguendo il consiglio @Hristo Iliev, in esecuzione numactl --interleave = tutti ./ipython non ha modificato i tempi (entro rumore), ma migliora i puri C runtime binari.

Answer 1

Sospetto che ciò sia dovuto a una sfortunata programmazione dei thread. Sono stato in grado di riprodurre un effetto simile al tuo. Python era in esecuzione a ~ 2,2 s, mentre la versione C mostrava enormi variazioni da 1.4-2.2 s.

Applicazione: KMP_AFFINITY=scatter,granularity=thread Ciò garantisce che i 28 thread siano sempre in esecuzione sullo stesso thread del processore.

Riduce entrambi i tempi di esecuzione a ~ 1.24 s più stabili per C e ~ 1.26 s per python.

Si tratta di un sistema Xeon E5-2680 v3 a 28 socket dual core.

È interessante notare che su un sistema Haswell a 24 socket dual core molto simile, sia Python sia C sono quasi identici anche senza affinità/pinning dei thread.

Perché Python influisce sulla pianificazione? Suppongo che ci sia più ambiente di runtime attorno ad esso. La linea di fondo è, senza bloccare la tua performance, i risultati saranno non deterministici.

Inoltre, è necessario considerare che il runtime Intel OpenMP genera un ulteriore thread di gestione che può confondere lo scheduler. Esistono più opzioni per il blocco, ad esempio KMP_AFFINITY=compact, ma per qualche motivo è completamente incasinato sul mio sistema. È possibile aggiungere ,verbose alla variabile per vedere come il runtime sta bloccando i thread.

likwid-pin è un'alternativa utile che offre un controllo più conveniente.

In generale, la precisione singola deve essere almeno altrettanto veloce della doppia precisione.La doppia precisione può essere più lenta perché:

• È necessaria più memoria/larghezza di banda della cache per la doppia precisione.
• È possibile creare ALU con througput più elevato per precisione singola, ma che di solito non si applica alle CPU ma piuttosto alle GPU.

Penso che una volta eliminata l'anomalia delle prestazioni, ciò si rifletterà nei numeri.

Quando si scala il numero di thread per MKL/* GEMM, prendere in considerazione la larghezza di banda/cache condivisa

• memoria può diventare un collo di bottiglia, che limita la scalabilità
modalità
• Turbo farà diminuire efficacemente la frequenza di core quando utilizzo crescente. Questo si applica anche quando si corre alla frequenza nominale: sui processori Haswell-EP, le istruzioni AVX imporranno una "frequenza di base AVX" inferiore - ma il processore può eccedere quello quando vengono utilizzati meno core/headroom termico disponibile e in generale anche più per un breve periodo. Se vuoi risultati perfettamente neutri, dovresti usare la frequenza di base AVX, che è 1.9 GHz per te. È documentato here e spiegato in one picture.

Non penso che ci sia un modo molto semplice per misurare il modo in cui la vostra applicazione è influenzata da una cattiva programmazione. Puoi esporlo con perf trace -e sched:sched_switch e c'è lo some software per visualizzarlo, ma questo avrà una curva di apprendimento elevata. E poi ancora: per l'analisi parallela delle prestazioni dovresti avere i thread bloccati comunque.