Implémentation De La Programmation Parallèle Avec CUDA

Tim inclus dans Technologie-Informatique Science-Des-Données-Apprentissage-Machine

2023-05-08 2024-12-25 environ 1500 mots 7 minutes - vues

Contenus

Cet article utilise des exemples comme la transposition de matrices pour montrer comment accélérer la résolution de problèmes à grande échelle en utilisant les nombreux cœurs du GPU via l’écriture de fonctions noyau CUDA.

0.1 Transposition de matrice

0.1.1 Processus de l’algorithme

Stocker la matrice à transposer dans la mémoire GPU.
Allouer de l’espace sur le GPU pour stocker la matrice transposée.
Définir une fonction noyau CUDA pour réaliser la transposition de la matrice. Cette fonction noyau doit utiliser les concepts de blocs de threads et de grilles de threads pour traiter tous les éléments de la matrice. Dans chaque bloc de threads, les threads peuvent utiliser la mémoire partagée pour traiter les données. Enfin, utiliser la mémoire globale pour écrire le résultat sur le GPU.
Appeler la fonction noyau CUDA pour exécuter la transposition de la matrice.
Copier la matrice transposée de la mémoire GPU vers la mémoire de l’hôte.
Libérer la mémoire GPU.

0.1.2 Implémentation du code

Sur la base du code traditionnel, utiliser la mémoire partagée pour optimiser l’accès à la mémoire globale, tout en utilisant l’opération de padding pour éviter les conflits de banque.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51


const int N = 10000;
const int TILE_DIM = 32;
const int grid_size_x = (N + TILE_DIM - 1) / TILE_DIM;
const int grid_size_y = grid_size_x;
const dim3 block_size(TILE_DIM, TILE_DIM);
const dim3 grid_size(grid_size_x, grid_size_y);
__global__ void transpose(const float *A, float *B, const int N){  
    __shared__ float S[TILE_DIM][TILE_DIM+1];
    int bx = blockIdx.x * TILE_DIM;
    int by = blockIdx.y * TILE_DIM;
    // Lecture séquentielle
    int nx1 = bx + threadIdx.x;
    int ny1 = by + threadIdx.y;
    if (nx1 < N && ny1 < N)
        S[threadIdx.y][threadIdx.x] = A[ny1 * N + nx1];
    __syncthreads();
    // Écriture séquentielle
    int nx2 = bx + threadIdx.y;
    int ny2 = by + threadIdx.x;
    if (nx2 < N && ny2 < N)
        B[nx2 * N + ny2] = S[threadIdx.x][threadIdx.y];
}

int main(){  
    int nn = N * N;
    int nBytes = nn * sizeof(float);
    float a = 1.1;
    float b = 2.2;
    // Allocation de l'espace CPU
    float* A = (float*)malloc(nBytes);
    float* B = (float*)malloc(nBytes);
    // Initialisation de la matrice
    for (int i = 1; i <= nn; i++) {
        if (i % 2 == 0)
            A[i-1] = a;
        else
            A[i-1] = b;
    }
    // Allocation de l'espace GPU (à définir d'abord)
    float* d_A, * d_B;
    cudaMalloc((void**)&d_A, nBytes);
    cudaMalloc((void**)&d_B, nBytes);
    cudaMemcpy(d_A, A, nBytes, cudaMemcpyHostToDevice);
    transpose<<<grid_size, block_size >>> (d_A, d_B, N);
    cudaMemcpy(B, d_B, nBytes, cudaMemcpyDeviceToHost);
    free(A);
    free(B);
    cudaFree(d_A);
    cudaFree(d_B);
    return 0;
}

0.2 Multiplication de matrices

0.2.1 Conception de la fonction noyau

Utilisation de la mémoire partagée : Il faut extraire une partie des deux matrices d’origine, la considérer comme une tuile, et la placer dans la mémoire partagée pour réduire l’accès à la mémoire globale, améliorant ainsi les performances du programme. Traitement des matrices non carrées : La multiplication de deux matrices non carrées nécessite d’ajuster la taille des tuiles dans la mémoire partagée et d’ajouter des instructions conditionnelles pour éviter les erreurs de dépassement de mémoire.

0.2.2 Implémentation du code

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29


__global__ void matrixmul(const float *A, const float *B, float *C,const int N,const int M,const int K){   //N*M x M*K
    __shared__ float A_S[TILE_DIM][TILE_DIM];
    __shared__ float B_S[TILE_DIM][TILE_DIM];
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int bx = blockIdx.x;
    int by = blockIdx.y;
    int row = by * TILE_DIM + ty;
    int col = bx * TILE_DIM + tx;
    float value = 0;
    for (int ph = 0; ph < M / TILE_DIM+1; ph++) {
        if (row < N && ph * TILE_DIM + tx < M)
            A_S[ty][tx] = A[row * M + ph * TILE_DIM + tx];
        else
            A_S[ty][tx] = 0.0;
        if (col < K && ph * TILE_DIM + ty < M) 
            B_S[ty][tx] = B[(ph * TILE_DIM + ty) * K + col];
        else
            B_S[ty][tx] = 0.0;
        __syncthreads();

        for (int k = 0; k < TILE_DIM; k++)
            value += A_S[ty][k] * B_S[k][tx];
        __syncthreads();
    }

    if (row < N && col < K)
        C[row*K+col]=value;
}

0.3 Histogramme statistique

