Comment extraire les indices de tous les éléments non nuls d'un tableau 0-1 en utilisant OpenMP ?

Question

Je suis un débutant en OpenMP. Je rencontre un tel problème.

J'ai un tableau de masques M avec la longueur N dont l'élément est either 0 or 1 . J'espère pouvoir extraire tous les indices i qui satisfait M[i]=1 et les stocker dans un nouveau tableau T .

Ce problème peut-il être accéléré par OpenMP ?

J'ai essayé le code suivant. Mais il n'est pas performant.

int count = 0;
#pragma omp parallel for
for(int i = 0; i < N; ++i) {
    if(M[i] == hashtag) {
        int pos = 0;
        #pragma omp critical (c1)
        pos = count++;
        T[pos] = i;
}

Answer 1

Je ne suis pas sûr à 100% que ce sera beaucoup mieux, mais vous pouvez essayer ce qui suit :

int count = 0;
#pragma omp parallel for
for(int i = 0; i < N; ++i) {
    if(M[i]) {
        #pragma omp atomic
        T[count++] = i;
    }
}

Si le tableau est assez clairsemé, les threads pourront parcourir un grand nombre de zéros sans attendre les autres. Mais vous ne pouvez mettre à jour qu'un seul indice à la fois. Le problème est que différents threads écrivent dans le même bloc de mémoire ( T ), ce qui signifie que vous rencontrerez des problèmes de mise en cache : chaque fois qu'un thread écrit dans le fichier T En effet, le cache de tous les autres cœurs est "sale", si bien que lorsqu'ils essaient de le modifier, il se passe beaucoup de choses en coulisse. Tout cela est transparent pour vous (vous n'avez pas besoin d'écrire de code pour le gérer) mais cela ralentit considérablement les choses - je soupçonne que c'est là votre véritable goulot d'étranglement. Si votre matrice est suffisamment grande pour que cela en vaille la peine, vous pouvez essayer de faire ce qui suit :

1. Créez autant de tableaux T qu'il y a de fils.
2. Chaque thread met à jour sa propre version de T
3. Combinez tous les tableaux T en un seul après la fin des boucles.

Cela pourrait être plus rapide (parce que les différents threads n'écrivent pas dans la même mémoire) - mais comme il y a si peu d'instructions dans la boucle, je pense que ce ne sera pas le cas.

EDITAR J'ai créé un programme de test complet, et j'ai constaté deux choses. Premièrement, le atomic ne fonctionne pas dans toutes les versions de omp et vous devrez peut-être utiliser T[count++] += i; pour qu'il compile (ce qui n'est pas grave puisque T peut être mis à zéro au départ) ; plus troublant, vous n'obtiendrez pas deux fois la même réponse si vous faites cela (la valeur finale de count change d'une passe à l'autre) ; si vous utilisez la fonction critical ce n'est pas le cas.

Une deuxième observation est que la vitesse du programme ralentit vraiment quand on augmente le nombre de threads, ce qui confirme ce que je soupçonnais à propos de la mémoire partagée (temps pour 10M d'éléments traités) :

threads  elapsed
  1        0.09s
  2        0.73s
  3        1.21s
  4        1.62s
  5        2.34s

Vous pouvez voir que c'est vrai en changeant la façon dont la matrice clairsemée M c'est que lorsque je crée M comme un tableau aléatoire, et tester pour M[i] < 0.01 * RAND_MAX (matrice dense à 0,1 %), les choses se déroulent beaucoup plus rapidement que si je la densifie à 10 % - ce qui montre que la partie située à l'intérieur de la matrice critical La section nous ralentit vraiment.

Cela étant, je ne pense pas qu'il y ait un moyen d'accélérer cette tâche dans OMP - le travail de consolidation des sorties de tous les threads en une seule liste à la fin va simplement consommer tout avantage de vitesse que vous auriez pu avoir, étant donné le peu de choses qui se passent à l'intérieur de la boucle interne. Donc, plutôt que d'utiliser plusieurs threads, je vous suggère de réécrire la boucle de la manière la plus efficace possible - par exemple :

for( i = 0; i < N; i++) {
  T[count] = i;
  count += M[i];
}

Dans mon benchmark rapide, cette solution était plus rapide que la solution OMP - comparable au threads = 1 solution. Encore une fois, cela est dû à la façon dont on accède à la mémoire ici. Notez que j'évite d'utiliser un if ce qui permet au code d'être aussi rapide que possible. Au lieu de cela, je profite du fait que M[i] est toujours égal à zéro ou un. À la fin de la boucle, vous devez éliminer l'élément T[count] parce qu'il sera invalide... les "bons éléments" sont T[0] ... T[count-1] . Un tableau M avec 10M d'éléments a été traité par cette boucle en ~ 0.02 sec sur ma machine. Cela devrait être suffisant pour la plupart des besoins ?

Answer 2

En me basant sur la fonction rapide de Floris, j'ai essayé de voir si je pouvais trouver un moyen de trouver une solution plus rapide avec OpenMP. J'ai trouvé deux fonctions foo_v2 y foo_v3 qui sont plus rapides pour les tableaux plus grands, foo_v2 est plus rapidement indépendant de la densité et foo_v3 est plus rapide pour les tableaux plus épars. La fonction foo_v2 crée essentiellement un tableau 2D avec une largeur N*nthreads ainsi qu'un tableau countsa qui contient les comptes pour chaque fil. Ceci est mieux expliqué avec du code. Le code suivant bouclerait sur tous les éléments écrits sur T.

for(int ithread=0; ithread<nthreads; ithread++) {
    for(int i=0; i<counta[ithread]; i++) {
        T[ithread*N/nthread + i]
    } 
}

La fonction foo_v3 crée un tableau 1D comme demandé. Dans tous les cas N doit être assez grande pour surmonter le surcoût de l'OpenMP. Le code ci-dessous prend par défaut 256MB avec une densité de M environ 10%. Les fonctions OpenMP sont toutes deux plus rapides d'un facteur 2 sur mon système Sandy Bridge à 4 cœurs. Si vous mettez la densité à 50% foo_v2 est encore plus rapide d'un facteur 2 environ mais foo_v3 n'est plus aussi rapide.

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

int foo_v1(int *M, int *T, const int N) {   
    int count = 0;
    for(int i = 0; i<N; i++) {
        T[count] = i;
        count += M[i];
    }
    return count;
}

int foo_v2(int *M, int *T, int *&counta, const int N) {
    int nthreads;
    #pragma omp parallel
    {
        nthreads = omp_get_num_threads();
        const int ithread = omp_get_thread_num();   
        #pragma omp single
        counta = new int[nthreads];
        int count_private = 0;      
        #pragma omp for
        for(int i = 0; i<N; i++) {
            T[ithread*N/nthreads + count_private] = i;
            count_private += M[i];
        }
        counta[ithread] = count_private;
    }
    return nthreads;
}

int foo_v3(int *M, int *T, const int N) {
    int count = 0;

    int *counta = 0;    
    #pragma omp parallel reduction(+:count)
    {
        const int nthreads = omp_get_num_threads();
        const int ithread = omp_get_thread_num();   
        #pragma omp single
        {
            counta = new int[nthreads+1];
            counta[0] = 0;
        }
        int *Tprivate = new int[N/nthreads];

        int count_private = 0;      
        #pragma omp for nowait
        for(int i = 0; i<N; i++) {
            Tprivate[count_private] = i;
            count_private += M[i];
        }
        counta[ithread+1] = count_private;
        count += count_private;

        #pragma omp barrier
        int offset = 0;
        for(int i=0; i<(ithread+1); i++) {
            offset += counta[i];
        }

        for(int i=0; i<count_private; i++) {
            T[offset + i] = Tprivate[i];
        }

        delete[] Tprivate;
    }
    delete[] counta;
    return count;
}

void compare(const int *T1, const int *T2, const int N, const int count, const int *counta, const int nthreads) {
    int diff = 0;
    int n = 0;
    for(int ithread=0; ithread<nthreads; ithread++) {
        for(int i=0; i<counta[ithread]; i++) {
            int i2 = N*ithread/nthreads+i;
            //printf("%d %d\n", T1[n], T2[i2]);
            int tmp = T1[n++] - T2[i2];
            if(tmp<0) tmp*=-1;
            diff += tmp;
        }
    }
    printf("diff %d\n", diff);
}

void compare_v2(const int *T1, const int *T2, const int count) {
    int diff = 0;
    int n = 0;
    for(int i=0; i<count; i++) {
        int tmp = T1[i] - T2[i];
        //if(tmp!=0) printf("%i %d %d\n", i, T1[i], T2[i]);
        if(tmp<0) tmp*=-1;
        diff += tmp;
    }
    printf("diff %d\n", diff);
}

int main() {
    const int N = 1 << 26;

    printf("%f MB\n", 4.0*N/1024/1024);
    int *M = new int[N];
    int *T1 = new int[N];
    int *T2 = new int[N];
    int *T3 = new int[N];

    int *counta;
    double dtime;
    for(int i=0; i<N; i++) {
        M[i] = ((rand()%10)==0);
    }

    //int repeat = 10000;
    int repeat = 1;
    int count1, count2;
    int nthreads;
    dtime = omp_get_wtime();
    for(int i=0; i<repeat; i++) count1 = foo_v1(M, T1, N);
    dtime = omp_get_wtime() - dtime;
    printf("time v1 %f\n", dtime);

    dtime = omp_get_wtime();
    for(int i=0; i<repeat; i++) nthreads = foo_v2(M, T2, counta, N);
    dtime = omp_get_wtime() - dtime;
    printf("time v2 %f\n", dtime);
    compare(T1, T2, N, count1, counta, nthreads);

    dtime = omp_get_wtime();
    for(int i=0; i<repeat; i++) count2 = foo_v3(M, T3, N);
    dtime = omp_get_wtime() - dtime;
    printf("time v2 %f\n", dtime);
    printf("count1 %d, count2 %d\n", count1, count2);
    compare_v2(T1, T3, count1); 

}

Answer 3

El critique L'opération doit être atomic au lieu de critical ; en fait, dans votre cas, vous devez utiliser l'option atomic capture clause :

int pos, count = 0;                     // pos declared outside the loop
#pragma omp parallel for private(pos)   // and privatized, count is implicitly
for(int i = 0; i < N; ++i) {            // shared by all the threads
    if(M[i]) {
        #pragma omp atomic capture
            pos = count++;
        T[pos] = i;
    }
}

Jetez un coup d'œil à cette réponse d'avoir une vue d'ensemble sur toutes les possibilités de atomic avec OpenMP.