Optimiser une boucle 1D avec des SIMD

[moomba]

Salut à tous,

Je galère depuis plusieurs jours, et je vous expose donc le problème, au cas où l'un d'entre vous maitrise la chose : je cherche à optimiser une simple boucle en C, qui représente une diffusion de chaleur dans un milieu 1D. La version 3D viendra quand celle 1D sera fonctionelle.

Problème : ma boucle avec SSE est pratiquement 2 fois plus lente que celle sans :

#include <stdio.h>#include <emmintrin.h>#include <time.h>//#include <vector>#define n1 1004#define niter 200000int i,j,t;double U0[n1] __attribute__ ((aligned(16)));double U1[n1] __attribute__ ((aligned(16)));double Dx,Dy,Lx,Ly,InvDxDx,Dt,alpha,totaltime,Stab,DtAlpha,DxDx;__m128d vmmx00;__m128d vmmx01;__m128d vmmx02;__m128d vmmx10;__m128d va;__m128d vb;__m128d vc;__m128d vd;clock_t time0,time1;FILE *f1;int main(){/* ---- GENERAL ---- */  alpha = 0.4;  totaltime = 1.0/100.0;  Dt = totaltime/((niter-1)*1.0);  Lx = 1.0;  Dx = Lx/((n1-1)*1.0);  InvDxDx = 1.0/(Dx*Dx);  DxDx = Dx*Dx;  Stab = alpha*Dt*(InvDxDx);  DtAlpha = Dt*alpha;/* Stability if result <= 0.5 */  printf("Stability factor : %f \n",Stab);  for( i = 0; i < n1; i++){U0[i] = 0.0;}  U0[1] = 1.0;  U0[2] = 1.0;  U0[3] = 1.0;  U0[n1-2] = 2.0;//    for ( i = 0; i < n1; i++) {//	   for ( j = i + 1; j < n2; j++) { //		  std::swap(U0[i][j], U0[j][i]);  //	 }   //}   va = _mm_set1_pd(-2.0);   vb = _mm_set1_pd(InvDxDx);   vd = _mm_set1_pd(DtAlpha);time0=clock();for( t = 0; t < niter; t++){   for( i = 2; i < n1-2; i+=2)   {	  //printf("%d  %d  \n",i,j);	  //fflush(stdout);	  vmmx00 = _mm_load_pd(&U0[i]);	  vmmx01 = _mm_loadu_pd(&U0[i+1]);	  vmmx02 = _mm_loadu_pd(&U0[i-1]);	  vmmx10 = _mm_mul_pd(va,vmmx00); // U1[i][j] = -2.0*U0[i][j];	  vmmx10 = _mm_add_pd(vmmx10,vmmx01); // U1[i][j] = U1[i][j] + U0[i+1][j];	  vmmx10 = _mm_add_pd(vmmx10,vmmx02); // U1[i][j] = U1[i][j] + U0[i-1][j];	  vmmx10 = _mm_mul_pd(vb,vmmx10); // U1[i][j] = U1[i][j] * InvDxDx;	  vmmx10 = _mm_mul_pd(vd,vmmx10); // U1[i][j] = U1[i][j] * DtAlpha;	  vmmx10 = _mm_add_pd(vmmx10,vmmx00); // U1[i][j] = U1[i][j] + U0[i][j];	  _mm_store_pd(&U1[i],vmmx10);	  // U1[i][j] = U0[i][j] + DtAlpha*( (U0[i+1][j]-2.0*U0[i][j]+U0[i-1][j])*InvDxDx   }   for( i = 2; i < n1-2; i+=2)   {	  //printf("%d  %d  \n",i,j);	  //fflush(stdout);	  vmmx00 = _mm_load_pd(&U1[i]);	  vmmx01 = _mm_loadu_pd(&U1[i+1]);	  vmmx02 = _mm_loadu_pd(&U1[i-1]);	  vmmx10 = _mm_mul_pd(va,vmmx00); // U0[i][j] = -2.0*U1[i][j];	  vmmx10 = _mm_add_pd(vmmx10,vmmx01); // U0[i][j] = U0[i][j] + U1[i+1][j];	  vmmx10 = _mm_add_pd(vmmx10,vmmx02); // U0[i][j] = U0[i][j] + U1[i-1][j];	  vmmx10 = _mm_mul_pd(vb,vmmx10); // U0[i][j] = U0[i][j] * InvDxDx;	  vmmx10 = _mm_mul_pd(vd,vmmx10); // U0[i][j] = U0[i][j] * DtAlpha;	  vmmx10 = _mm_add_pd(vmmx10,vmmx00); // U0[i][j] = U0[i][j] + U1[i][j];	  _mm_store_pd(&U0[i],vmmx10);	  // U1[i][j] = U0[i][j] + DtAlpha*( (U0[i+1][j]-2.0*U0[i][j]+U0[i-1][j])*InvDxDx   }}time1=clock();printf("Loop 0, total time : %f \n", (double) time1-time0);f1 = fopen ("out0.dat", "wt");for( i = 1; i < n1-1; i++){   fprintf (f1, "%d\t%f\n", i, U0[i]);}// REF  for( i = 0; i < n1; i++){U0[i] = 0.0;}  U0[1] = 1.0;  U0[2] = 1.0;  U0[3] = 1.0;  U0[n1-2] = 2.0;time0=clock();for( t = 0; t < niter; t++){   for( i = 2; i < n1-2; i++)   {   U1[i] = U0[i] + DtAlpha* (U0[i+1]-2.0*U0[i]+U0[i-1])*InvDxDx;   }   for( i = 2; i < n1-2; i++)   {   U0[i] = U1[i] + DtAlpha* (U1[i+1]-2.0*U1[i]+U1[i-1])*InvDxDx;   }}time1=clock();printf("Loop 0, total time : %f \n", (double) time1-time0);f1 = fopen ("outref.dat", "wt");for( i = 1; i < n1-1; i++){   fprintf (f1, "%d\t%f\n", i, U0[i]);}}

