tempo di gemm delle cblas dipende dai valori della matrice di input – Ubuntu 14.04

Questa è un’estensione della mia precedente domanda , ma la sto chiedendo separatamente perché mi sto davvero frustrando, quindi per favore non votarlo!

Domanda: quale potrebbe essere il motivo alla base di una chiamata cblas_sgemm che impiega molto meno tempo per le matrici con un numero elevato di zeri rispetto alla stessa chiamata cblas_sgemm per le matrici dense?

So che gemv è progettato per la moltiplicazione di matrice vettoriale, ma perché non posso usare gemm per la moltiplicazione della matrice vettoriale se impiega meno tempo, specialmente per le matrici sparse?

Di seguito è riportato un breve codice rappresentativo. Chiede di inserire un valore e quindi popola un vettore con quel valore. Sostituisce quindi ogni 32 ° valore con il suo indice. Quindi, se inseriamo ‘0’, otteniamo un vettore sparse ma per altri valori otteniamo un vettore denso.

#include  #include  #include  #include  using namespace std; int main() { const int m = 5000; timespec blas_start, blas_end; long totalnsec; //total nano sec double totalsec, totaltime; int i, j; float *A = new float[m]; // 1 xm float *B = new float[m*m]; // mxm float *C = new float[m]; // 1 xm float input; cout <> input; // enter 0 for sparse // input martix A: every 32nd element is non-zero, rest of the values = input for(i = 0; i < m; i++) { A[i] = input; if( i % 32 == 0) //adjust for sparsity A[i] = i; } // input matrix B: identity matrix for(i = 0; i < m; i++) for(j = 0; j < m; j++) B[i*m + j] = (i==j); clock_gettime(CLOCK_REALTIME, &blas_start); cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1.0f, A, m, B, m, 0.0f, C, m); clock_gettime(CLOCK_REALTIME, &blas_end); /* for(i = 0; i < m; i++) printf("%f ", C[i]); printf("\n\n"); */ // Print time totalsec = (double)blas_end.tv_sec - (double)blas_start.tv_sec; totalnsec = blas_end.tv_nsec - blas_start.tv_nsec; if(totalnsec < 0) { totalnsec += 1e9; totalsec -= 1; } totaltime = totalsec + (double)totalnsec*1e-9; cout<<"Duration = "<< totaltime << "\n"; return 0; } 

