Parallélisme : travaux dirigés et pratiques

4. OpenMP

4.1. Introduction

Le support d'OpenMP est pris en compte dans la presque totalité des compilateurs C++/Fortran. Par exemple GCC 6.1 supporte OpenMP 4.5.

On pourra consulter les très intéressantes vidéos de Tim Mattson (Senior Research Scientist chez Intel) qui constituent un tutoriel très didactique.

afin d'utiliser OpenMP :

on peut inclure le fichier entête #include <omp.h> si on désire utiliser les fonctions de la librairie
on doit réaliser l'édition de liens avec la librairie -fopenmp (gcc) ou -openmp selon les systèmes

4.2. Team of threads

The #pragma directive is the method specified by the C standard for providing additional information to the compiler, beyond what is conveyed in the language itself. Three forms of this directive (commonly known as pragmas) are specified by the 1999 C standard. A C compiler is free to attach any meaning it likes to other pragmas.

L'utilisation du parallélisme est réalisé grace aux pragma (directives de programmation). La pragma omp parallel indique à OpenMP de créer un ensemble de threads (qualifié de team) qui seront exécutés en parallèle.

CODE
omp_example_1

#include <omp.h>
#include <iostream>
#include <cstdlib>
using namespace std;
 
int main() {
 
    #pragma omp parallel
    {
            cout << "bonjour " << omp_get_thread_num() << endl;
    }
 
    return EXIT_SUCCESS;
}
 
 

Le résultat de l'exécution du programe précédent sur un Dual-Core avec Hyper-Threading est le suivant :

bonjour bonjour 03
bonjour 1

bonjour 2

On voit qu'il y a bien 4 threads 0,1,2,3 et que l'affichage est illisible étant donné que les threads partagent le même canal de communication cout.

Le modèle de fonctionnement est de type fork-join:

le thread qualifié de master exécute le code séquentiel
lorsqu'il rencontre une partie de code parallèle, il crée les threads nécessaires et exécute le code
la partie de code parallèle se termine lorsque tous les threads ont terminé leur exécution, il y a donc un join par défaut
puis le master thread reprend l'exécution séquentielle

4.3. Contrôle du nombre de threads

Plusieurs fonctions sont définies par OpenMP afin de manipuler le nombre total de threads :

omp_get_num_procs() retourne le nombre de processeurs du système; pour un quad-core (sans HyperThreading) ce nombre est égal à 4
omp_get_max_threads() retourne le nombre maximum de threads utilisables; pour un quad-core (sans HyperThreading) ce nombre est égal à 4
omp_get_num_threads() retourne le nombre de threads dans une équipe (team); pour une partie de code séquentiel, ce nombre est toujours égal à 1
omp_get_thread_num() retourne l'identifiant du thread dans la partie parallèle

CODE
omp_example_2

#include <omp.h>
#include <iostream>
#include <cstdlib>
using namespace std;
 
int main() {
 
    cout << "max procs = " << omp_get_num_procs() << endl;
    cout << "max threads = " << omp_get_max_threads() << endl;
    
    #pragma omp parallel
    {
            cout << "bonjour " << omp_get_thread_num() << "/" << omp_get_num_threads() << endl;
    }
 
    return EXIT_SUCCESS;
}
 
 

max procs = 4
max threads = 4
bonjour 1/4
bonjour 0/4
bonjour bonjour 2/4
3/4

4.4. Les sections critiques / synchronisation des threads

La synchronisation des threads est réalisée :

soit de manière critique (exclusion mutuelle) : #pragma omp critical, un thread à la fois exécute la zone critique
```
#pragma omp critical
{
	do_something(...);
}	
```
soit de manière atomique : #pragma omp atomic, un thread à la fois met à jour un emplacement mémoire
```
#pragma omp atomic
x = x + y;	
```
soit en attendant que l'ensemble des threads aient terminé leur travail #pragma omp barrier avant de continuer l'exécution du code à venir

soit en utilisant le mécanisme des verrous (lock) qui ne sera pas vu ici

