3. C++ et le monde extérieur



sommaire chapitre 3

3.2. Programmation parallèle

3.2.1. Introduction

Depuis 2005 avec l'arrivée des Pentium D (Smithfield = 2 Pentium 4 Prescott), les processeurs d'architecture x86 disposent aujourd'hui pour la plupart de plusieurs coeurs de calcul. En réalisant l'exécution du code sur plusieurs processeurs ou plusieurs coeurs on diminue généralement le temps de calcul, il faut pour cela paralléliser nos algorithmes.

On parle de programmation parallèle ou concurrentielle.

La programmation parallèle et les architectures parallèles remontent cependant aux années 50, mais son utilisation tend à décourager le néophyte bien qu'on puisse parfois en une seule ligne de code (telles les #pragma d'OpenMP) introduire du parallèlisme dans un programme.

On peut prendre par analogie la caisse de supermarché :

parallèlisme

Les solutions existantes reposent principalement sur l'utilisation de threads si on utilise un processeur multi-coeurs.

Un thread ou fil (d'exécution) ou tâche est similaire à un processus car tous deux représentent l'exécution d'un ensemble d'instructions du langage machine d'un processeur. Du point de vue de l'utilisateur, ces exécutions semblent se dérouler en parallèle. Toutefois, là où chaque processus possède sa propre mémoire virtuelle, les threads d'un même processus se partagent sa mémoire virtuelle. Par contre, tous les threads possèdent leur propre pile d’appel. (Wikipedia)

Pour information, fin 2014, Intel lance le Xeon E5-2699v3 doté de 18 coeurs (36 threads) avec un cache L3 de 45 Mo, TDP 145 W, $4115.

3.2.2. Les différents cas de parallèlisme

situations de parallèlisme

calculs indépendants

Dans l'exemple qui suit, chaque calcul est indépendant des autres, on peut donc réaliser l'ensemble des affectations et calculs (corps de la boucle) en parallèle.

for (int i=0; i<N; ++i) {
    a[i] = f(i);
}

synchronisation

La synchronisation de threads est un mécanisme qui vise à bloquer l'exécution des différents threads à des points précis du programme de manière à ce que tous les threads passent les étapes bloquantes au moment prévu par le programmeur.

La synchronisation dans un environnement multi-threads vise notamment à assurer la cohérence des données.

coopération

La coopération entre threads est l'une des choses les plus difficile à réaliser. Elle nécessite d'utiliser la synchronisation ou la communication entre threads.

3.2.3. Les problèmes liés au parallèlisme

L'exclusion mutuelle

L'exclusion mutuelle correspond au cas où un seul thread doit exécuter une série d'instructions alors que les autres threads doivent attendre leur tour pour exécuter la même série d'instructions.

C'est généralement le cas lors des affectations.

Dans l'exemple suivant chaque thread augmente la contenu d'une variable de 1.

exclusion mutuelle

Quand ne pas utiliser le parallèlisme

On peut considérer qu'il n'est pas nécessaire d'utiliser le parallèlisme dans les cas suivants :

Dans l'exemple suivant, on vérifie que :

  1. #include <iostream>
  2. using namespace std;
  3. #include <cmath>
  4. #include <cstdlib>
  5.  
  6. float function(float x) {
  7.     return sin(x*x/(1+x));
  8. }
  9.  
  10. int test_no = 1;
  11.  
  12. // ----------------------------------------------
  13. // corps de boucle avec traitement rapide
  14. // ici une affectation
  15. // l'exécution est séquentielle
  16. // ----------------------------------------------
  17. void test1() {
  18.     const int size = 10000000;
  19.     float *tab = new float [ size ];
  20.    
  21.     for (int i=0; i<size; ++i) {
  22.         tab[i] = i;
  23.     }
  24.    
  25.     delete [] tab;
  26. }
  27.  
  28. // ----------------------------------------------
  29. // corps de boucle avec traitement rapide
  30. // ici une affectation
  31. // l'exécution est effectuée en parallèle
  32. // ----------------------------------------------
  33. void test2() {
  34.     const int size = 10000000;
  35.     float *tab = new float [ size ];
  36.  
  37.     #pragma omp parallel for num_threads(4)
  38.     for (int i=0; i<size; ++i) {
  39.         tab[i] = i;
  40.     }
  41.    
  42.     delete [] tab;
  43. }
  44.  
  45. // ----------------------------------------------
  46. // corps de boucle avec traitement long
  47. // ici calcul d'une fonction
  48. // l'exécution est séquentielle
  49. // ----------------------------------------------
  50. void test3() {
  51.     const int size = 10000000;
  52.     float *tab = new float [ size ];
  53.  
  54.     for (int i=0; i<size; ++i) {
  55.         tab[i] = function(i);
  56.     }
  57.    
  58.     delete [] tab;
  59. }
  60.  
  61. // ----------------------------------------------
  62. // corps de boucle avec traitement long
  63. // ici calcul d'une fonction
  64. // l'exécution est effectuée en parallèle
  65. // ----------------------------------------------
  66. void test4() {
  67.     const int size = 10000000;
  68.     float *tab = new float [ size ];
  69.  
  70.     #pragma omp parallel for num_threads(4)
  71.     for (int i=0; i<size; ++i) {
  72.         tab[i] = function(i);
  73.     }
  74.    
  75.     delete [] tab;
  76. }
  77.  
  78. // ----------------------------------------------
  79. // main function
  80. // ----------------------------------------------
  81. int main(int argc , char *argv[]) {
  82.     if (argc > 1) {
  83.         test_no = atoi(argv[1]);
  84.     }
  85.    
  86.     clock_t start = clock();
  87.    
  88.     int loop = 1000;
  89.    
  90.     while (--loop) {
  91.         switch(test_no) {
  92.             case 1: test1(); break;
  93.             case 2: test2(); break;
  94.             case 3: test3(); break;
  95.             case 4: test4(); break;
  96.         };
  97.     }
  98.    
  99.     clock_t stop = clock();
  100.  
  101.     float sec = static_cast<float>(stop - start) / 1000000;
  102.     cerr << "elapsed = " << sec << endl;
  103.     return 0;
  104. }
  105.  
  106.  
 test   temps (s)   options de compilation   traitement 
 test 1   11.45   -O2   court en série 
 test 1   4.39   -O2 -ftree-vectorize   court en série 
 test 2   6.11   -O2 -fopenmp   court en parallèle 
 test 3   302.66   -O2 -ftree-vectorize   long en série 
 test 4   66.02   -O2 -ftree-vectorize -fopenmp   long en parallèle 
 test 4   66.60   -O2 -fopenmp   long en parallèle 
temps d'exécution des tests sur Core i5 4570

3.3. OpenMP

OpenMP (Open Multi-Processing) est une interface de programmation pour le calcul parallèle pour les langages de programmation C, C++ et Fortran.

Sur le principe OpenMP repose sur un partage des variables associées au threads.

L'intérêt d'OpenMP est qu'il permet à moindre coût de parallèliser le code en ajoutant des directives de programmation au code existant notamment sur les boucles. Le compilateur gcc supporte OpenMP depuis la version 4.2.

OpenMP.org

"We have designed modern OpenMP to be the best parallel language for the three general purpose languages C, C++ and Fortran", said Michael Wong, OpenMP CEO. "It is a language that enables you to access all capabilities of your machine without dropping to another language. This Challenge will convince you of that or show us where we need to improve. Either way, the industry benefits." (17 Nov. 2014)

Voici un premier exemple pour comprendre : on dispose d'une boucle for à l'intérieur de laquelle on affiche l'indice de la boucle. On veut montrer comment paralléliser la boucle avec différentes technologies.

3.3.0.a  Code séquentiel : non parallélisé

  1. #include <iostream>
  2. using namespace std;
  3.  
  4. int main() {
  5.     for (int i=0; i<10; ++i) {
  6.         cerr << "i = " << i << endl;
  7.     }
  8.     return 0;
  9. }
  10.  
  11.  

Le code est exécuté sur un seul coeur, on obtient donc le résultat séquentiel suivant :

i = 0
i = 1
i = 2
i = 3
i = 4
i = 5
i = 6
i = 7
i = 8
i = 9

3.3.0.b  Code parallélisé mais affichage désordonné

On modifie le programme notamment en ajoutant la directive #pragma omp parallel for qui indique que le contenu de la boucle for doit être exécuté en parallèle. Le code est exécuté par 4 threads différents (0 à 3) car on utilise un processeur Quad Core.

On compile avec l'option -fopenmp

  1. #include <iostream>
  2. using namespace std;
  3. #include <omp.h>
  4.  
  5. int main() {
  6.     #pragma omp parallel for
  7.     for (int i=0; i<10; ++i) {
  8.         cerr << "i = " << i << ", thread " << omp_get_thread_num() << endl;
  9.     }
  10.     return 0;
  11. }
  12.  
  13.  

Cependant l'affichage est perturbé car les quatre threads partagent le même flux de sortie.

i = 0, thread 0i = 6, thread 2
i = 7, thread 2
i = 
i = 1, thread 0
i = 2, thread 0
i = 8, thread 3
i = 9, thread 3
3, thread 1
i = 4, thread 1
i = 5, thread 1

3.3.0.c  Code parallélisé et affichage ordonné

On modifie le code en ajoutant la directive #pragma omp critical au niveau de la partie de code dédiée à l'affichage qui doit être exécutée sans être interrompue par d'autres threads.

  1. #include <iostream>
  2. using namespace std;
  3. #include <omp.h>
  4.  
  5. int main() {
  6.     #pragma omp parallel for
  7.     for (int i=0; i<10; ++i) {
  8.         #pragma omp critical
  9.         {
  10.             cerr << "i = " << i << ", thread " << omp_get_thread_num() << endl;
  11.         }
  12.     }
  13.     return 0;
  14. }
  15.  
  16.  
i = 0, thread 0
i = 1, thread 0
i = 8, thread 3
i = 9, thread 3
i = 2, thread 0
i = 3, thread 1
i = 4, thread 1
i = 5, thread 1
i = 6, thread 2
i = 7, thread 2

On voit bien à présent la répartition des tâches :

On remarque également que le traitement des threads n'est pas séquentiel : thread 3 exécuté avant thread 1 par exemple.

3.3.0.d  Version thread C++11

On utilise la classe thread du C++11. La gestion des threads est laissée à l'utilisateur, il est nécessaire :

Compiler avec -std=c++11 -pthread

  1. #include <iostream>
  2. using namespace std;
  3. #include <thread>
  4.  
  5. // thread to execute
  6. void code(int i) {
  7.     cerr << "i = " << i << endl;
  8. }
  9.  
  10. // main
  11. int main() {
  12.     // array of threads
  13.     thread *t[10];
  14.    
  15.     // create threads
  16.     for (int i=0; i<10; ++i) {
  17.         t[i] = new thread(code,i);
  18.     }
  19.     return 0;
  20. }
  21.  
  22.  

3.3.0.e  Version Pthread (POSIX)

La gestion des threads est laissée à l'utilisateur : il faut créer les threads ainsi qu'une structure de données qui contient les données liées au thread.

Compiler avec -lpthread


  1. #include <iostream>
  2. using namespace std;
  3. #include <pthread.h>
  4.  
  5. // data structure used to pass argument to thread
  6. typedef struct {
  7.     int value;
  8. } thread_data;
  9.  
  10. // thread to execute
  11. void *code(void *arg) {
  12.     // get data
  13.     thread_data *data = (thread_data *) arg;
  14.     cerr << "i = " << data->value << endl;
  15. }
  16.  
  17. int main() {
  18.     // array of threads
  19.     pthread_t t[10];
  20.     // array of data for threads
  21.     thread_data td[10];
  22.    
  23.     // create threads
  24.     for (int i=0; i<10; ++i) {
  25.         td[i].value = i;
  26.         pthread_create(&t[i], NULL, code, (void *)&td[i]);
  27.     }
  28.     return 0;
  29. }
  30.  
  31.  

3.3.1. Contrôle du nombre de threads exécutant le code

On peut controler le nombre de threads qui exécutent la partie parallèle du programme :

double tab[100];
#pragma omp parallel num_threads(4)
{
	ind thread_id = omp_get_thread_num();
	process(thread_id, tab);
}

3.3.2. Synchronization des threads

La synchronisation est réalisée :

3.3.3. Combinaison de parallélisme et partage des tâches

Les codes suivants sont équivalents :

double tab[MAX];
int i;
#pragma omp parallel
{
	#pragma omp for
	for (i = 0; i<MAX; i++) {
		tab[i] = f(i);
	}
}	

double tab[MAX];
int i;
#pragma omp parallel for
for (i = 0; i<MAX; i++) {
	tab[i] = f(i);
}

3.3.4. Réduction

La réduction correspond au fait qu'une variable est critique au coeur de la boucle : cela concerne le calcul de la somme des éléments d'un tableau, la recherche du minimum ou du maximum, ...

double tab[MAX];
double sum = 0.0;
int i;

for (i = 0; i<MAX; i++) {
	sum += tab[i];
}

Dans ce cas on utilise la directive reduction(operation : list-of-variables), où

double tab[MAX];
double sum = 0.0;
int i;

#pragma omp parallel for reduction (+:sum)
for (i = 0; i<MAX; i++) {
	sum += tab[i];
}

3.4. MPI

MPI (The Message Passing Interface), conçue en 1993-94, est une norme définissant une bibliothèque de fonctions, utilisable avec les langages C, C++ et Fortran. Elle permet d'exploiter des ordinateurs distants ou multiprocesseur par passage de messages (Wikipedia).

Le terme passage de message signifie que l'on échange des données, notamment on transmet aux autres programmes le contenu de certaines variables.

Lors de l'exécution le même programme est envoyé sur différentes machines/processeurs (éventuellement sur la même machine si elle dispose d'un processeur avec plusieurs coeurs) qui échangeront le contenu de variables via le réseau.