Lo eseguo come segue in Ubuntu 14.04 con Blas 3.0

     [email protected]:~/uas/stackoverflow$ g++ gemmcomp.cpp -o gemmcomp.o -lblas [email protected]:~/uas/stackoverflow$ ./gemmcomp.o Enter a value to populate the vector (0 for sparse) 5 Duration = 0.0291558 [email protected]:~/uas/stackoverflow$ ./gemmcomp.o Enter a value to populate the vector (0 for sparse) 0 Duration = 0.000959521 

    MODIFICARE

    Se sostituisco la chiamata gemm con la seguente chiamata gemv, la sparsità della matrice non ha importanza

     //cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1.0f, A, m, B, m, 0.0f, C, m); cblas_sgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0f, B, m, A, 1, 0.0f, C, 1); 

    risultati

     [email protected]:~/uas/stackoverflow$ g++ gemmcomp.cpp -o gemmcomp.o -lblas [email protected]:~/uas/stackoverflow$ ./gemmcomp.o Enter a value to populate the vector (0 for sparse) 5 Duration = 0.0301581 [email protected]:~/uas/stackoverflow$ ./gemmcomp.o Enter a value to populate the vector (0 for sparse) 0 Duration = 0.0299282 

    Ma il problema è che sto cercando di ottimizzare il codice di qualcun altro usando cublas e lui sta usando con successo ed efficientemente gemm per eseguire questa moltiplicazione di matrice vettoriale. Quindi, devo metterlo alla prova o dimostrare categoricamente che questa chiamata sia errata

    MODIFICARE

    Ho persino aggiornato la mia libreria blas oggi usando

     sudo apt-get install libblas-dev liblapack-dev 

    EDIT: eseguito i seguenti comandi come suggerito da diversi contributori

     [email protected]:~/uas/stackoverflow$ ll -d /usr/lib/libblas* /etc/alternatives/libblas.* lrwxrwxrwx 1 root root 26 مارچ 13 2015 /etc/alternatives/libblas.a -> /usr/lib/libblas/libblas.a lrwxrwxrwx 1 root root 27 مارچ 13 2015 /etc/alternatives/libblas.so -> /usr/lib/libblas/libblas.so lrwxrwxrwx 1 root root 29 مارچ 13 2015 /etc/alternatives/libblas.so.3 -> /usr/lib/libblas/libblas.so.3 lrwxrwxrwx 1 root root 29 مارچ 13 2015 /etc/alternatives/libblas.so.3gf -> /usr/lib/libblas/libblas.so.3 drwxr-xr-x 2 root root 4096 مارچ 13 2015 /usr/lib/libblas/ lrwxrwxrwx 1 root root 27 مارچ 13 2015 /usr/lib/libblas.a -> /etc/alternatives/libblas.a lrwxrwxrwx 1 root root 28 مارچ 13 2015 /usr/lib/libblas.so -> /etc/alternatives/libblas.so lrwxrwxrwx 1 root root 30 مارچ 13 2015 /usr/lib/libblas.so.3 -> /etc/alternatives/libblas.so.3 lrwxrwxrwx 1 root root 32 مارچ 13 2015 /usr/lib/libblas.so.3gf -> /etc/alternatives/libblas.so.3gf [email protected]:~/uas/stackoverflow$ ldd ./gemmcomp.o linux-gate.so.1 => (0xb76f6000) libblas.so.3 => /usr/lib/libblas.so.3 (0xb765e000) libstdc++.so.6 => /usr/lib/i386-linux-gnu/libstdc++.so.6 (0xb7576000) libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb73c7000) libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7381000) /lib/ld-linux.so.2 (0xb76f7000) libgcc_s.so.1 => /lib/i386-linux-gnu/libgcc_s.so.1 (0xb7364000) 

    Domanda: quale potrebbe essere il motivo alla base di una chiamata cblas_sgemm che impiega molto meno tempo per le matrici con un numero elevato di zeri rispetto alla stessa chiamata cblas_sgemm per le matrici dense?

    Sembra che l’implementazione BLAS fornita dal pacchetto libblas-dev predefinito per Ubuntu 14.04 (e probabilmente altre distribuzioni di Ubuntu) includa un’ottimizzazione per i casi in cui alcuni elementi della matrice sono zero.

    Per Ubuntu 14.04, il codice sorgente per l’implementazione / pacchetto BLAS (e cblas) può essere scaricato da qui .

    Dopo aver decompresso l’archivio, abbiamo una cblas/src che contiene l’API cblas e abbiamo un’altra directory src che contiene le implementazioni F77 di varie routine blas.

    Nel caso di cblas_sgemm, quando viene specificato il parametro CblasRowMajor , il cblas/src/cblas_sgemm.c chiamerà la routine fortran sottostante come segue:

     void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const float alpha, const float *A, const int lda, const float *B, const int ldb, const float beta, float *C, const int ldc) { ... } else if (Order == CblasRowMajor) ... F77_sgemm(F77_TA, F77_TB, &F77_N, &F77_M, &F77_K, &alpha, B, &F77_ldb, A, &F77_lda, &beta, C, &F77_ldc); 

    Si noti che per questa chiamata principale di riga, l’ordine delle matrici A e B viene invertito quando viene passato alla routine F77_sgemm . Questo è ragionevole ma non approfondirò perché qui. È sufficiente notare che A è diventato B nella chiamata / codice di Fortran e B è diventato A

    Quando ispezioniamo la routine fortran corrispondente in src/sgemm.f , vediamo la seguente sequenza di codice:

     * * Start the operations. * IF (NOTB) THEN IF (NOTA) THEN * * Form C := alpha*A*B + beta*C. * DO 90 J = 1,N IF (BETA.EQ.ZERO) THEN DO 50 I = 1,M C(I,J) = ZERO 50 CONTINUE ELSE IF (BETA.NE.ONE) THEN DO 60 I = 1,M C(I,J) = BETA*C(I,J) 60 CONTINUE END IF DO 80 L = 1,K IF (B(L,J).NE.ZERO) THEN ***OPTIMIZATION TEMP = ALPHA*B(L,J) DO 70 I = 1,M C(I,J) = C(I,J) + TEMP*A(I,L) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE 

    Quanto sopra è la sezione di codice che gestisce il caso in cui non viene indicata alcuna trasposizione di A e nessuna trasposizione di B (che è vero per questo caso di test di riga principale di cblas). L’operazione di moltiplicazione di righe / colonne di matrice viene gestita nei loop che iniziano dove ho aggiunto la nota ***OPTIMIZATION . In particolare, se l’elemento matrice B(L,J) è zero, allora la chiusura DO-loop alla linea 70 viene saltata. Ma ricorda che B corrisponde alla matrice A passata alla routine cblas_sgemm .

    Il saltare di questo do-loop consente alla funzione sgemm implementata in questo modo di essere sostanzialmente più veloce per i casi in cui vi è un numero elevato di zeri nella matrice A passati a cblas_sgemm quando viene specificato row-major.

    Sperimentalmente, non tutte le implementazioni di Blas hanno questa ottimizzazione. Test sulla stessa piattaforma ma usando libopenblas-dev invece di libblas-dev non offre tale accelerazione, cioè essenzialmente nessuna differenza di tempo di esecuzione quando la matrice A è per lo più zero, rispetto al caso in cui non lo è.

    Si noti che il codice di fortran (77) qui sembra essere simile o identico alle precedenti versioni pubblicate della routine sgemm.f come qui . Le versioni più recenti di questa routine fortran che ho trovato non contengono questa ottimizzazione, come qui .