Utilisation de la mémoire partagée : Créer un tableau temporaire dans la mémoire partagée pour réduire l’accès à la mémoire globale, améliorant ainsi les performances du programme. Boucle sur la grille : Pour de grandes quantités de données, un seul thread doit traiter plusieurs données. Opérations atomiques : Que ce soit pour les résultats en mémoire partagée ou en mémoire globale, il est nécessaire d’utiliser des opérations atomiques pour éviter les conflits.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


__global__ void histo_kernel(unsigned char *buffer, long size, unsigned int *histo){
	__shared__ unsigned int temp[256];
	temp[threadIdx.x] = 0;
	__syncthreads();
	int i = threadIdx.x + blockIdx.x * blockDim.x;
	int offset = blockDim.x *gridDim.x;
	while (i<size){
		atomicAdd(&temp[buffer[i]], 1);
		i += offset;
	}
	__syncthreads();
	atomicAdd(&(histo[threadIdx.x]), temp[threadIdx.x]);
}

0.4 Réduction et somme

Boucle sur la grille : Pour des quantités massives de données (>100w), même avec de nombreux cœurs GPU, il est impossible d’avoir autant de threads, donc un seul thread doit traiter plusieurs données. Optimisation de la mémoire partagée : Stocker le résultat de la somme de plusieurs éléments par chaque thread dans la mémoire partagée peut réduire l’accès à la mémoire globale. Méthode de couplage entrelacé : Utiliser la méthode de couplage entrelacé pour la réduction peut atténuer une partie de la divergence des faisceaux de threads. Opérations atomiques : L’accumulation des résultats dans la mémoire globale peut entraîner des conflits, nécessitant l’utilisation d’opérations atomiques.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


__global__ void _sum_gpu(int *input, int count, int *output)
{
    __shared__ int sum_per_block[BLOCK_SIZE];

    int temp = 0;
    for (int idx = threadIdx.x + blockDim.x * blockIdx.x;
         idx < count; idx += gridDim.x * blockDim.x
	)
        temp += input[idx];
    sum_per_block[threadIdx.x] = temp;
    __syncthreads();

    //**********étape de sommation en mémoire partagée***********
    for (int length = BLOCK_SIZE / 2; length >= 1; length /= 2)
    {
        int double_kill = -1;
	if (threadIdx.x < length)
	    double_kill = sum_per_block[threadIdx.x] + sum_per_block[threadIdx.x + length];
	__syncthreads();  //pourquoi avons-nous besoin de deux __syncthreads() ici, et,
	
	if (threadIdx.x < length)
	    sum_per_block[threadIdx.x] = double_kill;
	__syncthreads(); 
    } //la somme partielle par bloc est sum_per_block[0]

    if (threadIdx.x == 0) atomicAdd(output, sum_per_block[0]);

0.5 Problème TOP K

Le problème TOP K consiste à trouver les K plus grands ou plus petits éléments dans un ensemble de données. L’utilisation de la réduction CUDA peut résoudre efficacement le problème TOP K. Voici les étapes détaillées pour réaliser le tri et sélectionner les K plus grands/petits éléments en utilisant la réduction CUDA :

Définir un type de tuple binaire, contenant la valeur et l’indice, pour stocker les données d’origine et leur indice (afin de restaurer les données d’origine après le tri)
Parcourir les données d’origine et les stocker dans le tuple binaire
Effectuer une opération de réduction sur le tuple binaire pour obtenir les indices des K plus grands/petits éléments
Restaurer les données d’origine et trier selon les indices pour obtenir les K plus grands/petits éléments

Processus de mise en œuvre spécifique :

Copier les données dans la mémoire GPU float *d_data; cudaMemcpy(d_data, h_data, size, cudaMemcpyHostToDevice);
Stocker les données dans le tuple binaire

1
2
3
4
5


typedef struct {     float value;     int index; } Tuple;  
Tuple *d_tuples; 
int threadsPerBlock = 256; 
int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
initializeTuples<<<blocksPerGrid, threadsPerBlock>>>(d_data, d_tuples, n);

Effectuer une opération de réduction sur le tuple binaire pour obtenir les indices des K plus grands/petits éléments

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17


int *d_indices;
kReduceKernel<<<blocksPerGrid, threadsPerBlock>>>(d_tuples, d_indices, n, k);
__global__ void kReduceKernel(Tuple *input, int *output, int n, int k) {
    extern __shared__ Tuple shared[];
    int tid = threadIdx.x;
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    shared[tid] = (i < n) ? input[i] : Tuple{0, 0};
    __syncthreads();
    for (int s = blockDim.x / 2; s > 0; s >>= 1) {
        if (tid < s)
            shared[tid] = (shared[tid].value > shared[tid + s].value) ? shared[tid] : shared[tid + s];
        __syncthreads();
    }

    if (tid == 0)
        output[blockIdx.x] = shared[0].index;
}

Restaurer les données d’origine sur le CPU et trier selon les indices pour obtenir les K plus grands/petits éléments

1
2
3
4
5


cudaMemcpy(h_indices, d_indices, size, cudaMemcpyDeviceToHost);  
for (int i = 0; i < k; ++i) {     
    int index = h_indices[i];     
    h_result[i] = h_data[index]; }  
std::sort(h_result, h_result + k);

Buy me a coffee~

Faire un don

Alipay

PayPal

WeChat Pay