L'ensemble des fonctions peut être trouvé à cette adresse MPI Doc

3.4.1. installation

sudo apt-get install libopenmpi-dev openmpi-bin openmpi-doc

3.4.2. compilation

Voici deux exemples simples qui affichent un message avec un temps de latence (instruction sleep) pour chaque thread / processeur. Le premier utilise des fonctions C, le second utilise des espaces de noms et des objets

  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <unistd.h> // for sleep
  4. using namespace std;
  5. #include <mpi.h>
  6.  
  7. // ==============================================================
  8. // version C
  9. // ==============================================================
  10.  
  11. int main(int argc, char ** argv) {
  12.     int proc_id, n_procs;
  13.     char name[80];
  14.     int length;
  15.  
  16.     // initialization
  17.     MPI_Init(&argc, &argv);
  18.     // get id of this processor/core
  19.     MPI_Comm_rank(MPI_COMM_WORLD,&proc_id);
  20.     // get number of processors/cores that are used
  21.     MPI_Comm_size(MPI_COMM_WORLD,&n_procs);
  22.     // get name of processor
  23.     MPI_Get_processor_name(name, &length);
  24.  
  25.     sleep(proc_id);
  26.    
  27.     cerr << "running on " << name << " with id = " << proc_id << " / " << n_procs;
  28.     cerr << " on " << name << endl;
  29.    
  30.     MPI_Finalize();
  31.  
  32.     exit(EXIT_SUCCESS);
  33. }
  34.  
  35.  
  1. #include <iostream>
  2. #include <cstdlib>
  3. #include <unistd.h> // for sleep
  4. using namespace std;
  5. #include <mpi.h>
  6.  
  7. // ==============================================================
  8. // version C++
  9. // ==============================================================
  10.  
  11. int main(int argc, char ** argv) {
  12.     int proc_id, n_procs;
  13.     char name[MPI_MAX_PROCESSOR_NAME];
  14.     int length;
  15.  
  16.     MPI::Init(argc, argv);
  17.     n_procs = MPI::COMM_WORLD.Get_size();
  18.     proc_id = MPI::COMM_WORLD.Get_rank();
  19.  
  20.     memset(name, 0, MPI_MAX_PROCESSOR_NAME);
  21.     MPI::Get_processor_name(name,length);
  22.  
  23.     sleep(proc_id);
  24.    
  25.     cerr << "running on " << name << " with id = " << proc_id << " / " << n_procs;
  26.     cerr << " on " << name << endl;
  27.    
  28.     MPI::Finalize();
  29.  
  30.     exit(EXIT_SUCCESS);
  31. }
  32.  
  33.  
