Riduci le righe della matrice con CUDA

Windows 7, NVidia GeForce 425M. 

Ho scritto un semplice codice CUDA che calcola le somme di riga di una matrice. La matrice ha una rappresentazione unidimensionale (puntatore a un float).

La versione seriale del codice è sotto (ha 2 loop, come previsto):

 void serial_rowSum (float* m, float* output, int nrow, int ncol) { float sum; for (int i = 0 ; i < nrow ; i++) { sum = 0; for (int j = 0 ; j < ncol ; j++) sum += m[i*ncol+j]; output[i] = sum; } } 

All’interno del codice CUDA, chiamo la funzione del kernel che spazza la matrice per righe. Sotto, lo snippet di chiamata del kernel:

 dim3 threadsPerBlock((unsigned int) nThreadsPerBlock); // has to be multiple of 32 dim3 blocksPerGrid((unsigned int) ceil(nrow/(float) nThreadsPerBlock)); kernel_rowSum<<>>(d_m, d_output, nrow, ncol); 

e la funzione kernel che esegue la sum parallela delle righe (ha ancora 1 ciclo):

 __global__ void kernel_rowSum(float *m, float *s, int nrow, int ncol) { int rowIdx = threadIdx.x + blockIdx.x * blockDim.x; if (rowIdx < nrow) { float sum=0; for (int k = 0 ; k < ncol ; k++) sum+=m[rowIdx*ncol+k]; s[rowIdx] = sum; } } 

Fin qui tutto bene. I risultati seriali e paralleli (CUDA) sono uguali.

Il punto è che la versione CUDA impiega quasi il doppio di quella seriale per il calcolo, anche se cambio il parametro nThreadsPerBlock : ho testato nThreadsPerBlock da 32 a 1024 (il numero massimo di thread per blocco consentito per la mia scheda).

IMO, la dimensione della matrice è abbastanza grande da giustificare la parallelizzazione: 90,000 x 1,000 .

Di seguito, nThreadsPerBlock il tempo trascorso per le versioni seriali e parallele utilizzando diversi nThreadsPerBlock . Tempo riportato in msec su una media di 100 campioni:

Matrice: nrow = 90000 x ncol = 1000

Seriale: tempo medio trascorso per campione in msec ( 100 campioni): 289.18 .

CUDA ( 32 ThreadsPerBlock): tempo medio trascorso per campione in msec ( 100 campioni): 497.11 .

CUDA ( 1024 ThreadsPerBlock): tempo medio trascorso per campione in msec ( 100 campioni): 699.66 .

Per ogni evenienza, la versione con nThreadsPerBlock è la più veloce / lenta.

Capisco che c’è un tipo di sovraccarico durante la copia da host a dispositivo e viceversa, ma forse la lentezza è perché non sto implementando il codice più veloce.

Dal momento che sono lontano dall’essere un esperto CUDA:

Sto codificando la versione più veloce per questa attività? Come posso migliorare il mio codice? Posso sbarazzarmi del ciclo nella funzione del kernel?

Ogni pensiero è apprezzato.

MODIFICA 1

Sebbene descriva una rowSum standard, mi interessa l’operazione AND / OR delle righe che hanno valori (0;1} , come rowAND / rowOR . Detto questo, non mi permette di sfruttare il moltiplicarsi di cuBLAS di 1 ‘s COL vettore a colonna COL , come suggerito da alcuni commentatori.

MODIFICA 2

Come suggerito dagli utenti altri utenti e qui approvato:

DIMENTICARE DI PROVARE A SCRIVERE LE TUE FUNZIONI , usa invece la libreria Thrust e la magia arriva.

Dal momento che hai menzionato hai bisogno dell’algoritmo di riduzione generale diverso dalla sum. Proverò a dare 3 approcci qui. l’approccio kernel può avere le massime prestazioni. l’approccio di spinta è più facile da implementare. L’approccio cuBLAS funziona solo con sum e ha buone prestazioni.

Approccio al kernel

Ecco un ottimo documento che introduce come ottimizzare la riduzione parallela standard. La riduzione standard può essere divisa in 2 fasi.

  1. Blocchi di thread multipli ciascuno riduce una parte dei dati;
  2. Un blocco di thread si riduce dal risultato dello stage 1 all’elemento finale 1.

Per il tuo problema di riduzione multipla (ridurre le file di tappetini), solo la fase 1 è sufficiente. L’idea è di ridurre 1 riga per blocco di thread. Per ulteriori considerazioni come blocco di più righe per thread o 1 riga per blocchi di thread multipli, è ansible fare riferimento alla carta fornita da @Novak . Ciò potrebbe migliorare maggiormente le prestazioni, specialmente per le matrici con cattive condizioni.