omp_lock_t lock;
omp_init_lock(&lock);
#pragma omp parallel private (tmp, thread_id)
{
	thread_id = omp_get_thread_num();
	int tmp = f(thread_id);
	omp_set_lock(&lock);
	printf("%d %d", thread_id, tmp);
	omp_unset_lock(&lock);
}
omp_destroy_lock(&lock);

Comment utiliser OpenMP ?

Il suffit d'utiliser le parallélisme dans un des cadres suivants :

les boucles for
définir des sections de code parallèle
définir des tâches (task)

4.5. Les boucles for

Pour paralléliser les boucles on utilise #pragma omp for qui doit être introduite dans une section #pragma omp parallel. On peut également combiner les deux formes avec #pragma omp parallel for.

CODE
omp_example_3

#include <omp.h>
#include <iostream>
#include <cstdlib>
using namespace std;
 
int main() {
 
    cout << "max procs = " << omp_get_num_procs() << endl;
    cout << "max threads = " << omp_get_max_threads() << endl;
    
    const int SIZE = 100;
    int tab[SIZE];
    
    // initialize
    for (int i = 0; i < SIZE; ++i) {
        tab[i] = (i+1);
    }
    
    #pragma omp parallel
    {
        #pragma omp for
        for (int i = 0; i < SIZE; ++i) {
            #pragma omp critical
            {
                cout << "thread " << omp_get_thread_num() << " treats index " << i << endl;
            }
            tab[i] = 2 * i;
        }
    }
 
    return EXIT_SUCCESS;
}
 
 

4.5.1. Attributs de stockage des variables

On peut décrire la façon dont les variables sont accédées :