mpic++ -o c3_mpi_ex_1.exe c3_mpi_ex_1.cpp

3.4.3. exécution

Pour exécuter un programme compilé avec MPI il faut utiliser mpirun en spécifiant le nombre de processeurs grâce à l'option -n ou -np suivant les systèmes :

mpirun -n 4 ./c3_mpi_ex_1.exe

Voici le résultat à l'affichage si on lance sur une seule machine Intel Core i3-2375M CPU @ 1.50GHz (Dual Core + HyperThreading, soit 4 threads) :

running on inspiron with id = 3 / 4running on inspiron with id = 0 / 4
running on running on inspiron with id = 2 / 4

Sur un cluster on utilisera par exemple (-pe = parallel environment):

qsub -pe mpi 16 test_mpi.sh

Pour obtenir une section critique pour l'affichage on utilise la fonction MPI_Barrier qui bloque l'appelant kusqu'à ce que tous les autres programmes aient appelé cette fonction :

  1. #include <iostream>
  2. #include <cstdlib>
  3. using namespace std;
  4. #include <mpi.h>
  5.  
  6. int main(int argc, char ** argv) {
  7.     int proc_id, n_procs;
  8.     char name[80];
  9.     int length;
  10.  
  11.     MPI_Init(&argc, &argv);
  12.  
  13.     // get id of this processor/core
  14.     MPI_Comm_rank(MPI_COMM_WORLD,&proc_id);
  15.     // get number of processors/cores that are used
  16.     MPI_Comm_size(MPI_COMM_WORLD,&n_procs);
  17.    
  18.     MPI_Get_processor_name(name, &length);
  19.  
  20.     // display with critical section
  21.     int id = 0;
  22.     while (id < n_procs) {
  23.         if (id == proc_id) {
  24.             cerr << "running on " << name << " with id = " << proc_id << " / " << n_procs << endl;
  25.         }
  26.         ++id;
  27.         MPI_Barrier(MPI_COMM_WORLD);
  28.     }
  29.    
  30.     MPI_Finalize();
  31.     exit(EXIT_SUCCESS);
  32. }
  33.  
  34.  
  35.  
  36.  