Approccio di spinta

La multi-riduzione generale può essere effettuata da thrust::reduction_by_key in pochi minuti. Qui puoi trovare alcune discussioni Determinando l’elemento minimo e la sua posizione in ciascuna colonna della matrice con CUDA Thrust .

Comunque thrust::reduction_by_key non presuppone che ogni riga abbia la stessa lunghezza, quindi otterrai una penalizzazione delle prestazioni. Un altro post Come normalizzare le colonne di matrice in CUDA con prestazioni massime? fornisce un confronto di profiling tra thrust::reduction_by_key e cuBLAS approach sulla sum di righe. Potrebbe darti una comprensione di base della performance.

approccio CuBLAS

La sum di righe / colonne di una matrice A può essere vista come una moltiplicazione di matrice-vettore in cui gli elementi del vettore sono tutti. può essere rappresentato dal seguente codice MATLAB.

 y = A * ones(size(A,2),1); 

dove y è la sum delle righe di A.

La libreria cuBLAS fornisce una funzione di moltiplicazione di matrice vettoriale ad alte prestazioni cublasgemv() per questa operazione.

Il risultato del tempo mostra che questa routine è solo il 10 ~ 50% più lento del semplice leggere tutti gli elementi di A una volta, che può essere visto come il limite superiore teorico della prestazione per questa operazione.

Se questa è la misura (sumndo le righe) delle operazioni che devi fare con questi dati, non mi aspetto un notevole vantaggio dalla GPU. Hai esattamente una operazione aritmetica per elemento dati e per questo stai pagando il costo del trasferimento di quell’elemento dati alla GPU. E al di là di una certa dimensione del problema (qualunque cosa sia necessaria per mantenere la macchina occupata) non si ottiene alcun vantaggio dalle dimensioni dei problemi più grandi, perché l’intensità aritmetica è O (n).

Quindi questo non è un problema particolarmente eccitante da risolvere sulla GPU.

Ma come hanno indicato i talonmie, hai un problema di coalescenza nel modo in cui l’hai realizzato, il che rallenterà ulteriormente le cose. Diamo un’occhiata ad un piccolo esempio:

  C1 C2 C3 C4 R1 11 12 13 14 R2 21 22 23 24 R3 31 32 33 34 R4 41 42 43 44 

Sopra è un semplice esempio pittorico di una piccola porzione della matrice. L’archiviazione dei dati macchina è tale che gli elementi (11), (12), (13) e (14) sono memorizzati in posizioni di memoria adiacenti.

Per l’accesso a coalescenza , vogliamo un modello di accesso tale che le posizioni di memoria adiacenti siano richieste dalla stessa istruzione, eseguite attraverso la distorsione.

Dobbiamo pensare all’esecuzione del codice dal punto di vista di un warp , ovvero 32 thread in esecuzione in lock-step. Cosa sta facendo il tuo codice? Quali elementi sta recuperando (chiedendo) ad ogni passo / istruzione? Diamo un’occhiata a questa riga di codice:

  sum+=m[rowIdx*ncol+k]; 

I fili adiacenti nel warp hanno valori adiacenti (cioè consecutivi) per rowIdx quando hai creato quella variabile. Quindi, quando k = 0, quale elemento di dati viene richiesto da ciascun thread quando proviamo a recuperare il valore m[rowIdx*ncol+k] ?

Nel blocco 0, il thread 0 ha un rowIdx di 0. Il thread 1 ha un rowIdx di 1, ecc. Quindi i valori richiesti da ciascun thread in questa istruzione sono:

 Thread: Memory Location: Matrix Element: 0 m[0] (11) 1 m[ncol] (21) 2 m[2*ncol] (31) 3 m[3*ncol] (41) 

Ma questo non è un accesso coalescente! Gli elementi (11), (21), ecc. Non sono adiacenti in memoria. Per l’accesso a coalescenza, vorremmo che la riga dell’elemento Matrix venisse letta in questo modo:

 Thread: Memory Location: Matrix Element: 0 m[?] (11) 1 m[?] (12) 2 m[?] (13) 3 m[?] (14) 

Se poi lavori all’indietro per determinare qual è il valore di ? dovrebbe essere, ti verrà in mente un’istruzione simile a questa:

  sum+=m[k*ncol+rowIdx]; 

