Comparaison de vecteur long

**Gouyon** · 25/10/2022, 09h11

Bonjour

Je suis en train de faire des comparaisons deux à deux de vecteurs très long (plusieurs millions de points). Il s'agit de savoir si ces deux vecteurs ont des points en commun et où sont ces points en commun dans chaque vecteur.
Les vecteurs sont des ensemble de coordonnées géographiques (longitude, latitude).
Avant de passer au vecteur réel j'ai fait un petit test sur des vecteurs plus petit mais ça ne fonctionne pas car seul les premier élément sont testé et je n'arriva pas bien à savoir pourquoi.
Voici mon code

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
 
 
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
 
#include <stdio.h>
 
cudaError_t compVect(unsigned char *c, const int *lona, const int *lata, const int *lonb, const int *latb, unsigned int sizeA, unsigned int sizeB);
 
__global__ void compKernel(unsigned char *c, const int *lona, const int *lata, const int *lonb, const int *latb, int nbA, int nbB)
{
	int i = threadIdx.x + blockIdx.x*blockDim.x;
	int j = threadIdx.y + blockIdx.y*blockDim.y;
	int k = i + j * nbA;		
	while (k < nbA*nbB)
	{
		if (lona[i] == lonb[j] && lata[i] == latb[j])
			c[k] = 1;
		else
			c[k] = 0;
		i += gridDim.x*blockDim.x;
		j += gridDim.y*blockDim.y;
		k = i + j * nbA;
	}
}
 
int main()
{
	int *lon[2];
	int *lat[2];
	int nb[2];
#ifdef TEST
	lon[0] = new int[15];
	lat[0] = new int[15];
	lon[1] = new int[7];
	lat[1] = new int[7];
	nb[0] = 15;
	nb[1] = 7;
	for (int i = 0; i < 15; i++)
	{
		lon[0][i] = i + 1;
		lat[0][i] = 2 * i;
 
	}
 
	for (int i = 0; i < 7; i++)
	{
		lon[1][i] = 3 * i + 1;
		lat[1][i] = 4 * i;
 
	}
 
	lon[0][14] = lon[1][6];
	lat[0][14] = lat[1][6];
	for (int i = 0; i < 15; i++)
	{
		printf("L=%i l=%i   ", lon[0][i], lat[0][i]);
	}
	printf("\n");
	for (int i = 0; i < 7; i++)
	{
		printf("L=%i l=%i   ", lon[1][i], lat[1][i]);
	}
	unsigned char *c = new unsigned char[nb[0] * nb[1]];
	// Add vectors in parallel.
	cudaError_t cudaStatus = compVect(c, lon[0], lat[0], lon[1], lat[1], nb[0], nb[1]);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "addWithCuda failed!");
		return 1;
	}
	for (int i = 0; i < nb[0] * nb[1]; i++)
	{
		if (c[i] == 1)
		{
			int y = i / nb[0];
			int x = i % nb[0];
			printf("A=%i B=%i", x, y);
		}
	}
 
 
	// cudaDeviceReset must be called before exiting in order for profiling and
	// tracing tools such as Nsight and Visual Profiler to show complete traces.
	cudaStatus = cudaDeviceReset();
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaDeviceReset failed!");
		return 1;
	}
 
	return 0;
}
 