running on inspiron with id = 0 / 4
running on inspiron with id = 1 / 4
running on inspiron with id = 2 / 4
running on inspiron with id = 3 / 4

3.4.4. échange de donneés Send, Recv

Pour transmettre des données entre les processeurs on utilise deux fonctions MPI_Send et MPI_Recv :

int MPI_Send(
    void* data,
    int count,
    MPI_Datatype datatype,
    int destination,
    int tag,
    MPI_Comm communicator)
    
int MPI_Recv(
    void* data,
    int count,
    MPI_Datatype datatype,
    int source,
    int tag,
    MPI_Comm communicator,
    MPI_Status* status)

3.4.5. Réduction

  1. # include <cstdlib>
  2. # include <iostream>
  3. # include <iomanip>
  4. # include <ctime>
  5.  
  6. using namespace std;
  7.  
  8. # include "mpi.h"
  9.  
  10. int main ( int argc, char *argv[] );
  11. void timestamp ( );
  12.  
  13. //****************************************************************************80
  14.  
  15. int main ( int argc, char *argv[] )
  16.  
  17. //****************************************************************************80
  18. //
  19. //  Purpose:
  20. //
  21. //    MAIN is the main program for SUM.
  22. //
  23. //  Discussion:
  24. //
  25. //    SUM is an example of using the MPI message passing interface library.
  26. //
  27. //  Licensing:
  28. //
  29. //    This code is distributed under the GNU LGPL license.
  30. //
  31. //  Modified:
  32. //
  33. //    01 May 2003
  34. //
  35. //  Author:
  36. //
  37. //    John Burkardt
  38. //
  39. //  Reference:
  40. //
  41. //    Forrest Hoffman,
  42. //    Message Passing with MPI and PVM,
  43. //    LINUX Magazine,
  44. //    Volume 4, Number 4, April 2002, pages 38-41, 63.
  45. //
  46. //    William Gropp, Ewing Lusk, Anthony Skjellum,
  47. //    Using MPI: Portable Parallel Programming with the
  48. //    Message-Passing Interface,
  49. //    Second Edition,
  50. //    MIT Press, 1999,
  51. //    ISBN: 0262571323.
  52. //
  53. {
  54. # define N 100
  55.  
  56.   double array[N];
  57.   int i;
  58.   int id;
  59.   int p;
  60.   double PI = 3.141592653589793238462643;
  61.   double seed;
  62.   MPI::Status status;
  63.   double sum;
  64.   double sum_all;
  65.   int tag;
  66. //
  67. //  Initialize MPI.
  68. //
  69.   MPI::Init ( argc, argv );
  70. //
  71. //  Get the number of processes.
  72. //
  73.   p = MPI::COMM_WORLD.Get_size ( );
  74. //
  75. //  Determine the rank of this process.
  76. //
  77.   id = MPI::COMM_WORLD.Get_rank ( );
  78. //
  79. //  Say hello.
  80. //
  81.   if ( id == 0 )
  82.   {
  83.     timestamp ( );
  84.     cout << "\n";
  85.     cout << "SUM - Master process:\n";
  86.     cout << "  C++ version\n";
  87.     cout << "\n";
  88.     cout << "  An MPI example program.\n";
  89.     cout << "  The master process computes some coefficients,\n";
  90.     cout << "  sends them to each worker process, which sums them.\n";
  91.     cout << "\n";
  92.     cout << "  Compiled on " << __DATE__ << " at " << __TIME__ << ".\n";
  93.     cout << "\n";
  94.     cout << "  The number of processes available is " << p << "\n";
  95.   }
  96. //
  97. //  The master process initializes the array.
  98. //
  99.   if ( id == 0 )
  100.   {
  101.     seed = 1.2345;
  102.  
  103.     for ( i = 0; i < N; i++ )
  104.     {
  105.       array[i] = ( double ) i * seed * PI;
  106.     }
  107.   }
  108. //
  109. //  The master process broadcasts the computed initial values
  110. //  to all the other processes.
  111. //
  112.   MPI::COMM_WORLD.Bcast ( array, N, MPI::DOUBLE, 0 );
  113. //
  114. //  Each process adds up its entries.
  115. //
  116.   sum = 0.0;
  117.   for ( i = 0; i < N; i++ )
  118.   {
  119.     sum = sum + array[i] * ( double ) id;
  120.   }
  121.  
  122.   cout << "\n";
  123.   cout << "SUM - Process " << id << ":\n";
  124.   cout << "  My contribution to the sum is " << sum << "\n";
  125. //
  126. //  Each worker process sends its sum back to the master process.
  127. //
  128.   if ( id != 0 )
  129.   {
  130.     MPI::COMM_WORLD.Send ( &sum, 1, MPI::DOUBLE, 0, 1 );
  131.   }
  132.   else
  133.   {
  134.     sum_all = sum;
  135.     for ( i = 1; i < p; i++ )
  136.     {
  137.       tag = 1;
  138.       MPI::COMM_WORLD.Recv ( &sum, 1, MPI::DOUBLE, MPI::ANY_SOURCE, tag,
  139.         status );
  140.  
  141.       sum_all = sum_all + sum;
  142.     }
  143.   }
  144.   if ( id == 0 )
  145.   {
  146.     cout << "\n";        
  147.     cout << "SUM - Master process:\n";        
  148.     cout << "  The total sum is " << sum_all << "\n";
  149.   }
  150. //
  151. //  Terminate MPI.
  152. //
  153.   MPI::Finalize ( );
  154. //
  155. //  Terminate.
  156. //
  157.   if ( id == 0 )
  158.   {
  159.     cout << "\n";        
  160.     cout << "SUM - Master process:\n";      
  161.     cout << "  Normal end of execution.\n";
  162.     cout << "\n";
  163.     timestamp ( );      
  164.   }
  165.   return 0;
  166.  
  167. # undef N
  168. }
  169. //****************************************************************************80
  170.  
  171. void timestamp ( )
  172.  
  173. //****************************************************************************80
  174. //
  175. //  Purpose:
  176. //
  177. //    TIMESTAMP prints the current YMDHMS date as a time stamp.
  178. //
  179. //  Example:
  180. //
  181. //    31 May 2001 09:45:54 AM
  182. //
  183. //  Licensing:
  184. //
  185. //    This code is distributed under the GNU LGPL license.
  186. //
  187. //  Modified:
  188. //
  189. //    08 July 2009
  190. //
  191. //  Author:
  192. //
  193. //    John Burkardt
  194. //
  195. //  Parameters:
  196. //
  197. //    None
  198. //
  199. {
  200. # define TIME_SIZE 40
  201.  
  202.   static char time_buffer[TIME_SIZE];
  203.   const struct std::tm *tm_ptr;
  204.   size_t len;
  205.   std::time_t now;
  206.  
  207.   now = std::time ( NULL );
  208.   tm_ptr = std::localtime ( &now );
  209.  
  210.   len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );
  211.  
  212.   std::cout << time_buffer << "\n";
  213.  
  214.   return;
  215. # undef TIME_SIZE
  216. }
  217.  

