Je veux optimiser cette boucle courte

Je voudrais optimiser cette boucle simple:

unsigned int i; while(j-- != 0){ //j is an unsigned int with a start value of about N = 36.000.000 float sub = 0; i=1; unsigned int c = j+s[1]; while(c < N) { sub += d[i][j]*x[c];//d[][] and x[] are arrays of float i++; c = j+s[i];// s[] is an array of unsigned int with 6 entries. } x[j] -= sub; // only one memory-write per j } 

La boucle a un temps d’exécution d’environ une seconde avec un bulldozer AMD à 4000 MHz. J’ai pensé à SIMD et à OpenMP (que j’utilise normalement pour obtenir plus de vitesse), mais cette boucle est récursive.

Aucune suggestion?

Je suis d’accord avec la transposition pour une meilleure mise en cache (mais voir mes commentaires à ce sujet à la fin), et il y a encore beaucoup à faire, voyons ce que nous pouvons faire avec la fonction complète …

Fonction originale, pour référence (avec quelques rangements pour ma santé mentale):

 void MultiDiagonalSymmesortingcMasortingx::CholeskyBackSolve(float *x, float *b){ //We want to solve LD Lt x = b where D is a diagonal masortingx described by Diagonals[0] and L is a unit lower sortingagular masortingx described by the rest of the diagonals. //Let D Lt x = y. Then, first solve L y = b. float *y = new float[n]; float **d = IncompleteCholeskyFactorization->Diagonals; unsigned int *s = IncompleteCholeskyFactorization->StartRows; unsigned int M = IncompleteCholeskyFactorization->m; unsigned int N = IncompleteCholeskyFactorization->n; unsigned int i, j; for(j = 0; j != N; j++){ float sub = 0; for(i = 1; i != M; i++){ int c = (int)j - (int)s[i]; if(c < 0) break; if(c==j) { sub += d[i][c]*b[c]; } else { sub += d[i][c]*y[c]; } } y[j] = b[j] - sub; } //Now, solve x from D Lt x = y -> Lt x = D^-1 y // Took this one out of the while, so it can be parallelized now, which speeds up, because division is expensive #pragma omp parallel for for(j = 0; j < N; j++){ x[j] = y[j]/d[0][j]; } while(j-- != 0){ float sub = 0; for(i = 1; i != M; i++){ if(j + s[i] >= N) break; sub += d[i][j]*x[j + s[i]]; } x[j] -= sub; } delete[] y; } 

En raison du commentaire sur la division parallèle donnant un gain de vitesse (bien que n’étant que O (N)), je suppose que la fonction elle-même s’appelle beaucoup. Alors, pourquoi allouer de la mémoire? x simplement x comme __ressortingct__ et changez y en x partout ( __ressortingct__ est une extension GCC, tirée de C99. Vous pouvez utiliser une define pour elle. Peut-être que la bibliothèque en a déjà une).

De même, bien que je suppose que vous ne pouvez pas modifier la signature, vous pouvez faire en sorte que la fonction prenne un seul paramètre et le modifie. b n’est jamais utilisé lorsque x ou y ont été définis. Cela signifie également que vous pouvez vous débarrasser de la twig dans la première boucle qui exécute ~ N * M fois. Utilisez memcpy au début si vous devez avoir 2 parameters.

Et pourquoi d un tableau de pointeurs? Doit-il être? Cela semble trop profond dans le code d’origine, donc je ne le toucherai pas, mais s’il y a une possibilité d’aplatir le tableau stocké, ce sera un gain de vitesse même si vous ne pouvez pas le transposer (multipliez, ajoutez, la déréférencement est plus rapide que déréférencement, append, déréférencement).

Donc, nouveau code:

 void MultiDiagonalSymmesortingcMasortingx::CholeskyBackSolve(float *__ressortingct__ x){ // comments removed so that suggestions are more visible. Don't remove them in the real code! // these definitions got long. Feel free to remove const; it does nothing for the optimiser const float *const __ressortingct__ *const __ressortingct__ d = IncompleteCholeskyFactorization->Diagonals; const unsigned int *const __ressortingct__ s = IncompleteCholeskyFactorization->StartRows; const unsigned int M = IncompleteCholeskyFactorization->m; const unsigned int N = IncompleteCholeskyFactorization->n; unsigned int i; unsigned int j; for(j = 0; j < N; j++){ // don't use != as an optimisation; compilers can do more with < float sub = 0; for(i = 1; i < M && j >= s[i]; i++){ const unsigned int c = j - s[i]; sub += d[i][c]*x[c]; } x[j] -= sub; } // Consider using processor-specific optimisations for this #pragma omp parallel for for(j = 0; j < N; j++){ x[j] /= d[0][j]; } for( j = N; (j --) > 0; ){ // changed for clarity float sub = 0; for(i = 1; i < M && j + s[i] < N; i++){ sub += d[i][j]*x[j + s[i]]; } x[j] -= sub; } } 

Eh bien, cela semble plus ordonné, et le manque d’allocation de mémoire et le nombre réduit de twigs, si ce n’est rien d’autre, sont un coup de pouce. Si vous pouvez modifier s pour inclure une valeur UINT_MAX supplémentaire à la fin, vous pouvez supprimer davantage de twigs (les deux contrôles i , qui exécutent à nouveau ~ N * M fois).

Maintenant, nous ne pouvons plus créer de boucles parallèles, ni combiner des boucles. Comme nous l’avons suggéré dans l’autre réponse, il faut maintenant réorganiser le boost. Sauf que… le travail requirejs pour réorganiser d présente exactement les mêmes problèmes de cache que le travail de création de la boucle. Et il faudrait de la mémoire allouée. Pas bon. Les seules options pour optimiser davantage sont les suivantes: modifier la structure de IncompleteCholeskyFactorization->Diagonals elle-même, ce qui impliquera probablement beaucoup de changements, ou trouver un algorithme différent qui fonctionne mieux avec les données dans cet ordre.

Si vous voulez aller plus loin, vos optimisations devront impacter une bonne partie du code (ce qui n’est pas une mauvaise chose; à moins que Diagonals ait une bonne raison de Diagonals un tableau de pointeurs, il semble que cela pourrait être utile avec un refactor).

pensez que vous souhaiterez peut-être transposer la masortingce d – signifie que vous la stockez de manière à pouvoir échanger les indices – faites de i l’indice extérieur:

  sub += d[j][i]*x[c]; 

au lieu de

  sub += d[i][j]*x[c]; 

Cela devrait améliorer les performances du cache.

Je souhaite donner une réponse à ma propre question: La mauvaise performance est due à des erreurs de conflit de cache dues au fait que Win7 (au moins) aligne de gros blocs de mémoire sur la même limite. Dans mon cas, pour tous les tampons, les adresses avaient le même alignement (bufferadress% 4096 était identique pour tous les tampons), de sorte qu’elles tombaient dans le même ensemble de caches que le cache L1. J’ai changé l’allocation de mémoire pour aligner les tampons sur différentes limites afin d’éviter les erreurs de conflit de cache et j’ai obtenu une accélération du facteur 2. Merci pour toutes les réponses, en particulier celles de Dave!