// Helper function for using CUDA to add vectors in parallel.
cudaError_t compVect(unsigned char *c, const int *lona, const int *lata, const int *lonb, const int *latb, unsigned int sizeA, unsigned int sizeB)
{
	int *dev_lona = 0;
	int *dev_lonb = 0;
	int *dev_lata = 0;
	int *dev_latb = 0;
 
	unsigned char *dev_c = 0;
	cudaError_t cudaStatus;
 
	// Choose which GPU to run on, change this on a multi-GPU system.
	cudaStatus = cudaSetDevice(0);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");
		goto Error;
	}
 
	// Allocate GPU buffers for three vectors (two input, one output)    .
	cudaStatus = cudaMalloc((void**)&dev_c, sizeA*sizeB * sizeof(unsigned char));
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMalloc failed!");
		goto Error;
	}
 
	cudaStatus = cudaMalloc((void**)&dev_lona, sizeA * sizeof(int));
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMalloc failed!");
		goto Error;
	}
	cudaStatus = cudaMalloc((void**)&dev_lata, sizeA * sizeof(int));
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMalloc failed!");
		goto Error;
	}
	cudaStatus = cudaMalloc((void**)&dev_lonb, sizeB * sizeof(int));
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMalloc failed!");
		goto Error;
	}
	cudaStatus = cudaMalloc((void**)&dev_latb, sizeB * sizeof(int));
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMalloc failed!");
		goto Error;
	}
	// Copy input vectors from host memory to GPU buffers.
	cudaStatus = cudaMemcpy(dev_lona, lona, sizeA * sizeof(int), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMemcpy failed!");
		goto Error;
	}
	cudaStatus = cudaMemcpy(dev_lata, lata, sizeA * sizeof(int), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMemcpy failed!");
		goto Error;
	}
	cudaStatus = cudaMemcpy(dev_lonb, lonb, sizeB * sizeof(int), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMemcpy failed!");
		goto Error;
	}
	cudaStatus = cudaMemcpy(dev_latb, latb, sizeB * sizeof(int), cudaMemcpyHostToDevice);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMemcpy failed!");
		goto Error;
	}
	// Launch a kernel on the GPU with one thread for each element.
	dim3 grid(2, 2);
	compKernel << <grid, grid >> > (dev_c, dev_lona, dev_lata, dev_lonb, dev_latb,sizeA,sizeB);
 
	// Check for any errors launching the kernel
	cudaStatus = cudaGetLastError();
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "addKernel launch failed: %s\n", cudaGetErrorString(cudaStatus));
		goto Error;
	}
 
	// cudaDeviceSynchronize waits for the kernel to finish, and returns
	// any errors encountered during the launch.
	cudaStatus = cudaDeviceSynchronize();
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching addKernel!\n", cudaStatus);
		goto Error;
	}
 
	// Copy output vector from GPU buffer to host memory.
	cudaStatus = cudaMemcpy(c, dev_c, sizeA * sizeB * sizeof(unsigned char), cudaMemcpyDeviceToHost);
	if (cudaStatus != cudaSuccess) {
		fprintf(stderr, "cudaMemcpy failed!");
		goto Error;
	}
 
Error:
	cudaFree(dev_c);
	cudaFree(dev_lona);
	cudaFree(dev_lonb);
	cudaFree(dev_lata);
	cudaFree(dev_latb);
	return cudaStatus;
}

Je pense que le problème se situe dans cette boucle mais je n'arrive pas à comprendre pourquoi

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 
__global__ void compKernel(unsigned char *c, const int *lona, const int *lata, const int *lonb, const int *latb, int nbA, int nbB)
{
	int i = threadIdx.x + blockIdx.x*blockDim.x;
	int j = threadIdx.y + blockIdx.y*blockDim.y;
	int k = i + j * nbA;		
	while (k < nbA*nbB)
	{
		if (lona[i] == lonb[j] && lata[i] == latb[j])
			c[k] = 1;
		else
			c[k] = 0;
		i += gridDim.x*blockDim.x;
		j += gridDim.y*blockDim.y;
		k = i + j * nbA;
	}
}

Toute aide est la bienvenue

**Marlan** · 26/10/2022, 06h42

Bonjour,
Je vois une erreur au niveau de la définition des éléments de la grille et du bloc : il est écrit que vous choisissez une grille dimensionnée en 2x2.
Sur le calcul GPU, la grille correspond à une grille de blocs de threads. La grille comme les blocs pouvant être 1, 2 ou 3D (ce qui donne une description 6D potentielle de votre problème).
Puisque vous traitez ici le cas de vecteurs, il est logique de donner une description 1D du problème (donc une grille 1D de blocs 1D de threads).
Pour ce faire, vous pouvez déclarer :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 
dim3 block,grid;
block.x=256; //en général, ceci est la valeur recommandée par Nvidia
grid.x=(N+block.x-1)/block.x;//le problème de taille N est découpé en grid.x blocs de block.x threads.

Dans votre noyau, l'indice du thread correspondant devient alors

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
const int tid=threadIdx.x+blockIdx.x*blockDim.x;

Dans votre cas, l'idée est maintenant de décrire le problème en prenant le plus grand des deux vecteurs (afin de maximiser le nombre de threads utilisés) puis de réaliser vos comparaisons en chargeant le second vecteur. Notez que vous pouvez ici utiliser la mémoire partagée pour partager de l'information au sein d'un même bloc de threads (ce qui vous évitera que tous les threads chargent le même bloc de donnée).

Bonne journée,

Marlan

**Gouyon** · 26/10/2022, 22h04