3.5. Intel TBB (Threading Building Block)

Intel Threading Building Block est un framework C++ qui aide à la parallèlisation du code.

Voir le tutoriel

3.5.1. Hello world ! avec TBB

  1. #include "tbb/tbb.h"
  2. #include <iostream>
  3. using namespace tbb;
  4. using namespace std;
  5.  
  6. class first_task : public task {
  7.     public:
  8.     task* execute( ) {
  9.        cout << "Hello World!\n";
  10.        return NULL;
  11.     }
  12. };
  13.  
  14. int main( )
  15. {
  16.     task_scheduler_init init(task_scheduler_init::automatic);
  17.     first_task& f1 = *new(tbb::task::allocate_root()) first_task( );
  18.     tbb::task::spawn_root_and_wait(f1);
  19.    
  20.     return 0;
  21. }
  22.  
  23.  

3.5.1.a  Application d'une fonction sur les éléments d'un tableau : parallel_for

  1. #include "tbb/tbb.h"
  2. #include <iostream>
  3. using namespace tbb;
  4. using namespace std;
  5.  
  6. // ----------------------------------------------
  7. // task that multiplies by 2 the value of A[i]:
  8. // A[i] = 2*A[i] for all i
  9. // ----------------------------------------------
  10. class TaskMultiplyBy2 {
  11.     float *array;
  12. public:
  13.     // function to apply
  14.     void operator()(const blocked_range<size_t>& r ) const {
  15.         float *a = array;
  16.         for( size_t i=r.begin(); i!=r.end(); ++i ) {
  17.             a[i] *= 2;
  18.             cout << "** " << i << endl;
  19.          }
  20.     }
  21.    
  22.     // constructor
  23.     TaskMultiplyBy2(float *a) : array(a) {
  24.     }
  25. };  
  26.  
  27.  
  28. // ----------------------------------------------
  29. // main function
  30. // ----------------------------------------------
  31. int main( ) {
  32.     // create array of 10 integers
  33.     const int N = 10;
  34.     float *A = new float[ N ];
  35.     for (size_t i=0; i<N; ++i) A[i] = i;
  36.    
  37.     // perform computation in parallel
  38.     parallel_for(blocked_range<size_t>(0,N), TaskMultiplyBy2(A));
  39.    
  40.     // print result
  41.     for (size_t i=0; i<N; ++i) {
  42.         cout << A[i] << " ";
  43.     }
  44.     return 0;
  45. }
  46.  
  47.  
