Discussion :

[moomba]

Bonjour

Je cherche à optimiser une opération que je répète un très grand nombre de fois, et qui correspond à une convolution (c'est en fait ici une dérivée seconde centrée discrétisée).

D'une manière simple, le code est le suivant :

Tous les tableaux sont des réels double précision.
DF est de taille 512*512 (1 à 512, 1 à 512).
F est de taille 514*514 (0 à 513, 0 à 513).
Le tableau M des coefficients est de taille 512*512*5 (1 à 512, 1 à 512, 1 à 5).

Code :

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

1
2
3
4
5
6
7
for ( j=0 ; j<=512 ; j++ )
{
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = M[i][j][1] * F[i][j] + M[i][j][2] * F[i-1][j] + M[i][j][3] * F[i+1][j] + M[i][j][4] * F[i][j-1] + M[i][j][5] * F[i][j+1] ;
   }
}

Il est possible que je me soit trompé dans l'ordre des boucles (i avant j, ou l'inverse, pour optimiser le cache). En principe, je code en Fortran, donc je ne connais pas les optimisations du C++.

Ici, je cherche plutôt à utiliser la bibliothèque BLAS qui devrait me permettre de gagner un bon facteur 2 ou 3 sur cette opération. Cependant, après avoir lus la doc, je ne sais pas trop quel routine utiliser... Il doit falloir dérouler les boucles pour clarifier la chose.

Quelqu'un as-t-il déjà utilisé BLAS sur ce genre d'opérations ?

En vous remerciant d'avance

Moomba

[Invité]

Pour faire une convolution le plus rapidement possible, le plus simple est de passer dans le domaine fréquentiel. Avec une FFT, la complexité est en (n*log(n))...Par rapport à N^2 de la convolution temporelle, y'a meêm pas à se poser de question.

Donc mieux vaut t'intéresser à une bibliothèque qui fait des FFT rapidement comme fftw.

[moomba]

Merci pour cette réponse.

Je me suis intéressé aux FFT, mais il y a deux problèmes majeurs :

Les calculs de chimie de mon code doivent se faire dans le domaine non fréquentielle, donc je doit constamment passer d'un domaine à l'autre...

Autre problème, je calcul avec une mémoire distribuée (MPI) : ici, la convolution se fait de F vers DF, en utilisant les ghosts de F qui ont été échangé avec les processus voisins. Or la FFT n'est pas distribuable en mémoire...

C'est pour cela que je m'intéresse à BLAS

[Astraya]

Il est possible que je me soit trompé dans l'ordre des boucles (i avant j, ou l'inverse, pour optimiser le cache)

Si ton cache est un image tu doit effectivement faire les lignes j puis les colonnes i,. Il y auras moins de saut en ram et un accès beaucoup plus rapide à la mémoire. Pour le reste tout dépend de l'ordre de lecture des données.

BLAS connais pas

[moomba]

Merci pour cette précision
Comme ça je saurais pour le coup suivant.

Après recherches, il semble que BLAS ne puisse pas répondre à mes attentes. Je cherche donc à optimiser simplement ma boucle.

Les gars d'INTEL m'ont conseillé l'astuce suivante :


Avant :

Code :

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

1
2
3
4
5
6
7
for ( j=0 ; j<=512 ; j++ )
{
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = M[i][j][1] * F[i][j] + M[i][j][2] * F[i-1][j] + M[i][j][3] * F[i+1][j] + M[i][j][4] * F[i][j-1] + M[i][j][5] * F[i][j+1] ;
   }
}

Après :

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
 
for ( j=0 ; j<=512 ; j++ )
{
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = M[i][j][1] * F[i][j];
   }
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = DF + M[i][j][2] * F[i-1][j];
   }
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = DF + M[i][j][3] * F[i+1][j];
   }
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = DF + M[i][j][4] * F[i][j-1];
   }
   for ( i=0 ; i<=512 ; i++ )
   {
      DF[i][j] = DF + M[i][j][5] * F[i][j+1];
   }
}

Le comportement est différent suivant le compilateur :

Sous le compilateur INTEL, la seconde boucle vas 1.7 fois plus vite en optimisation -O4, et vas 2.3 fois plus vite si on rajoute le SSE4.2. Pas mal du tout donc.

Sous GCC, la première boucle reste un peu plus rapide que la seconde, avec les mêmes options que le compilateur INTEL (SSE4.2, et -O4 et même avec le tuning pour mon XEON).

Donc c'est super sur compilateur INTEL, mais pourri sur GCC.


Auriez vous d'autres idées pour optimiser la boucle ?

Sachant que je m'intéresse à 2 cas :
1 - les coefficients M ne dépendent pas de i et j donc M n'est plus qu'un tableau de 5 valeurs. (Dérivée différence finie)
2 - les coefficients M dépendent de i et j. (Laplacien, Solveur de Poisson)

Merci d'avance pour vos idées

[fezvez]

Bonjour, c'est un problème intéressant, et je souhaiterais avoir un petit peu plus d'informations.

Tout d'abord, je me demande pourquoi on utilise le tableau M au lieu de juste 5 variables. Je veux dire, si on a en effet affaire à une dérivée seconde centrée discrétisée, alors, je ne vois pas pourquoi on utiliserai des coefficients différents selon la position dans l'image.

Ensuite, il y a quelque chose de plus "trivial". DF[i][j] ne dépend que de M et de F. Donc, on peut paralléliser le calcul. Avec la démocratisation des processeurs multi-coeur, ce n'est vraiment pas à négliger.

(En gros, l'idée est de dire : si j'ai N processeurs, je créé N thread
- pour j allant de 0 à 512/N, et i allant de 0 à 512, calcul sur le processeur 1
- pour j allant de 512/N +1 à 2*(512/N), et i allant de 0 à 512, calcul sur le processeur 2
...

C'est affreusement mal présenté ici, mais il y a de très bons tutos sur la parallélisation. Et j'insiste, les opérations matricielles sont très très bien parallélisables, les GPU nous le prouvent.)

Edit : Oups, je m'excuse, je n'avais pas lu la fin du précédent message sur M constant

Re-edit : Je me disais bien que j'avais une idée si M est constant!
On a une "astuce" si certains des coefficients de M sont égaux, et je pense qu'il y en a peut-être, vu qu'avec un peu de chance, on a

Code :

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

1
2
3
   0  1  0
M= 1 -4  1
   0  1  0

Bref, on suppose qu'il y a au moins 2 coefficients identiques parmi les 5 coefficients de M. On va supposer que ce l'index de ces coefficients sont a et b ( on a donc a et b compris entre 0 et 4)
Soit C la valeur du coefficient. On a donc M[a] = M[b] = C

On va calculer la matrice Fc, qui est tout simplement la matrice F où chaque coefficient a été multiplié par C

Puis, il suffit de remplacer la multiplication M[a]*F[i][j] par Fc[i][j] ET M[b]*F[i][j] par Fc[i][j]

En un mot, on évite des calcules redondants, mais ça a un coût mémoire (logique)

[méphistopheles]

Je ne sais pas si je suis le seul, mais quelque chose m'étonne dans la code :
alors que les boucles vont de 0 à 512, on trouve des

Code :

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

[...] F[i-1][j] [...] F[i][j-1] [...]

donc des acces en -1.... je suis étonné qu'il n'y ai pas de segmentation fault.....

est-ce que c'est juste moi ?