shared(var) : par défaut les variables sont partagées par tous les threads
private(var) désigne une variable privée de la section parallèle (à noter que la valeur originale déclarée dans la partie non parallèle du programme n'est pas recopiée)
firstprivate(var) semblable à private mais la valeur privée est initialisée à la dernière valeur dans la partie non parallélisée qui précède
lastprivate(var) à la fin de la section parallèle, la valeur que prend var pour le dernier thread sera affectée dans la section non parallèle qui suit

CODE
omp_example_4

#include <omp.h>
#include <iostream>
#include <cstdlib>
using namespace std;
 
int main() {
 
    int a[4], b , c, d;
    
    b = 1234;
    c = 5678;
    d = 91011;
    
    #pragma omp parallel for num_threads(4) shared(a) private(b) firstprivate(c) lastprivate(d)
    for (int i=0; i<4; ++i) {
        #pragma omp critical 
        {
            a[i] = (i + 1) * 5;
            b += a[ i ];
            c += a[ i ];
            d += a[ i ];
            cout << "execute thread " << omp_get_thread_num() << endl;
            cout << "b=" << b << endl;
            cout << "c=" << c << endl;
            cout << "d=" << d << endl; 
        }
    }
 
    cout << "after parallel section :" << endl;
    for (int i=0; i<4; ++i) {
        cout << "a[" << i << "]=" << a[i] << endl;
    }
    cout << "b=" << b << endl;
    cout << "c=" << c << endl;
    cout << "d=" << d << endl; 
    return EXIT_SUCCESS;
}
 
 

execute thread 0
b=-1097711931
c=5683
d=32772
execute thread 3
b=20
c=5698
d=20
execute thread 2
b=15
c=5693
d=15
execute thread 1
b=10
c=5688
d=10
after parallel section :
a[0]=5
a[1]=10
a[2]=15
a[3]=20
b=1234
c=5678
d=20

4.5.2. la sous-directive collapse

Specifies how many loops in a nested loop should be collapsed into one large iteration space and divided according to the schedule clause. The sequential execution of the iterations in all associated loops determines the order of the iterations in the collapsed iteration space.

Dans l'exemple suivant on a $5 x 10$ éléments à traiter suivant deux boucles imbriquées :

CODE
omp_collapse_1

#include <iostream>
#include <cstdlib>
using namespace std;
#include <omp.h>
 
int main() {
    const int MAX_I = 5;
    const int MAX_J = 10;
    const int MAX_THREADS = 4;
    int nbr_processed[MAX_THREADS];
    
    for (int i=0; i<MAX_THREADS; ++i) nbr_processed[i] = 0;
    
    #pragma omp parallel for 
    for (int i=0; i<MAX_I; ++i) {
        for (int j=0; j<MAX_J; ++j) {
            #pragma omp critical
            {
                cout << i << " " << j << " " << omp_get_thread_num() << endl;
                ++nbr_processed[omp_get_thread_num() ];
            }
        }
    }
    
    for (int i=0; i<MAX_THREADS; ++i) {
        cout << "thread " << i << " processed " <<  nbr_processed[i] << " elements" << endl;
    }
    
    return EXIT_SUCCESS;
}
 

Le résultat à l'exécution est le suivant : resultat. On dispose de 4 threads d'exécution et donc l'un thread va traiter $2$ fois $10$ éléments, soit deux fois plus que les autres.

thread 0 processed 20 elements
thread 1 processed 10 elements
thread 2 processed 10 elements
thread 3 processed 10 elements

On voit que la répartition de la charge de travail est déséquilibrée. Afin de mieux la répartir on peut considérer que les $5 x 10$ éléments à traiter ne font partie que d'une seule boucle. Dans ce cas, on utilise la sous-directive collapse :

CODE
omp_collapse_2

#include <iostream>
#include <cstdlib>
using namespace std;
#include <omp.h>
 
int main() {
    const int MAX_I = 5;
    const int MAX_J = 10;
    const int MAX_THREADS = 4;
    int nbr_processed[MAX_THREADS];
    
    for (int i=0; i<MAX_THREADS; ++i) nbr_processed[i] = 0;
    
    #pragma omp parallel for collapse(2)
    for (int i=0; i<MAX_I; ++i) {
        for (int j=0; j<MAX_J; ++j) {
            #pragma omp critical
            {
                cout << i << " " << j << " " << omp_get_thread_num() << endl;
                ++nbr_processed[omp_get_thread_num() ];
            }
        }
    }
    
    for (int i=0; i<MAX_THREADS; ++i) {
        cout << "thread " << i << " processed " <<  nbr_processed[i] << " elements" << endl;
    }
    
    return EXIT_SUCCESS;
}
 

thread 0 processed 13 elements
thread 1 processed 13 elements
thread 2 processed 12 elements
thread 3 processed 12 elements

Attention cependant avec des boucles du type suivant, collapse ne fonctionne pas, car dans la seconde boucle le nombre d'éléments à traiter n'est pas constant :


for (int i=1; i<MAX_I; ++i) {
	for (int j=0; j<i; ++j) {
		...
	}
}

4.6. Les sections

Les sections permettent d'exécuter du code en parallèle non lié à une boucle. Dans l'exemple on définit 3 sections indépendantes. Attention on distingue sections (au pluriel) qui engloble section (au singulier).


#pragma omp parallel sections
{
	#pragma omp section
	double a = compute_1();
	#pragma omp section
	double b = compute_2();
	#pragma omp section
	double c = compute_3();
}
double s = a + b + c;

4.7. Les taches

Le problème des sections est qu'elles ont tendance à ralentir l'exécution du code. En particulier si on dispose de 8 threads d'exécution et que seules deux sections sont exécutées en parallèle, alors les 6 autres threads sont en attente de la fin de l'exécution des sections. Pour remédier à ce problème on introduit les tâches avec OpenMP 3.0

La notion de tâche permet de paralléliser certains traitements qui étaient difficiles à paralléliser par le passé (boucle non contrainte, algorithme récursif). Elle repose sur le fait qu'un thread peut générer des tâches qui deviennent des threads pour être traités dans les boucles while par exemple. On introduit la notion de parentée entre les tâches liée à la synchronisation (voir ci-après).

Les tâches peuvent être exécutées immédiatement ou en différé.

La synchronisation des tâches est réalisée par :

#pragma omp barrier qui attend la fin de toutes les tâches
#pragma omp taskwait qui attend que les tâches filles (descendants directs) soient terminées avant de continuer l'exécution de la tâche parent


LinkedList *l = ...;

while (l != nullptr) {
	process(l);
	l = l->next;
}


LinkedList *l = ...;

#pragma omp parallel
{
	#pragma omp single
	{
		while (l != nullptr) {
			#pragma omp task firstprivate(l)
			{
				process(l);
			}
			l = l->next;
		}
	}
}

Ici firstprivate permet de rendre le pointeur l local à chaque thread.

4.7.1. Trop de tâches tuent le parallélisme

Attention, il faut pouvoir gérer les tâches de manière intelligente sous peine de perdre tout le bénéfice de la parallélisation.

Dans ce premier exemple beaucoup trop de tâches seront générées. La conséquence est un temps d'exécution anormalement long.

CODE
omp_fib_1

// Use of OpenMP tasks to compute fibonacci(45)
#include <iostream>
using namespace std;
 
 
int fib(int n) {
    int x, y;
    
    if (n < 2) return n;
    #pragma omp task shared(x)
    x = fib(n-1);
    #pragma omp task shared(y)
    y = fib(n-2);
    #pragma omp taskwait
    return x + y;
}
 
int main() {
    int res, n = 45;
 
    // parallel section but one thread will
    // lauch the function fib() 
    #pragma omp parallel 
    {
        #pragma omp single
        {
            res = fib(n);
        }
    }
    
    cout << "res=" << res << endl;
    return 0;
}
 

Il faut donc limiter le nombre de tâches :

CODE
omp_fib_2

// OpenMP task version of finoacci(45) with less tasks
#include <iostream>
using namespace std;
 
int fib(int n) {
    int x, y;
    
    if (n < 2) return n;
    if (n < 30) return fib(n-1) + fib(n-2);
    
    #pragma omp task shared(x)
    x = fib(n-1);
    #pragma omp task shared(y)
    y = fib(n-2);
    #pragma omp taskwait
    return x + y;
}
 
int main() {
    int res, n = 45;
    
    #pragma omp parallel 
    {
        #pragma omp single
        {
            res = fib(n);
        }
    }
    
    cout << "res=" << res << endl;
    return 0;
}
 

g++ -o omp_fib_1.exe omp_fib_1.cpp -O2

# exécution séquentielle sans OpenMP
time ./omp_fib_1.exe
res=1134903170

real	0m4.148s
user	0m4.135s
sys	0m0.010s

g++ -o omp_fib_1.exe omp_fib_1.cpp -O2 -fopenmp

# exécution parallèle avec OpenMP avec des milliers de tâches
time ./omp_fib_1.exe
???? plus de 56 minutes sans résultat


g++ -o omp_fib_2.exe omp_fib_2.cpp -O2 -fopenmp

# exécution parallèle avec OpenMP version améliorée
time ./omp_fib_2.exe
res=1134903170

real	0m3.336s
user	0m7.085s
sys	0m0.010s

Avec la première version (omp_fib_1.cpp) trop de tâches sont créées, par exemple pour fib(30), on a 2.692.537 tâches.

4.8. La vectorisation

La vectorisation a été introduite avec OpenMP 4.0. On pourra consulter à ce sujet

ce document
ou ce guide

CODE
omp_parallel_simd

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <xmmintrin.h>
using namespace std;
 
typedef float real;
 
const int TEN_MILLION = 10000000;
 
void process(real *x, real *y, real *z, int size) {
    #pragma omp parallel for simd schedule(static,10) 
    for (int i=0; i<size; ++i) {
        z[i] = (sin(x[i]) - cos(y[i])) / (sqrt(x[i] * x[i] + y[i] * y[i]));
    }
}
 
int main() {
    real *a, *b, *c;
    
    a = new real [TEN_MILLION];
    b = new real [TEN_MILLION];
    c = new real [TEN_MILLION];
    
    for (int i=0; i<TEN_MILLION; ++i) {
        a[i] = i % 1000;
    }
    for (int i=0; i<TEN_MILLION; ++i) {
        b[i] = i % 133;
    }
    
    process(a,b,c,TEN_MILLION);
    
    _mm_free(a);
    _mm_free(b);
    _mm_free(c);
    
    return EXIT_SUCCESS;
}
 
 

g++ -o omp_parallel_simd.exe omp_parallel_simd.cpp -O2
time ./omp_parallel_simd.exe 

real	0m1.541s
user	0m1.524s
sys	0m0.016s

g++ -o omp_parallel_simd.exe omp_parallel_simd.cpp -O2 -fopenmp
time ./omp_parallel_simd.exe 

real	0m0.482s
user	0m1.708s
sys	0m0.020s

Autre exemple:

CODE
omp_parallel_simd_prod

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <xmmintrin.h>
using namespace std;
#include <getopt.h>
 
struct option long_options[] {
    {"method", required_argument, 0, 'm'},
    {"size", required_argument, 0, 's'},
    {"iterations", required_argument, 0, 'i'},
    {0,0,0,0}
};
 
int method = 1;
int vsize = 10*1024*1024;
int iterations = 1;
    
    
    
typedef float real;
 
real prod1(real * x, real * y, real * z, int size) {
    real total = 0.0;
    for (int i=0; i<size; ++i) {
        z[i] = (z[i] + x[i]) / y[i];
        total += z[i];
//      cout << "z[" << i << "]=" << z[i] << " " << total << endl;
    }
    
    return total;
}
 
 
real prod2(real * x, real * y, real * z, int size) {
    real total = 0.0;
 
    #pragma omp parallel for reduction(+:total) 
    for (int i=0; i<size; ++i) {
        z[i] = (z[i] + x[i]) / y[i];
        total += z[i];
    }
    return total;
}
 
real prod3(real * x, real * y, real * z, int size) {
    real total = 0.0;
 
    #pragma omp simd reduction(+:total) 
    for (int i=0; i<size; ++i) {
        z[i] = (z[i] + x[i]) / y[i];
        total += z[i];
    }
    return total;
}
 
#pragma omp declare simd aligned(x,y,z:16) simdlen(4)
real prod4(real * x, real * y, real * z, int size) {
    real total = 0.0;
 
    #pragma omp parallel for simd reduction(+:total)
    for (int i=0; i<size; ++i) {
        z[i] = (z[i] + x[i]) / y[i];
        total += z[i];
    }
    
    return total;
}
 
typedef real (*ProductType)(real *x, real *y, real *z, int size);
 
extern "C" {
real prod_asm(real * x, real * y, real * z, int size);
real prod_asm_unroll2(real * x, real * y, real * z, int size);
real prod_asm_unroll4(real * x, real * y, real * z, int size);
}
 
ProductType products[] = {
    nullptr,
    prod1,
    prod2,
    prod3,
    prod4,
    prod_asm,
    prod_asm_unroll2,
    prod_asm_unroll4
};
 
int main(int argc, char *argv[]) {
    int cc, index;
    
    while ((cc = getopt_long(argc, argv, "m:s:i:", long_options, &index)) != -1) {
        switch (cc) {
            case 'm':
                method = atoi(optarg);
                break;
            case 's':
                vsize = atoi(optarg);
                break;
            case 'i':
                iterations = atoi(optarg);
                break;  
            default: /* '?' */
                exit(EXIT_FAILURE);
        }
    }
    real *a, *b, *c;
    
 
    a = (real *) _mm_malloc(vsize * sizeof(real), 16);
    b = (real *) _mm_malloc(vsize * sizeof(real), 16);
    c = (real *) _mm_malloc(vsize * sizeof(real), 16);
    
    for (int i=0; i<vsize; ++i) {
        //a[i] = 1 + (i % 10);
        a[i] = 0.2;
    }
    for (int i=0; i<vsize; ++i) {
        //b[i] = 1 + (i % 133);
        b[i] = 5;
    }
    for (int i=0; i<vsize; ++i) {
        c[i] = 0;
    }
    
    real result = 0;
 
    cout << "method=" << method << endl;    
    for (int k=0; k<iterations; ++k) {  
        cout << ".";
        result = (*products[method])(a,b,c,vsize);
    }
    cout << endl;
    cout << "result=" << std::fixed << result << endl;
    
    _mm_free(a);
    _mm_free(b);
    _mm_free(c);
    
    return EXIT_SUCCESS;
}
 
 

On dispose de 7 implantations de la méthode de traitement, à savoir :

méthode 1 : aucune optimisation
méthode 2 : réduction en parallèle avec OpenMP
méthode 3 : réduction vectorielle avec OpenMP
méthode 4 : traitement parallèle avec réduction vectorielle OpenMP
méthode 5 : vectorisation SSE en assembleur 64 bits
méthode 6 : vectorisation SSE assembleur 64 bits, dépliage de boucle de 2
méthode 7 : vectorisation SSE assembleur 64 bits, dépliage de boucle de 4

Intel(R) Core(TM) i5-7400 CPU @ 3.00GHz
Ubuntu 16.04.2 LTS (xenial)
g++ version 4.9.4
1000 itérations

float sse4.2
method ; result ; time
1 574328.500000 49.25
2 510569.562500 16.77
3 510569.562500 12.53
4 527646.750000 6.33
5 510569.562500 12.25
6 510569.562500 12.40
7 510569.562500 12.44

float avx2
method ; result ; time
1 574328.500000 13.23
2 510569.562500 6.32
3 529619.062500 6.43
4 522841.937500 6.04
5 510569.562500 12.31
6 510569.562500 12.35
7 510569.562500 12.33


double sse4.2
method ; result ; time
1 524287.999914 57.43
2 524288.000022 15.37
3 524287.999950 29.13
4 524288.000013 11.27
5 524287.999950 28.81
6 524287.999950 29.11
7 524287.999950 29.53

double avx2
method ; result ; time
1 524287.999914 16.38
2 524288.000022 12.08
3 524288.000022 13.79
4 524287.999995 12.14
5 524287.999950 29.39
6 524287.999950 29.46
7 524287.999950 29.12

On remarquera des résultats avec des erreurs d'arrondi et notamment pour la méthode 1 en float.

En outre le fait de compiler avec -msse4.2 ou -mavx2 gènère des temps de calcul très différents à l'exécution.

4.8.1. directives gérant la vectorisation

safelen(length) : nombre maximum d'itérations exécutées en parallèle (en pratique taille du vecteur, càd pour un registre SSE de 128 bits, il faut prendre 4)
aligned(list[:alignment]) : spécifie pour un tableau si son adresse est multiple d'une puissance de 2, en général on prend 16
simdlen(length) taille du registre vectoriel à utiliser

Vectorisation et réduction :

#pragma omp simd for reduction(+:sum)

4.9. Exercices

Exercice 4.1

Evaluer $π$ par la méthode de MonteCarlo qui consiste en une estimation statistique.

On considère le carré de côté 2 unités et le cercle de rayon 1 unité centré dans le carré. On génère des points dans le carré, seuls ceux dans le cercle seront utilisés pour évalué $π$ :

$$ \π = 4 × \text"points dans le cercle" / \text"points dans le carré" $$

Utiliser le code qui suit et écrire un programme mono thread. Mesurez le temps d'exécution si on estime $π$ avec $2 × 10^9$ points dans le carré.

CODE
omp_pi_montecarlo.cpp

 
#include <random>
#include <cstdint>
using namespace std;
 
typedef uint64_t u64;
 
// Fonction pour générer des points aléatoires dans le carré
// et compter ceux qui sont dans le cercle
u64 monte_carlo(u64 max_points_inside_square)
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> dis(-1.0, 1.0);
 
    u64 points_inside_circle = 0;
 
    for (u64 i = 0; i < max_points_inside_square; ++i)
    {
        // génère un point dans le carré
        double x = dis(gen);
        double y = dis(gen);
 
        // est-il dans le cercle ?
        if (x * x + y * y <= 1.0)
        {
            ++points_inside_circle;
        }
    }
 
    return points_inside_circle;
}

Exercice 4.2

Ecrire une version OpenMP afin d'évaluer $π$ en parallèle.

Mesurer le temps d'exécution avec 2, 4, 8 puis 16 threads.

Page de Jean-Michel Richer