* ** 50
** 1
** 2
** 3
** 6
** 8
** 9

** 4
** 7
0 2 4 6 8 10 12 14 16 18

3.5.2. Somme des éléments d'un tableau : parallel_reduce

  1. #include "tbb/tbb.h"
  2. #include <iostream>
  3. using namespace tbb;
  4. using namespace std;
  5. #include <time.h>
  6. #include <math.h>
  7.  
  8. // ----------------------------------------------
  9. // Time
  10. // ----------------------------------------------
  11. double ts_start, ts_stop;
  12.  
  13. void chrono_start() {
  14.     struct timespec ts;
  15.     clock_gettime(CLOCK_REALTIME, &ts);
  16.     ts_start = ts.tv_sec + 1e-9*ts.tv_nsec;
  17. }
  18.  
  19. void chrono_stop() {
  20.     struct timespec ts;
  21.     clock_gettime(CLOCK_REALTIME, &ts);
  22.     ts_stop = ts.tv_sec + 1e-9*ts.tv_nsec;
  23. }
  24.  
  25.  
  26. double chrono_get_seconds() {
  27.     return (ts_stop - ts_start);
  28. }
  29.  
  30. float function(float x) {
  31.     return sin(x*x/(1+x));
  32. }
  33.  
  34. // ----------------------------------------------
  35. // Task
  36. // ----------------------------------------------
  37. class TaskSum {
  38. protected:
  39.     float *array;
  40.     float sum;
  41. public:
  42.     // initial constructor with array
  43.     TaskSum(float *a) : array(a), sum(0.0) {
  44.     }
  45.    
  46.     // constructor for subtask, use split and assign 0 to sum
  47.     // and distinguish from copy constructor
  48.     TaskSum(const TaskSum& t, split) : array(t.array), sum(0) {
  49.     }
  50.    
  51.     // function to apply
  52.     void operator()(const blocked_range<size_t>& r) {
  53.         float *a = array;
  54.         float local_sum = sum; // important see page 23 of tbb_toturial.pdf
  55.        
  56.         for(size_t i=r.begin(); i!=r.end(); ++i) {
  57.             local_sum += function(a[i]);  
  58.          }  
  59.          sum = local_sum;
  60.     }
  61.      
  62.     // join threads at the end
  63.     void join(const TaskSum& t) {
  64.         sum += t.sum;
  65.     }
  66.    
  67.     // return result
  68.     float get_sum() {
  69.         return sum;
  70.     }
  71. };  
  72.  
  73. // ----------------------------------------------
  74. // main function
  75. // ----------------------------------------------
  76. int main(int argc, char *argv[]) {
  77.     const int TEN_MILLION = 10000000;  
  78.     int N = TEN_MILLION;
  79.    
  80.     if (argc > 2) {
  81.         N = atoi(argv[2]);
  82.     }
  83.    
  84.     // create and fill array
  85.     float *A = new float[ N ];
  86.     for (size_t i=0; i<N; ++i) A[i] = 1.0;
  87.    
  88.     // perform treatment
  89.     if (atoi(argv[1]) == 1) {
  90.         // without TBB
  91.         float sum = 0.0;
  92.         chrono_start();
  93.         for (size_t i=0; i<N; ++i) {
  94.             sum += function(A[i]);
  95.         }
  96.         chrono_stop();
  97.         cerr << "without TBB:" << endl;
  98.         cerr << "N = " << N << endl;
  99.         cerr << "elapsed = " << chrono_get_seconds() << endl;
  100.         cerr << "sum = " << sum << endl;
  101.    
  102.     } else {
  103.         // with TBB
  104.         chrono_start();
  105.         TaskSum task(A);
  106.         parallel_reduce(blocked_range<size_t>(0,N), task);
  107.         chrono_stop();
  108.         cerr << "with TBB: " << endl;
  109.         cerr << "N = " << N << endl;
  110.         cout << "elapsed=" << chrono_get_seconds() << endl;
  111.         cout << "sum=" << task.get_sum() << endl;
  112.     }
  113.     return 0;
  114. }
  115.  
  116.  
