1 pièce(s) jointe(s)
performance : tableau a plusieurs dimensions
Salut,
J'ai une quantité bi-dimensionnelle : f(i,j) que je veux accéder en lecture et en écriture, rapidement. Le domaine sur lequel f est défini contient beaucoup d'éléments (possiblement plusieurs millions de 'double').
A chaque élément f(i,j), je vais avoir besoin d'accéder rapidement (en lecture surtout) aux élements "autour", c'est à dire :
Code:
1
2
3
4
5
|
f(i+1,j)
f(i-1,j)
f(i,j+1)
f(i,j-1) |
Je me des questions quant à l'implémentation la plus efficace de ces données du point de vue rapidité.
J'ai trois solutions (peut-on en imaginer d'autres ?)
La première consiste a faire un "vrai" tableau 2D dynamique :
Code:
1
2
3
4
5
6
7
8
9
|
typedef unsigned int uint;
typedef struct array2D
{
double **a;
uint nx;
uint ny;
} Array2D; |
que je construit comme ceci :
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
30
31
32
33
34
35
|
Array2D *array2D_new(uint nx, uint ny)
{
Array2D *this;
uint i;
uint k;
this = calloc(1, sizeof *this);
if (this == NULL) return NULL;
this->a = calloc(nx, sizeof *this->a);
if (this->a == NULL)
{
free(this);
return NULL;
}
for (i=0; i < nx; i ++)
{
this->a[i] = calloc(ny, sizeof *this->a[i]);
if (this->a[i] == NULL)
{
for (k=0; k <= i; k ++)
free (this->a[k]);
free(this);
return NULL;
}
}
this->nx = nx;
this->ny = ny;
return this;
} |
détruit ainsi :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
|
void
array2D_delete(Array2D *this)
{
uint i;
if (this != NULL)
{
for (i=0; i < this->nx; i++)
{
if (this->a[i]) free(this->a[i]);
}
}
free(this);
} |
et pour lequel j'accède aux éléments de cette manière :
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
|
void array2D_stencil(const Array2D *grid)
{
uint i,j;
double tmp;
uint nx, ny;
nx = grid->nx;
ny = grid->ny;
for (i=1; i < nx-1; i ++)
{
for (j = 1; j < ny-1; j++)
{
// do some stuff with the four
// neighbouring points
tmp = 0.25*(grid->a[i-1][j ] +
grid->a[i ][j+1] +
grid->a[i+1][j ] +
grid->a[i ][j-1]);
}
}
} |
La seconde méthode, consiste à faire un tableau dynamique 1D, et que l'utilisateur calcule lui même l'indice 1D à partir de deux boucles imbriquées sur i et j.
Ca se déclare comme ceci :
Code:
1
2
3
4
5
6
7
|
typedef struct _array2D1D
{
double *a;
uint nx;
uint ny;
} Array2D1D; |
se construit ainsi :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
|
Array2D1D *array2D1D_new(uint nx, uint ny)
{
Array2D1D *this;
this = calloc(1, sizeof *this);
if (!this) return NULL;
this->a = calloc(nx*ny, sizeof *this->a);
if (! this->a)
{
free(this);
return NULL;
}
this->nx = nx;
this->ny = ny;
return this;
} |
se détruit ainsi
Code:
1
2
3
4
5
6
7
8
9
10
11
|
void
array2D1D_delete(Array2D1D *this)
{
if (this)
{
if (this->a) free (this->a);
free(this);
}
} |
et on accède aux éléments ainsi :
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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
|
void array2D1D_stencil_a(const Array2D1D *grid)
{
uint i,j;
uint ij, ij1, ij2, ij3, ij4;
double tmp;
double *a;
uint nx, ny;
nx = grid->nx;
ny = grid->ny;
a = grid->a;
// +ij2
// |
// |
// ij1 | ij3
// +--------+ij------+
// |
// |
// |
// +
// ij4
for (i=1; i < nx-1; i ++)
{
for (j = 1; j < ny-1; j++)
{
ij = (i ) + (j )*(nx);
ij1 = (i-1) + (j )*(nx);
ij2 = (i ) + (j+1)*(nx);
ij3 = (i+1) + (j )*(nx);
ij4 = (i ) + (j-1)*(nx);
// do some stuff with the four
// neighbouring points
tmp = 0.25*(a[ij1] + a[ij2] + a[ij3] + a[ij4]);
}
}
} |
J'ai chronométré les fonctions *_stencil() pour les deux objets, sur des tableaux de 8192*8192 éléments, et ça 10fois de suite, afin de m'assurer que le chronométrage ne dépend pas d'autres processus tournant sur ma machine.
j'ai été très surpris de constater que la méthode 1 est BIEN plus rapide. En effet, voici les temps, en seconde des deux méthodes :
Code:
1
2
3
4
5
6
7
8
9
10
11
12
|
méthode 1 méthode 2
1.10430538 8.25240617
0.967257762 7.72914395
0.963137127 7.71300082
0.963082132 7.70107212
0.963722737 7.69799643
0.969022776 7.80377674
0.976541885 7.83630202
0.963828863 7.70916793
0.964121218 7.71000535
0.963386301 7.70990839 |
En fait, j'aurai dit que la méthode 1 serait la plus lente. Car la "seconde dimension" du tableau est allouée à un endroit différent dans la mémoire que la première dimension... et chaque ligne de la seconde dimension est allouée à un endroit différent. En allant chercher des valeur aux points (i+1,j), (i-1,j), (i,j+1), (i,j-1) j'aurais donc pensé qu'on allait charger des zones différentes en cache, et perdre en rapidité.
La méthode 2 en revanche, présente en effet un surcout car elle demande, pour charger un élément à la ligne supérieure (j+1), de charger toutes les données i de la ligne j. Je pensais que ça ne serait pas efficace, mais tout de même plus que la méthode 1, car même si on doit sauter beaucoup de point, les données restent contigues en mémoire.
Devant le résultat, j'ai pensé que le surcout venait peut-être du calcul d'indice que je fais à chaque fois :
Code:
1
2
3
4
5
|
ij1 = (i-1) + (j )*(nx);
ij2 = (i ) + (j+1)*(nx);
ij3 = (i+1) + (j )*(nx);
ij4 = (i ) + (j-1)*(nx); |
J'ai donc légèrement modifié mon objet, de sorte à ce que les indices des points autour de (i,j) de soient plus calculés à chaque fois, mais calculés une bonne fois pour toute lors de la création de l'objet, pour chaque point (i,j) :
Code:
1
2
3
4
5
6
7
8
|
typedef struct _array2D1D
{
double *a;
uint nx;
uint ny;
uint *stencil[4];
} Array2D1D; |
le "stencil" est donc un ensemble de 4 indices par point de tableau 'a', il est créé dans le constructeur précédent, un peu modifié pour la cause :
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
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
|
Array2D1D *array2D1D_new(uint nx, uint ny)
{
Array2D1D *this;
uint i,j,ij;
uint ij1, ij2, ij3, ij4;
this = calloc(1, sizeof *this);
if (!this) return NULL;
this->a = calloc(nx*ny, sizeof *this->a);
if (! this->a)
{
free(this);
return NULL;
}
this->nx = nx;
this->ny = ny;
// allocation pas verifiee.. a l'arrache
this->stencil[0] = calloc(nx*ny, sizeof *this->stencil[0]);
this->stencil[1] = calloc(nx*ny, sizeof *this->stencil[1]);
this->stencil[2] = calloc(nx*ny, sizeof *this->stencil[2]);
this->stencil[3] = calloc(nx*ny, sizeof *this->stencil[3]);
// +ij2
// |
// |
// ij1 | ij3
// +-----+ij------+
// |
// |
// |
// +
// ij4
for (i=0; i < nx; i++)
{
for (j=0; j < ny; j++)
{
ij = (i ) + (j )*(nx);
ij1 = (i-1) + (j )*(nx);
ij2 = (i ) + (j+1)*(nx);
ij3 = (i+1) + (j )*(nx);
ij4 = (i ) + (j-1)*(nx);
this->stencil[0][ij] = ij1;
this->stencil[1][ij] = ij2;
this->stencil[2][ij] = ij3;
this->stencil[3][ij] = ij4;
}
}
return this;
} |
Cette méthode, je pensais, permet de s'affranchir du cout des calculs d'indices, au prix d'un cout mémoire additionnel que je suis prêt à payer.
l'accès au tableau se fait à présent ainsi :
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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
void array2D1D_stencil_b(const Array2D1D *grid)
{
uint i,j;
double tmp;
double *a;
uint *ij1;
uint *ij2;
uint *ij3;
uint *ij4;
uint nx, ny;
nx = grid->nx;
ny = grid->ny;
a = grid->a;
ij1 = grid->stencil[0];
ij2 = grid->stencil[1];
ij3 = grid->stencil[2];
ij4 = grid->stencil[3];
// +ij2
// |
// |
// ij1 | ij3
// +--------+ij------+
// |
// |
// |
// +
// ij4
for (i=1; i < nx-1; i ++)
{
for (j = 1; j < ny-1; j++)
{
// do some stuff with the four
// neighbouring points
tmp = 0.25*(a[ij1[i + j*nx]] + a[ij2[i+j*nx]] + a[ij3[i+j*nx]] + a[ij4[i+j*nx]]);
}
}
} |
Je relance mon code de mesure de temps, et à mes deux colonnes précédentes s'ajoute une troisième, mesurant le temps en secondes passé dans la routine array2D1D_stencil_b()
Code:
1
2
3
4
5
6
7
8
9
10
11
12
|
methode 1 méthode 2 méthode 3
1.10430538 8.25240617 26.6270934
0.967257762 7.72914395 26.6097884
0.963137127 7.71300082 26.6851487
0.963082132 7.70107212 26.6141512
0.963722737 7.69799643 26.7747686
0.969022776 7.80377674 26.8523572
0.976541885 7.83630202 26.7336635
0.963828863 7.70916793 26.6687285
0.964121218 7.71000535 26.6810654
0.963386301 7.70990839 26.7819276 |
Alors ça ! c'est encore beaucoup plus lent que la seconde méthode !!
en gros, je pensais que les temps seraient :
méthode3 < méthode2 < méthode 1
et finalement c'est tout l'inverse !
Quelqu'un peut-il m'expliquer ça ? Il semble finalement qu'il faille privilégier la méthode 1, et de loin....
J'ai mis le code dans un zip en pièce jointe. Ca devrait compiler sans probleme avec un simple 'make' sous linux (ça ne tourne pas sous mac, et je ne sais pas pour windows, à cause de la bibliothèque 'lrt').
Merci :)
