Décomposition de Cholesky avec OpenMP

J’ai un projet dans lequel nous résolvons l’inverse de grandes masortingces denses définies positives positives (plus de 3000 x 3000) à l’aide de la décomposition de Cholesky . Le projet est en Java et nous utilisons la bibliothèque CERN Colt BLAS . Le profilage du code montre que la décomposition de Cholesky est le goulot d’étranglement.

J’ai décidé d’essayer de paralléliser la décomposition de Cholesky à l’aide d’OpenMP et de l’utiliser comme DLL en Java (avec JNA). J’ai commencé avec le code de décomposition Cholesky en C de Rosetta Code .

Ce que j’ai remarqué, c’est que les valeurs d’une colonne, à l’exception de l’élément diagonal, sont indépendantes. J’ai donc décidé de calculer les éléments diagonaux en série et le rest des valeurs de la colonne en parallèle. J’ai également échangé l’ordre des boucles afin que la boucle interne passe par-dessus les lignes et la boucle externe par-dessus les colonnes. La version série est légèrement plus lente que celle de RosettaCode mais la version parallèle est six fois plus rapide que la version RosettaCode sur mon système à 4 coeurs (8 HT). L’utilisation de la DLL en Java accélère également nos résultats six fois. Voici le code:

double *cholesky(double *A, int n) { double *L = (double*)calloc(n * n, sizeof(double)); if (L == NULL) exit(EXIT_FAILURE); for (int j = 0; j <n; j++) { double s = 0; for (int k = 0; k < j; k++) { s += L[j * n + k] * L[j * n + k]; } L[j * n + j] = sqrt(A[j * n + j] - s); #pragma omp parallel for for (int i = j+1; i <n; i++) { double s = 0; for (int k = 0; k < j; k++) { s += L[i * n + k] * L[j * n + k]; } L[i * n + j] = (1.0 / L[j * n + j] * (A[i * n + j] - s)); } } return L; } 

Vous pouvez trouver le code complet pour le tester sur http://coliru.stacked-crooked.com/a/6f5750c20d456da9

Au départ, je pensais que le faux partage poserait un problème lorsque les éléments restants d’une colonne étaient petits comparés au nombre de threads, mais cela ne semble pas être le cas. j’ai essayé

 #pragma omp parallel for schedule(static, 8) // a cache line is 8 doubles 

Je n’ai pas trouvé d’exemples clairs sur la manière de paralléliser la décomposition de Choleskey. Je ne sais pas si ce que j’ai fait est idéal. Par exemple, cela fonctionnera-t-il bien sur un système NUMA?

Peut-être une approche basée sur les tâches est-elle meilleure en général? Les diapositives 7 à 9, à l’ adresse http://courses.engr.illinois.edu/cs554/fa2013/notes/07_cholesky.pdf , présentent un exemple de décomposition parallèle de Cholesky à l’aide de “tâches à granularité fine”. Je ne sais pas encore comment mettre cela en œuvre.

J’ai deux questions, spécifiques et générales. Avez-vous des suggestions sur la façon d’améliorer mon implémentation de Cholesky Decomposition with OpenMP? Pouvez-vous suggérer une implémentation différente de Cholesky Decomposition avec OpenMP, par exemple avec des tâches?

Edit: comme demandé ici est la fonction AVX que j’ai utilisée pour calculer s . Cela n’a pas aidé

 double inner_sum_AVX(double *li, double *lj, int n) { __m256d s4; int i; double s; s4 = _mm256_set1_pd(0.0); for (i = 0; i < (n & (-4)); i+=4) { __m256d li4, lj4; li4 = _mm256_loadu_pd(&li[i]); lj4 = _mm256_loadu_pd(&lj[i]); s4 = _mm256_add_pd(_mm256_mul_pd(li4, lj4), s4); } double out[4]; _mm256_storeu_pd(out, s4); s = out[0] + out[1] + out[2] + out[3]; for(;i<n; i++) { s += li[i]*lj[i]; } return s; } 

J’ai réussi à faire en sorte que SIMD fonctionne avec la décomposition de Cholesky. Je l’ai fait en utilisant la mosaïque de boucle comme je l’ai déjà utilisée dans la multiplication masortingcielle. La solution n’était pas anodine. Voici les heures pour une masortingce 5790×5790 sur mon système à 4 cœurs / 8 Ivy Bridge (eff = GFLOPS / (pic GFLOPS)):

 double floating point peak GFLOPS 118.1 1 thread time 36.32 s, GFLOPS 1.78, eff 1.5% 8 threads time 7.99 s, GFLOPS 8.10, eff 6.9% 4 threads+AVX time 1.36 s, GFLOPS 47.64, eff 40.3% 4 threads MKL time 0.68 s, GFLOPS 95.14, eff 80.6% // from LAPACKE_dpotrf single floating point peak GFLOPS 236.2 1 thread time 33.88 s, GFLOPS 1.91, eff 0.8% 8 threads time 4.74 s, GFLOPS 13.64, eff 5.8% 4 threads+AVX time 0.78 s, GFLOPS 82.61, eff 35.0% 

La nouvelle méthode est 25 fois plus rapide en double et 40 fois plus rapide en simple. L’efficacité est d’environ 35-40% du pic de FLOPS maintenant. Avec la multiplication masortingcielle, j’obtiens jusqu’à 70% avec AVX dans mon propre code. Je ne sais pas à quoi s’attendre de la décomposition de Cholesky. L’algorithme est partiellement série (lors du calcul du bloc diagonal, appelé sortingangle dans mon code ci-dessous), contrairement à la multiplication masortingcielle.

Mise à jour: Je suis dans un facteur pour 2 de la MKL. Je ne sais pas si je devrais être fier de cela ou gêné par cela, mais apparemment, mon code peut encore être amélioré de manière significative. J’ai trouvé une thèse de doctorat sur ce sujet qui montre que mon algorithme de bloc est une solution commune, j’ai donc réussi à réinventer la roue.

J’utilise des carreaux 32×32 pour les doubles et 64×64 pour les flotteurs. Je réordonne également la mémoire pour que chaque tuile soit contiguë et soit sa transposée. J’ai défini une nouvelle fonction de production de masortingce. La multiplication masortingcielle est définie comme suit:

 C_i,j = A_i,k * B_k,j //sum over k 

J’ai réalisé que dans l’algorithme de Cholesky, il y avait quelque chose de très similaire

 C_j,i = A_i,k * B_j,k //sum over k 

En écrivant la transposition des tuiles, j’ai pu utiliser presque exactement ma fonction optimzied pour la multiplication de masortingces (je n’avais à changer qu’une ligne de code). Voici la fonction principale:

 reorder(tmp,B,n2,bs); for(int j=0; j 

Voici les autres fonctions. J'ai six fonctions de produit pour SSE2, AVX et FMA, chacune avec une version double et une version flottante. Je montre seulement celui pour AVX et double ici:

 template  void sortingangle(Type *A, int n) { for (int j = 0; j < n; j++) { Type s = 0; for(int k=0; k void block(Type *A, Type *B, int n) { for (int j = 0; j  void reorder(Type *A, Type *B, int n, int bs) { int nb = n/bs; int ssortingde = bs*bs; //printf("%d %d %d\n", bs, nb, ssortingde); #pragma omp parallel for schedule(static) for(int i=0; i void reorder_inverse(Type *A, Type *B, int n, int bs) { int nb = n/bs; int ssortingde = bs*bs; //printf("%d %d %d\n", bs, nb, ssortingde); #pragma omp parallel for schedule(static) for(int i=0; i