> ./tbb_example_3.exe 1
without TBB:
N = 10000000
elapsed = 0.0130236
sum = 4.77843e+06

> ./tbb_example_3.exe 2
with TBB: 
N = 10000000
elapsed=0.00438952
sum=4.79802e+06

Sur le graphique suivant on peut voir que le parallèlisme devient intéressant à mesure que la taille des données augmente.

tbb example 3 results

temps d'exécution (s) en fonction de la taille des données

3.5.3. Un cas concret avec la phylogénie

La SBL (Standard Bioinformatics Library) est un package logiciel écrit par Jean-Michel Richer et qui offre un ensemble de classes C++ pour la gestion des séquences et la reconstruction phylogénétique avec Maximum de Parcimonie (Maximum Parsimony).

On utilise pour résoudre le problème du Maximum de Parcimonie des méta-heuristiques (Recherche Locale, Algorithmes génétiques, mémétiques, Recuit Simulé, ...) car il s'agit d'un problème d'Optimisation Combinatoire :

Le but du problème est de trouver le ou les arbres pour lesquels la fonction objectif est minimale. C'est à dire l'arbre le plus parcimonieux ou encore l'abre qui possède le moins de mutations.

Le voisinage SPR consiste à dégrapher un sous-arbre $S$ de l'arbre initial $T$ et à le regrapher sur un autre noeud de l'arbre résultant du dégraphe $(T-S)$ de manière à ce que le nouvel arbre obtenu ait un score inférieur à $T$.

Cette opération peut être réalisée en parallèle : pour être plus efficace il ne faut pas regrapher $S$ sur $T-S$ mais réaliser une évaluation du regraphe.

La modification du code est réduite puisqu'il suffit d'ajouter un #pragma omp parallel for devant la boucle à paralléliser.

$T^{⋆} = T$ // meilleur arbre de score minimum

1. pour chaque sous-arbre $S$ de l'arbre $T$ faire
    $T^' = T-S$
    #pragma omp parallel for
    2. pour chaque sous-arbre $U$ de $T^'$ faire
    	$T^{⋄} = T^' +  U ⋈ S$ // regraphe de $S$ sur $U$
    	#pragma omp critical
        si $ score(T^{⋄}) < score(T^{⋆})$ alors
            $T^{⋆} = T^{⋄}$
        fin si
    fin pour
fin pour