Ciò fornirà un accesso coalescente, ma non ti darà la risposta corretta, perché ora stiamo sumndo le colonne della matrice invece delle righe della matrice. Possiamo risolvere questo problema riorganizzando la memorizzazione dei dati in ordine di colonna principale anziché in ordine di riga principale. (Dovresti essere in grado di google per idee, giusto?) Concettualmente, questo equivale a trasporre la tua matrice m . Se questo è conveniente per te o non lo fa è al di fuori della portata della tua domanda, come la vedo io, e non proprio un problema CUDA. Potrebbe essere una cosa semplice da fare mentre stai creando la matrice sull’host o trasferendo la matrice da un host all’altro. Ma in sintesi, non conosco un modo per sumre le righe della matrice con accesso a coalescenza al 100%, se la matrice è memorizzata in ordine di riga maggiore. (Potresti ricorrere a una sequenza di riduzioni di riga, ma mi sembra doloroso.)

Non è raro che, quando stiamo pensando ai modi per accelerare il codice sulla GPU, prendiamo in considerazione la possibilità di riorganizzare la memorizzazione dei dati per facilitare la GPU. Questo è un esempio.

E, sì, quello che sto delineando qui mantiene ancora un ciclo nel kernel.

Come commento aggiuntivo, suggerirei di sincronizzare separatamente le porzioni di copia dei dati e le parti del kernel (calcolo). Non posso dire dalla tua domanda se stai programmando solo il kernel o l’intera operazione (GPU), incluse le copie dei dati. Se cronometrate le copie dei dati separatamente, potreste scoprire che solo il tempo di copia dei dati supera il tempo della CPU. Qualsiasi tentativo di ottimizzare il codice CUDA non influirà sul tempo di copia dei dati. Questo potrebbe essere un punto dati utile prima di dedicare molto tempo a questo.

La riduzione delle righe di una matrice può essere risolta utilizzando CUDA Thrust in tre modi (potrebbero non essere gli unici, ma affrontare questo punto è fuori ambito). Come riconosciuto anche dallo stesso OP, l’uso di CUDA Thrust è preferibile per un tale tipo di problema. Inoltre, è ansible un approccio usando cuBLAS .

APPROCCIO # 1 – reduce_by_key

Questo è l’approccio suggerito in questa pagina di esempio di Thrust . Include una variante usando make_discard_iterator .

APPROCCIO # 2 – transform

Questo è l’approccio suggerito da Robert Crovella su CUDA Thrust: reduce_by_key solo su alcuni valori in un array, basandosi su valori in una matrice “chiave” .

APPROCCIO # 3 – inclusive_scan_by_key

Questo è l’approccio suggerito da Eric su Come normalizzare le colonne di matrice in CUDA con prestazioni massime? .

APPROCCIO # 4 – cublasgemv

Usa cuBLAS gemv per moltiplicare la matrice pertinente per una colonna di 1 ‘s.

IL CODICE COMPLETO

