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

[ElGat0]

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)

[moomba]

Merci pour cette réponse

En fait, il y a deux endroits dans le code qui ont un comportement "similaire" :

• Les dérivées avec 5 coefficients ou plus (avec un stencil plus grand si plus précises, donc des sauts en mémoire énormes) qui sont beaucoup utilisées.
• Le solveur de Poisson (laplacien) qui n'est autre qu'un BICGSTAB, une méthode itérative utilisant à outrance la boucle avec M variable (dans sa boucle principale, mais aussi dans le pré-conditionnement). La résolution de Poisson est la partie la plus couteuse du code si l'on met de côté la chimie.

Pour ce qui est de la parallélisation, et bien j'ai déjà parallélisé le code J'utilise MPI aussi bien en local que sur le cluster. Cependant, je n'utilise aucunes librairies externes pour le moment. Mon directeur de thèse tient à ce que le code soit portable, ou du moins dépendant uniquement des librairies classiques (BLAS, LAPACK, ESSL, FFTW, etc).


L'idée du champ Fc est futée.
Cela vas poser problème lorsque le pas d'espace sera variable, mais pour le moment, si l'on considère que 2 valeurs sont identiques (sur x, puis sur y), le gain pourrait être intéressant. Cependant j'ai peur que l'ajout d'un nouveau tableau engorge le cache et réduise les perfs.
Il faut que je teste pour savoir
Je regarde dés que j'ai un peu de temps, et je poste les résultats.

J'avais plutôt pensé à répartir le calcul sur des blocs de mémoire (par exemple ne prendre que la partie 128*128 du tableau 512*512) afin que toutes les données soient dans la cache, puis faire le calcul du bloc, et prendre le bloc suivant, etc. On évite ainsi les appels à un morceau hors cache et donc la latence induite.

[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 ?

[fezvez]

Une ch'tite réponse rapide (parce que le week-end arrive!)

Il a été dit que :

F est de taille 514*514

En gros, On entoure F d'une "bande de protection" justement pour éviter les segfaults. Et traiter les effets de bord!

[moomba]

Une réponse plus complète (mais la réponse de fezvez est la bonne ).

En CFD, on utilise des champs avec des "ghosts", surtout quand on découpe le champ physique sur plusieurs processeurs (MPI).
En bref, lorsque l'on calcule la dérivée sur les points, pour calculer la dérivée en i,j par exemple, il faut connaitre les points en i-2,i-1,i,i+1,i+2 pour une dérivée à l'ordre 4, et plus si on augmente la précision.

Or, comment calculer les points aux bords du domaine ?

Soit on diminue l'ordre et on décentre la dérivée (pour les conditions aux limites "physiques", on utilise par exemple i,i+1,i+2,i+3), soit on vas piocher dans le ghost. Ce fameux ghost correspond aux valeurs du processeur voisin qui calcule une autre partie du domaine physique. En fait, juste avant le calcul de la dérivée, il y a des communications entre processeurs. Pour un cas 1D avec 512 points (donc -1 à 514 si on considère une DF4) le proc 0 vas donner ses points 511 et 512 au proc 1 qui vas les mettre dans son ghost -1 et 0. Le 1 lui envoie en retour ses points 1 et 2, et 0 les stocks dans son ghost 513 et 514. Il reste à calculer la dérivée de 1 à 512

En espérant que c'est pas trop fouillis

[Joel F]

Les FFT ne servent à rien pour des images < 1024*1024 ou des masques < 15*15

Ensuite, le vrai probleme est ici :
- de maximiser la localité spatiale et temporelle du cache
- vectoriser le m-add de la convolution

Pour la GPU, pas la peine, le probleme est memory-bound, t'auras un gain ridicule (pour rester poli).

Quelques idées ici :
http://www.hipeac.net/system/files?f...ec_revised.pdf

le papier est ciblé CELL mais les optimsiations présentées sont valables sur toutes machines avec un cache.

[moomba]

Désolé du retard, j'étais en formation

fezvez
Bon, alors j'ai testé l'astuce de calculer les coefs qu'une seule fois, en supposant qu'ils ne dépendent pas de i et j.
J'ai mis le temps d'exécution moyen pour 5000 IT :

Cas normal : 3.2
Cas x /= y : 5.0
Cas x == Y : 3.0

Le derniers cas est très improbable car on ne simule que rarement des domaines carrés ou cubiques, et le gain n'est pas probant.

Donc comme on pouvais s'y attendre, le cache ne supporte pas l'ajout d'un tableau supplémentaire.

Joel F
Dans notre cas, le domaine dépasse rarement 512*512 sur un proc, ou alors 128*128*128. Le stencil quand à lui est de l'ordre de 6*6 à 9*9. Donc pas de FFT.

J'ai lue l'article, et il semble que l'optimisation porte surtout sur le "chaining" des opérations, et pas sur les opérations proprement dites. A moins que je me trompe bien sûr.
Ils cherchent en fait à cumuler certaines opérations pour conserver les valeurs en mémoire et saturer les pipes entre procs. Mais dans le code de calcul que j'utilise, il n'est pas possible de fusionner les dérivées avec autre chose... (même si c'était possible, c'est interdit)


J'ai testé l'astuce d'INTEL sur les POWER 6 du super calculateur, et malheureusement ces derniers ne supportent pas le SSE, donc pas de gain.

On se retrouve au poins de départ
Il faut trouver une autre astuce... Si vous avez d'autres idées, je suis prêt à les tester

[Joel F]

sur Power6, y a Altivec qui a l abonne idée de forunri un vec_madd en 2 cycles/vecteurs

et l'astuce que je preconisais n'est pas de pipeliner mais de reduire les chargements dans le cache via l'optimsiation par rotation de registres

[moomba]

reduire les chargements dans le cache via l'optimsiation par rotation de registres

Heuuu... Moi y'en a pas comprendre
Tu veut dire retourner le tableau afin que la mémoire reste continue en lecture ? Mais retourner les registres, ca doit couter un max non ?

Je jette un coup d'oeil du côté d'Altivec sinon

[Joel F]

Non. Relis le papier sur la partie optimisation des filtres.

[moomba]

Vu

Je teste, et je te tient au courant
C'est un peu long car je doit aussi optimiser la chimie du code (avec MKL et MASS), et ca me prend du temps

A noter pour les intéressés que l'astuce d'INTEL donnée plus haut fonctionne en fait aussi sur POWERPC pourvus que l'on utilise le -O4. On a dont un facteur 4 sur les machines INTEL avec IFORT (à cause du xsse4.2) et un facteur 2 sur les machines IBM. Je n'ai pas encore testé avec Altivec.

[moomba]

Désolé du retard.

On a un gros bug bien méchant dans le code, et ça fait plusieurs jours que je cours après...

Je propose de "freezer" le topic temporairement, le temps que je puisse reprendre l'optim. Je t'enverrais un MP dés que j'ai testé la méthode de retournement des registres