Nous avons effectué des tests de performance pour une descente stricte avec implantation des arbres sous forme statique sur le problème m1038.sg (297 taxa de 2021 résidus) pour les processeurs suivants :

  • un AMD Phenom II 1090T Hexa core @ 3.20GHz
  • un Intel i5 4570 Quad core @ 3.20GHz
  • un Intel i7 4790 Quad Core + HyperThreading @ 3.60GHz , soit 4 coeurs + 4 HT
  • un Intel Xeon E5 2670 v2 @ 2.50GHz Deca Core (doté de 10 coeurs) avec HyperThreading, soit potentiellement 20 threads (10C + 10HT). Le noeud du cluster sur lequel on a réalisé les calculs est un bi-socket, soit deux Xeon E5, ce qui fait potentiellement 40 threads (20 + 20).

m1038.sg results

temps d'exécution en fonction du nombre de threads

3.5.3.a  AMD Phenom II 1090T

L'AMD Phenom II 1090T est le moins puissant des processeurs étudiés, il dispose cependant de 6 coeurs.

 #threads   1   2   4   6 
 temps (s)   328   186   112   91 
 facteur   -   1.76   2.92   3.06 
 pourcentage   -   -43%    -66%   -72% 
Temps d'exécution (s) sur AMD Phenom II 1090T

3.5.3.b  Intel i5 4570

 #threads   1   2   4 
 temps (s)   174   98   58 
 facteur   -   1.77   3.00 
 pourcentage   -   -44%   -67% 
Temps d'exécution (s) sur i5 4570

3.5.3.c  Intel i7 4790

Sur l'Intel i7-4790, l'utilisation de 8 threads est à proscrire car 4 threads sont des threads HyperThreading et ont tendance à dégrader les performances, il faut donc se limiter à 4 threads

 #threads   1   2   4   8 
 temps (s)   160   90   52   74 
 facteur   -   1.77   3.07   2.16 
 pourcentage   -   -44%   -68%   -54% 
Temps d'exécution (s) sur i7 4790

3.5.3.d  Intel Xeon E5 2670

Sur le Xeon E5 2670, la situation est très complexe avec des fluctuations importantes entre deux exécutions consécutives si le nombre de threads est supérieur à 10. La répartition des coeurs est la suivante :

 id   Cpu   Core 
 0   0   0 
 1   1   0 
 2   0   1 
 3   1   1 
 ...       
 20   0 (HT)   0 
 21   1 (HT)   0 
 ...       
Identification des coeurs sur Xeon E5 2670

En d'autres termes, les threads se répartissent :

 #threads   1   2   4   8   10   12   14   16   18   20   40 
 temps (s)   207   120   67   41   36   36   39   43   47   49   60 
 facteur   -   1.72   3.08   5.04   5.75   -   -   -   -   -   - 
 pourcentage   -   - 42%   - 68%   - 80%   - 83%   - 83%   - 81%   - 79%   - 77%   - 76%   - 71% 
Temps d'exécution (s) sur Xeon E5 2670

Que peut on déduire des résultats ?

  • jusqu'à 10 (voire 12) threads on note une diminution du temps de calcul
  • que l'on soit sur un seul CPU ou sur deux CPU (le tout sans HT) on obtient les mêmes résultats :
    • 8 threads : 41 secondes avec taskset -c 0-7 ou taskset -c 0,2,4,6,8,10,12,14
    • 10 threads : 36 secondes avec taskset -c 0-9 ou taskset -c 0,2,4,6,8,10,12,14,16,18
  • si on utilise l'HyperThreading on dégrade les performances :
    • 8 threads : 1m08s avec taskset -c 0-3,20-23 (4C + 4HT sur le même CPU)
    • 10 threads : 1m06s avec taskset -c 0-4,20-24 (5C + 5HT sur le même CPU)

3.5.3.e  Loi de Amdhal

Si on suit la loi de Amdhal, $R = 1/((1-S) + S/N)$ le gain obtenu devrait être le suivant sachant que :

Loi de Amdhal, $R$ (Gain) en fonction de $N$ (nombre de threads) pour $S=0.92$
$N$ $R$ % $Δ$ %
21.85 4646
43.23 6923
85.13 8011
105.81 833
126.38 841
207.94 873

Notons que plus $S$ est proche de 1, plus $R$ est grand quand $N$ augmente, cf. l'image ci-dessous..

Amdhal's Law

#threadsS=0.5S=0.6S=0.7S=0.8S=0.9S=0.92S=0.95S=0.99
21.33 1.43 1.54 1.67 1.82 1.85 1.90 1.98
41.60 1.82 2.11 2.50 3.08 3.23 3.48 3.88
81.78 2.11 2.58 3.33 4.71 5.13 5.93 7.48
101.82 2.17 2.70 3.57 5.26 5.81 6.90 9.17
121.85 2.22 2.79 3.75 5.71 6.38 7.74 10.81
201.90 2.33 2.99 4.17 6.90 7.94 10.26 16.81
Tableau de valeurs pour la loi de Amdhal

Par exemple avec 99% du temps d'exécution parallélisable, on peut diviser le temps d'exécution par un facteur 16.81 lorsque l'on utilise 20 threads/processeurs.