Ecco il codice che condensa i due approcci. I file Utilities.cu e Utilities.cuh sono mantenuti qui e omessi qui. TimingGPU.cu e TimingGPU.cuh vengono mantenuti qui e vengono omessi.

 #include  #include  #include  #include  #include  #include  #include  #include  #include  #include  #include "Utilities.cuh" #include "TimingGPU.cuh" // --- Required for approach #2 __device__ float *vals; /**************************************************************/ /* CONVERT LINEAR INDEX TO ROW INDEX - NEEDED FOR APPROACH #1 */ /**************************************************************/ template  struct linear_index_to_row_index : public thrust::unary_function { T Ncols; // --- Number of columns __host__ __device__ linear_index_to_row_index(T Ncols) : Ncols(Ncols) {} __host__ __device__ T operator()(T i) { return i / Ncols; } }; /******************************************/ /* ROW_REDUCTION - NEEDED FOR APPROACH #2 */ /******************************************/ struct row_reduction { const int Ncols; // --- Number of columns row_reduction(int _Ncols) : Ncols(_Ncols) {} __device__ float operator()(float& x, int& y ) { float temp = 0.f; for (int i = 0; i struct MulC: public thrust::unary_function { TC; __host__ __device__ MulC(T c) : C(c) { } __host__ __device__ T operator()(T x) { return x * C; } }; /********/ /* MAIN */ /********/ int main() { const int Nrows = 5; // --- Number of rows const int Ncols = 8; // --- Number of columns // --- Random uniform integer distribution between 10 and 99 thrust::default_random_engine rng; thrust::uniform_int_distribution dist(10, 99); // --- Matrix allocation and initialization thrust::device_vector d_matrix(Nrows * Ncols); for (size_t i = 0; i < d_matrix.size(); i++) d_matrix[i] = (float)dist(rng); TimingGPU timerGPU; /***************/ /* APPROACH #1 */ /***************/ timerGPU.StartCounter(); // --- Allocate space for row sums and indices thrust::device_vector d_row_sums(Nrows); thrust::device_vector d_row_indices(Nrows); // --- Compute row sums by summing values with equal row indices //thrust::reduce_by_key(thrust::make_transform_iterator(thrust::counting_iterator(0), linear_index_to_row_index(Ncols)), // thrust::make_transform_iterator(thrust::counting_iterator(0), linear_index_to_row_index(Ncols)) + (Nrows*Ncols), // d_matrix.begin(), // d_row_indices.begin(), // d_row_sums.begin(), // thrust::equal_to(), // thrust::plus()); thrust::reduce_by_key( thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index(Ncols)), thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index(Ncols)) + (Nrows*Ncols), d_matrix.begin(), thrust::make_discard_iterator(), d_row_sums.begin()); printf("Timing for approach #1 = %f\n", timerGPU.GetCounter()); // --- Print result for(int i = 0; i < Nrows; i++) { std::cout << "[ "; for(int j = 0; j < Ncols; j++) std::cout << d_matrix[i * Ncols + j] << " "; std::cout << "] = " << d_row_sums[i] << "\n"; } /***************/ /* APPROACH #2 */ /***************/ timerGPU.StartCounter(); thrust::device_vector d_row_sums_2(Nrows, 0); float *s_vals = thrust::raw_pointer_cast(&d_matrix[0]); gpuErrchk(cudaMemcpyToSymbol(vals, &s_vals, sizeof(float *))); thrust::transform(d_row_sums_2.begin(), d_row_sums_2.end(), thrust::counting_iterator(0), d_row_sums_2.begin(), row_reduction(Ncols)); printf("Timing for approach #2 = %f\n", timerGPU.GetCounter()); for(int i = 0; i < Nrows; i++) { std::cout << "[ "; for(int j = 0; j < Ncols; j++) std::cout << d_matrix[i * Ncols + j] << " "; std::cout << "] = " << d_row_sums_2[i] << "\n"; } /***************/ /* APPROACH #3 */ /***************/ timerGPU.StartCounter(); thrust::device_vector d_row_sums_3(Nrows, 0); thrust::device_vector d_temp(Nrows * Ncols); thrust::inclusive_scan_by_key( thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index(Ncols)), thrust::make_transform_iterator(thrust::make_counting_iterator(0), linear_index_to_row_index(Ncols)) + (Nrows*Ncols), d_matrix.begin(), d_temp.begin()); thrust::copy( thrust::make_permutation_iterator( d_temp.begin() + Ncols - 1, thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC(Ncols))), thrust::make_permutation_iterator( d_temp.begin() + Ncols - 1, thrust::make_transform_iterator(thrust::make_counting_iterator(0), MulC(Ncols))) + Nrows, d_row_sums_3.begin()); printf("Timing for approach #3 = %f\n", timerGPU.GetCounter()); for(int i = 0; i < Nrows; i++) { std::cout << "[ "; for(int j = 0; j < Ncols; j++) std::cout << d_matrix[i * Ncols + j] << " "; std::cout << "] = " << d_row_sums_3[i] << "\n"; } /***************/ /* APPROACH #4 */ /***************/ cublasHandle_t handle; timerGPU.StartCounter(); cublasSafeCall(cublasCreate(&handle)); thrust::device_vector d_row_sums_4(Nrows); thrust::device_vector d_ones(Ncols, 1.f); float alpha = 1.f; float beta = 0.f; cublasSafeCall(cublasSgemv(handle, CUBLAS_OP_T, Ncols, Nrows, &alpha, thrust::raw_pointer_cast(d_matrix.data()), Ncols, thrust::raw_pointer_cast(d_ones.data()), 1, &beta, thrust::raw_pointer_cast(d_row_sums_4.data()), 1)); printf("Timing for approach #4 = %f\n", timerGPU.GetCounter()); for(int i = 0; i < Nrows; i++) { std::cout << "[ "; for(int j = 0; j < Ncols; j++) std::cout << d_matrix[i * Ncols + j] << " "; std::cout << "] = " << d_row_sums_4[i] << "\n"; } return 0; } 

RISULTATI DEL TEMPO (testati su un Kepler K20c)

 Matrix size #1 #1-v2 #2 #3 #4 #4 (no plan) 100 x 100 0.63 1.00 0.10 0.18 139.4 0.098 1000 x 1000 1.25 1.12 3.25 1.04 101.3 0.12 5000 x 5000 8.38 15.3 16.05 13.8 111.3 1.14 100 x 5000 1.25 1.52 2.92 1.75 101.2 0.40 5000 x 100 1.35 1.99 0.37 1.74 139.2 0.14 

Sembra che l'approccio n. 1 e n. 3 superi l'approccio n. 2, tranne nei casi di piccoli numeri di colonne. L'approccio migliore, tuttavia, è l'approccio n. 4, che è significativamente più conveniente rispetto agli altri, a condizione che il tempo necessario per creare il piano possa essere ammortizzato durante il calcolo.