Une idée du problème ? Compilé avec Gcc : gcc -msse2 -O3, sous windows avec un code 2.

[RinDman]

Pourquoi tu dois compiler avec SSE ? Vu que je n'étais qu'en deuxième année d' info, je ne vois pas l'intérêt de compiler avec le flag sse

Je te conseille de lancer le programme avec GDB, tu pourras voir l’entaille du code qui s’exécute et tu pourras avancer ligne PAR ligne

Il y a surement une boucle de merde quelque part ^^

Pour moi, tu peux mettre des conditions dans le FOR afin de ne pas à avoir parcourir toute la boucle ... Après vu le manque de commentaire ( description de l'algorythme ), je ne vois qu'une simple boucle

[moomba]

Le compilateur refuse la compilation sans les SSE si je les utilises. Pas le choix.

Le code s’exécute bien, je l'ai passé sous valgrind, mais ca reste lent.

L'algorithme est assez simple : ont itère sur une simple dérivée mathématique d'ordre 2, que l'on a modélisé (discrétisation) sous la forme : DDF = F[i-1]*A + F*B + F[i+1]*C. Le jeu est d’appliquer cette dérivée à outrance afin de converger vers un résultat stable, qui sera la solution à l'équilibre du système, qui n'est autre qu'un barreau froid, auquel on applique des extrémités chaudes. Il vas donc se réchauffer, et à la fin, on peut observer le profile de température dans le barreau. (l'exemple donné dans le code n'est pas physique, mais le principe y est).

Cette opération appliquée à des tableaux 3D de 512**3 prend un temps considérable, d'où un besoin d'optimiser. Je suis toujours dessus, mais ça prend du temps :)

Le but ici est d'avoir la boucle la plus rapide possible.

[Latios]

Tente avec gcc -funroll-loops ou -funroll-all-loops, tu devrais gagner du temps.

Sur mon laptop, sans :

Loop 0, total time : 1800000.000000

Loop 0, total time : 1890000.000000

avec -funroll-all-loops :

Loop 0, total time : 1200000.000000

Loop 0, total time : 1820000.000000

Ça marche bien pour la version SSE visiblement.

Après, tu peux essayer de paralléliser l'algorithme pour qu'il s'exécute sur plusieurs CPU (OpenMP peut peut-être aider là dessus). Je vois pas trop ce qu'on peut faire de plus.

Par contre, au total tu fais 2*niter itérations, non ? Puisque tu recalcules U0 et U1 à chaque tour de boucle.

EDIT : SSE c'est du SIMD normalement, j'ai pas l'impression que ton code en tire partie en fait...

[RinDman]

http://fr.wikipedia....SIMD_Extensions

En regardant les lignes de code faisant des additions vectorielles, je dirais pas non

Sinon au lieu d'utiliser les boucles simples, tu peux tenter une boucle récursives Mais ça risque de bouffer trop de ressources et ça peut empirer ...

[moomba]

Merci pour vos réponses. Je vais regarde l'impact du unroll. J'ai trouvé ce site aussi : http://www.virtualdub.org/blog/pivot/entry.php?id=307

Il y a ici un cas très similaire au mien, je vais donc partir de cette base, et apprendre à jongler avec les SIMD. Comme ça peut intéresser pas mal de monde, je vais utiliser ce topic comme "fil" de mon développement, et y mettre les étapes : de la découverte des instructions au code définitif.

Je suis certain qu'il y a moyen d'optimiser ces boucles, à voir cette article (et d'autres).

Par contre, au total tu fais 2*niter itérations, non ? Puisque tu recalcules U0 et U1 à chaque tour de boucle.

Exactement. C'est une technique de mise à jours des valeurs simples. U1 est en fait un simple tableau tampon.

Sinon au lieu d'utiliser les boucles simples, tu peux tenter une boucle récursives

Ca c'est gore . N'ayant jamais fait d'études en info, j'ai jamais osé y toucher.

Je commence en fin d'aprem a mettre les étapes.

[RinDman]

Non seulement ça peut être gore, mais ça peut en faire chialer plus d'un si on ne respecte pas les conditions d'arrêts et les cas non gérés

C'est utile dans certains cas spécifique