En fait je suis parti sur un tableau à 2 dimensions la première est celle du vecteur A et la seconde du vecteur B. Car je dois comparer le premier point du vecteur A avec tous les points du vecteur B, puis le second point du vecteur A avec tous les points du vecteur B etc etc

Autrement dit l'opération que je fais est la suivante

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
 
for(int i=0;i<nbA;i++)
{
 for(int j=0;j<nbB;j++)
   {
      k=j+i*nbA;
      if (vecA[i]==vecB[j])
       c[k]=1;
     else
       c[k]=0;
  }
}

Si je met

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
dim3 grid(4, 4);

Ça fonctionne bien car j'ai 4 bloc de 4 4 thread pour chaque dimension et ça couvre mon petit exemple. Mais ce n'est pas viable car les données réelles son beaucoup plus grande et dépassent les limites de la carte.

En fait je me suis inspiré de cet exemple en 1D qui est une addition de 2 vecteurs

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
 
__global__ void add(int *a,int*b,int*c)
{
  tid=threadIdx.x+blockIdX.x+blockDim.x;
  while(tid<N)
  {
   c[tid]=a[tid]+b[tid];
  tid+=blockDim.x*gridDim.x;
  }
}
 
//appel à add
add<<<128,128>>>(dev_a,dev_b,dev_c);

Mais manifestement je n'ai pas tout compris

**Marlan** · 07/11/2022, 05h48

Bonjour,
L'exemple que vous avez choisi pour l'addition est étrange à mon sens car l'idée serait plutôt de réaliser l'addition en affectant un thread à chaque coefficient.
Personnellement, je ferai quelque chose comme :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
 
void
__global__ add_vectors(const int N,const double * __restrict__ a, const double * __restrict__ b,double * __restrict__ c)
{
  int tid=threadIdx.x+blockDim.x*blockIdx.x;
  if(tid<N)
  {
    c[tid]=a[tid]+b[tid];
  }
}

Puis un appel par :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
dim3 grid,block;
block.x=256;//Valeur suggérée par NVidia pour ce type de problème (mais devrait très bien marche avec 128 aussi)
grid.x=(N+block.x-1)/block.x;
add_vectors<<<grid,block,0,0>>>(N,a,b,c);//réalise c=a+b, N est la taille du vecteur

En utilisant cette approche, la boucle while disparait au profit de l'emploi de threads directement.
Cet algorithme n'est alors limité que par les capacités mémoire RAM de la carte (suivant la génération, plusieurs Go).
Si vous deviez traiter un vecteur dépassant ces capacités, il faudrait alors changer d'approche en découpant votre vecteur en plusieurs parties puis en transférant une partie du cpu vers le gpu, réaliser l'opération et transférer le résultat puis passer à une autre partie du vecteur et ainsi de suite.
Le noyau appelé en revanche restera le même quelle que soit la portion de vecteur considérée.
En d'autres termes, cela donnerait qqch comme :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
//Calcul du nombre de blocs à réaliser pour les très grands vecteurs
const int sub_vecotr_size=67108864;//Sous vecteur de taille 67M éléments
const int nb_blocks=N/sub_vecotr_size;//Nb de sous vecteurs
const int rem=N-nb_blocks*sub_vecotr_size;//Calcul du reste
for(int bid=0;bid<nb_blocks;++bid)
{
  //Transfert CPU->GPU de la portion de vecteur à traiter sur GPU
  //Calcul C=A+B avec comme taille de vecteur sub_vector_size
  //Transfert du résultat GPU->CPU
}
//Traitement du dernier bloc correspondant au reste
if(rem>0)
{
  ...
}

Vis à vis de votre problème de base, vous risquez par contre d'avoir du mal à passer à l'échelle sur votre vrai jeu de données : la taille na*nb risque d'être rédhibitoire en terme de consommation mémoire.

Dans le cas contraire, l'algorithme sur lequel je partirai serait :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
void 
__global__ kernel(const int na, const int nb, const * __restrict__ a, const * __restrict__ b, char *__restrict__ c)
{
  const int tid=threadIdx.x+blockIdx.x*blockDim.x;
  if(tid<na)
  {
    const int vala=a[tid];
    for(int idx=0;idx<nb;++idx)
    {
      if(vala==b[idx])
      {
        c[tid+idx*nb]=1;
      }
      else
      {
        c[tidx+idx*nb]=0;
      }
    }
  }
}

avec c un vecteur dimensionné à na*nb (mais de nouveau, cette taille peut être prohibitive).
Bonne journée,

Comparaison de vecteur long

CUDA

Discussions similaires

Partager

Partager