L'objet de ce chapitre est de familiariser l'utilisateur avec les techniques d'optimisation concernant :

• alignement
• l'accès mémoire
• les unités x87, MMX et SSE

Vous trouverez beaucoup d'informations sur le site d'Agner Fog, notamment dans les sections Optimizing subroutines in assembly language: An optimization guide for x86 platforms et The microarchitecture of Intel, AMD and VIA CPUs: An optimization guide for assembly programmers and compiler makers.

Trouver la valeur minimum entre 2 valeurs non signées :

sub eax,  ebx     ; = a-b
sbb edx,  edx     ; = (b > a) ? 0xFFFFFFFF : 0
and edx,  eax     ; = (b > a) ? a-b : 0
add ebx,  edx     ; Result is in ebx

Trouver la valeur minimum entre 2 valeurs signées :

sub eax,  ebx     ; Will not work if overflow here
cdq               ; = (b > a) ? 0xFFFFFFFF : 0
and edx,  eax     ; = (b > a) ? a-b : 0
add ebx,  edx     ; Result is in ebx

if (a < 0) d = b; else d = c;

test eax, eax
mov   edx, ecx
cmovs edx, ebx    ; = (a < 0) ? b : c

Nous avons déjà évoqué dans le chapitre 3, les problèmes liés à l'alignement des données en mémoire. Nous savons qu'il est préférable d'aligner les données sur une adresse multiple de 2, 4, 8 ou 16.

Plusieurs possibilités d'alignement existent :

• utilisation des options en ligne de commande du compilateur pour aligner le code :

  gcc -falign-functions=n  -falign-jumps=n -falign-labels=n  -falign-loops=n
  

  où n représente une puissance de 2.
  
  
• utilisation des directives du compilateur pour aligner les données :
  • gcc / g++ :

    int j __attribute__((aligned(16));
    

  • icc / icpc / icl :

    int j _declspec(align(16));
    

• lors de l'allocation mémoire, les données sont alignées mais pas forcément sur la puissance de 2 désirée, notamment pour les données qui seront utilisées par les unités SSE. On peut réaliser un alignement à la main comme suit :

  #include <stdio.h>
  #include <stdlib.h>
  
  #define malloc_align(p,q,T,size,align)     \
    printf("allocate %d bytes\n", size); \
    p = (T *) malloc(size + align-1); \
    q = (T *) (((int)p+align-1) & ~(align-1)); \
    printf("unaligned @ %p\n", p); \
    printf("aligned   @ %p\n", q);
  
  int main() {
    int i;
    int *tab_unaligned, *tab;
  
    malloc_align(tab_unaligned, tab, int, 10*sizeof(int), 16);
    for (i=0; i<10; ++i) tab[i] = i;
    free(tab_unaligned);
  
    return 0;
  }
  

  le résultat en sortie de ce programme est (par exemple) :

  allocate 40 bytes
  unaligned @ 0x804a008
  aligned   @ 0x804a010
  

Cette technique que je qualifie de multi-allocation, consiste à éviter de perdre du temps lors de nombreuses allocation/déallocation de la mémoire. Elle repose sur les principes suivants :

• les objets sont alloués par blocs et non pas un à un
• les objets ne sont pas libérés mais placés dans une liste chaînées lorsqu'ils ne sont plus utilisés

Cette technique possède des avantages et des inconvénients :

• + lors de l'allocation, on réalise une allocation de N objets plutôt que N allocations de 1 objet
• + lors de la désallocation, on place l'élément au début de la liste, on ne perd pas de temps à faire un appel système
• - les objets sont toujours présents en mémoire

Voici deux exemples, l'un en C l'autre en C++ :

• exemple en C : multi_alloc.h et multi_alloc.c
• exemple en C++ : multi_alloc.tgz

Sur un Pentium 4 le cache de données possède une taille de 8 ko et est organisé sous la forme d'un cache 4-way associative. Chaque ligne de cache occupe 64 octets.

Une rapide analyse montre que les 6 bits les moins significatifs (bits 0 à 5) d'une adresse permettent d'adresser l'un des 64 octets de la ligne de cache. En outre :

les 5 bits suivants (bits 6 à 10) permettent de sélectionner l'une des 32 lignes parmi les 4 blocs (4-way) possibles.

En conclusion, si les bits 6 à 10 de deux adresses différentes sont égaux, alors il seront placés dans le même espace de cache. Etant donné qu'il y a 4 blocs (4-way), un bloc de cache permet de stocker 8 ko / 4 = 2 ko.

• application au produit de matrices

for (int i = 0; i < n; i++) {
   // (loop body)
}

Dont la traduction basique est :

    mov ecx, n    ; Load n
    xor eax, eax  ;i=0
for_begin:
    cmp eax, ecx  ;i= n
    ; (loop body) ; Loop body goes here
    add eax, 1    ; i++
    jmp for_begin ; Jump back
for_end:

On préfère utiliser add eax,1 plutôt que inc eax,1, car l'instruction inc a un problème pour accèder et mettre à jour une partie des flags, ce qui la rend moins efficace.

L'autre problème est le fait que l'on trouve 2 instructions de saut jge et jmp. On écrira donc :

    mov  edx, n    ; Load n
    test edx, edx ; Test n
    jz   for_end   ; Skip if n <= 0
    xor  ecx, ecx  ; i = 0
for_begin:
    ; (loop body) ; Loop body goes here
    add  ecx, 1    ; i++
    cmp  ecx, edx  ; i < n
    jl   for_begin  ; Loop back if i < n
for_end:

comment générer des constantes entières pour les registres SSE :

L'instruction LEA est très utile car elle permet de regrouper :

• un décalage (1,2 ou 3 bits)
• deux additions
• et un déplacement

Par exemple, la série d'instructions suivantes :

mov  eax,ecx
shl  eax,3
add  eax,ebx
add  eax,200

peut être remplacée par :

lea eax,[ebx+ecx*8+200]

Attention : les symboles [] ne signifient pas qu'il faut charger dans eax la donnée pointée par l'expression entre les crochets, mais calculer ebx+ecx*8+200 et mettre ce résultat dans eax.

De même :

imul eax,5

peut être remplacée par :

lea eax,[eax+eax*4]

Les instructions inc et dec ne modifient pas le drapeau de retenue (carry flag) mais sont susceptibles de modifier d'autres drapeaux, cela correspond à ne modifier qu'une partie du registre des flags et cela se traduit par un coût d'une micro-opération supplémentaire sur Pentium 4 et des délais d'attente sur Penitum II, III et M si une autre instruction tente de lire le registre de flags.

La division entière et en flottant est généralement très pénalisante, bien plus que la multiplication. Il existe cependant quelques astuces qui permettent de simplifier la division.

On sait que diviser X par 2^N consiste à décaler la représentation binaire de X de N rangs vers la droite :

mov  eax,X
shr  eax,N

Dans le cas où on divise par une constante, il est possible de remplacer la division par une multiplication. Pour calculer q = X / d pour une représentation de l'information sur w bits, on calcule la valeur :

on multiplie ensuite X par f et on décale à droite de r positions.

l'algorithme est le suivant :

• b = (nombre de bits significatifs de d) - 1
• r = w + b
• f = 2^r / d

1. si f est en entier alors d est une puissence de 2 :
   result = x shr b
2. sinon si la partie fractionnaire de f < 0.5 alors
   arrondir f à l'entier k (k < f) le plus proche
result = ((X+1) × f) shr r
3. sinon (f > 0.5) alors :
   arrondir f à l'entier k (k > f) le plus proche
result = (X × f) shr r

example : on désire diviser par 5 sur 32 bits :

• 5₁₀ = 101₂
• b = 3 - 1 = 2
• r = 32 + 2 = 34
• f = 2³⁴ / 5 = 3435973836.8 = CCCCCCCC,CCCCCCCCD₁₆
• k = CCCCCCCD₁₆

L'utilisation des options en ligne de commande des compilateurs permet parfois de simplifier l'optimisation. Par exemple, le compilateur C Intel (icc ou icpc pour le C++), offre différents types d'optimisations :

• optimisation interprocédurale
• optimisation guidée par le profilage
• optimisation des boucles (vectorisation)

sur les architectures de type IA-32, elle concerne :

• le passage d'arguments de sous-programmes dans des registres plutôt que dans la pile
• le déplacement d'invariants en dehors des boucles

ce type d'optimsation se déroule en 3 phases :

1. compiler en utilisant l'option -prof_gen
2. exécuter le programme, un fichier .dyn est alors généré qui contient des informations qui sont utilisées dans la phase suivante
3. recompiler avec -prof-use qui prend en compte les informations générées durant l'exécution

Le compilateur Intel est capable de vectoriser certaines boucles. Par exemple, le code suivant :

i=0;
while (i < n) {
  a[i] = b[i] + c[i];
  ++i;
}

sera vectorisé en :

i=0;
while (i < (n - (n%4)) {
  a[i:i+3] = b[i:i+3] + c[i:i+3];
  i+=4;
}
while (i < n) {
  a[i] = b[i] + c[i];
  ++i;
}

Soit le code suivant :

#define N 1024
int __attribute__ ((__aligned__(16))) a[N];
int __attribute__ ((__aligned__(16))) b[N];
int __attribute__ ((__aligned__(16))) c[N];

void sum_vec() {
  for (int i=0;i<N;++i) {
    a[i]=b[i]+c[i];
  }
}

Si on compile avec g++, comme suit :

g++ -S -masm=intel a.cpp -ftree-vectorize -ftree-vectorizer-verbose=4 -O2 -msse2

on obtient le résultat suivant :

a.cpp:11: note: LOOP VECTORIZED.
a.cpp:11: note: vectorized 1 loops in function.

        push    %ebp
        xor     %eax, %eax
        mov     %ebp, %esp
        .p2align 4,,7
.L2:
        movdqa  %xmm0, XMMWORD PTR c[%eax]
        paddd   %xmm0, XMMWORD PTR b[%eax]
        movdqa  XMMWORD PTR a[%eax], %xmm0
        add     %eax, 16
        cmp     %eax, 4096
        jne     .L2
        pop     %ebp
        ret

Il est également possible de déplier la boucle :

g++ -S -masm=intel a.cpp -ftree-vectorize -ftree-vectorizer-verbose=4 -O2 -msse2
 -funroll-loops --param max-unroll-times=4

        push    %ebp
        xor     %edx, %edx
        mov     %ebp, %esp

        .p2align 4,,7
.L2:
        lea     %eax, [%edx+16]
        movdqa  %xmm3, XMMWORD PTR c[%edx]
        movdqa  %xmm2, XMMWORD PTR c[%eax]
        paddd   %xmm3, XMMWORD PTR b[%edx]
        paddd   %xmm2, XMMWORD PTR b[%eax]
        lea     %ecx, [%edx+32]
        movdqa  XMMWORD PTR a[%edx], %xmm3
        movdqa  %xmm1, XMMWORD PTR c[%ecx]
        paddd   %xmm1, XMMWORD PTR b[%ecx]
        movdqa  XMMWORD PTR a[%eax], %xmm2
        lea     %eax, [%edx+48]
        add     %edx, 64
        movdqa  %xmm0, XMMWORD PTR c[%eax]
        paddd   %xmm0, XMMWORD PTR b[%eax]
        cmp     %edx, 4096
        movdqa  XMMWORD PTR a[%ecx], %xmm1
        movdqa  XMMWORD PTR a[%eax], %xmm0
        jne     .L2
        pop     %ebp
        ret

mais, le code suivant ne sera pas vectorisé car le compilateur gcc ne peut déterminer si les variables a, b, c sont des tableaux alignés et qu'elle est leur taille. Par contre icc est capable de vectoriser cette boucle.

#define N 1024

void sum_vec(int *a, int *b, int *c) {
  for (int i=0;i<N;++i) {
    a[i]=b[i]+c[i];
  }
}

Avec le compilateur, icc, on obtient :

icc -O3 -S a.c -fsource-asm -use-asm -msse2 -axT
a.c(12): (col. 3) remark: LOOP WAS VECTORIZED.

sum_vec:
..B1.1:                         # Preds ..B1.0
;;; void sum_vec() {
;;;   for (int i=0;i<N;++i) {
        xorl      %eax, %eax                                    #12.3
                                # LOE eax ebx ebp esi edi
..B1.2:                         # Preds ..B1.2 ..B1.1
;;;     a[i]=b[i]+c[i];
        movdqa    b(%eax), %xmm0                                #13.10
        paddd     c(%eax), %xmm0                                #13.15
        movdqa    16+b(%eax), %xmm1                             #13.10
        paddd     16+c(%eax), %xmm1                             #13.15
        movdqa    32+b(%eax), %xmm2                             #13.10
        paddd     32+c(%eax), %xmm2                             #13.15
        movdqa    48+b(%eax), %xmm3                             #13.10
        paddd     48+c(%eax), %xmm3                             #13.15
        movdqa    64+b(%eax), %xmm4                             #13.10
        paddd     64+c(%eax), %xmm4                             #13.15
        movdqa    80+b(%eax), %xmm5                             #13.10
        paddd     80+c(%eax), %xmm5                             #13.15
        movdqa    96+b(%eax), %xmm6                             #13.10
        paddd     96+c(%eax), %xmm6                             #13.15
        movdqa    112+b(%eax), %xmm7                            #13.10
        paddd     112+c(%eax), %xmm7                            #13.15
        movdqa    %xmm0, a(%eax)                                #13.5
        movdqa    %xmm1, 16+a(%eax)                             #13.5
        movdqa    %xmm2, 32+a(%eax)                             #13.5
        movdqa    %xmm3, 48+a(%eax)                             #13.5
        movdqa    %xmm4, 64+a(%eax)                             #13.5
        movdqa    %xmm5, 80+a(%eax)                             #13.5
        movdqa    %xmm6, 96+a(%eax)                             #13.5
        movdqa    %xmm7, 112+a(%eax)                            #13.5
        addl      $128, %eax                                    #12.3
        cmpl      $4096, %eax                                   #12.3
        jb        ..B1.2        # Prob 99%                      #12.3
                                # LOE eax ebx ebp esi edi
..B1.3:                         # Preds ..B1.2

;;;   }
;;; } 

        ret       

Les intrinsics sont des fonctions codées en assembleur qui peuvent être insérées dans un programme C ou C++ et qui permettent de faire appel aux instructions SSE. Sous gcc elles peuvent être utilisées en faisant appel à <xmmintrin.h>.

Les registres SSE ne sont pas directement spécifiables, on utilise donc de nouveaux types de données :

• __m128 : représentation sous forme de 4 float
• __m128d : représentation sous forme de 2 double
• __m128i : représentation sous forme entière (16 bytes ou 8 words ou 4 double words ou 2 quad words)

La plupart des noms de fonctions se table sur le format _mm_<intrin_op>_<suffix> :

• <intrin_op> : représente l'opération à réaliser
• <suffix> : définit le type de donnée sur lequel on applique l'opération :
  • s : flottant simple précision
  • d : flottant double précision
  • i128 : entier sur 128 bits signé
  • i64 : entier sur 64 bits signé
  • u64 : entier sur 64 bits non signé
  • i32 : entier sur 32 bits signé
  • u32 : entier sur 32 bits non signé
  • i16 : entier sur 16 bits signé
  • u16 : entier sur 16 bits non signé
  • i8 : entier sur 8 bits signé
  • u8 : entier sur 8 bits non signé

exemple : intrinsics_1.c

guide de référence Intel des Intrinsics

On désire vectoriser le code C suivant :

void xor(char *t0, char *t1, char*t2) {
   for (int i=0;i<N;++i) {
     t2[i] = t0[i] ^ t1[i]; // partie à vectoriser
  }
}

la traduction assembleur vetorisée qui consiste à réaliser un XOR entre les deux vecteurs t0 et t1 et stocker le résultat dans le troisième vecteur t2 est :

movdqa   xmm0,[esi] ; avec esi stocke l'adresse de t0
movdqa   xmm1,[edi] ; avec edi stocke l'adresse de t1
pxor     xmm0,xmm1
movdqa   [ebx],xmm0 ; avec ebx stocke l'adresse de t2

Son équivalent en intrinsics serait :

__m128i __attribute__((__aligned__(16))) x0, x1;

x0 = _mm_load_si128((const __m128i *) t0);
x1 = _mm_load_si128((const __m128i *) t1);
x0 = _mm_xor_si128(x0,x1);
_mm_store_si128((__m128i *) t2,x0);

qui sera traduit en :

movdqa  xmm0,[t0]
movdqa  [x0],xmm0

movdqa  xmm0,[t1]
movdqa  [x1],xmm0

movdqa  xmm0,[x0]
movdqa  xmm1,[x1]
pxor    xmm0,xmm1
movdqa  [x0],xmm0

mov     xmm0,[x0]
mov     [t2],xmm0

ce qui est sans aucun doute moins efficace, mais il est intéressant de pouvoir utiliser les intrinsics soit lorsqu'on ne dispose pas de programme assembleur ou si on désire assurer une certaine compatibilité entre architectures différentes.

Un aspect souvent négligé en programmation est la mise au point de tests afin de vérifier les fuites mémoires (memory leaks) ou les accès en déhors de zones mémoires autorisées à être accèdées.

Voici un morceau de code qui peut mener à un arrêt brutal de votre programme (la fameuse segmentation fault sous Unix) :

// a.cpp
int main() {
int *tab = new int[10];
for (i=0; i!=11; ++i) tab[i]=0;
return 0;
}

On accède ici à tab[10], ce qui n'est normalement pas autorisé. Votre programme peut alors posséder deux comportements différents :

• malgrè l'erreur, le programme fonctionne correctement, car la donnée que vous modifiez n'a pas d'influence sur d'autres programmes
• le programme produit une erreur de segmentation car vous modifiez une donnée en mémoire qui est vitale

Pour déterminer les erreurs de ce type vous pouvez utiliser ElectricFence. ElectricFence est une librairie qui intercepte les allocations mémoires et détermine quand vous procédez à un accès en dehors de zones allouées. Pour utiliser ElectricFence vous pouvez, soit réaliser l'édition de lien avec la version statique de la librairie :

> g++ -o a.exe a.cpp -lefence -ggdb
> gdb a.exe
(gdb) run
Starting program: /home/richer/tmp/a.exe 
[Thread debugging using libthread_db enabled]
[New Thread -1211124016 (LWP 7616)]

  Electric Fence 2.1 Copyright (C) 1987-1998 Bruce Perens.

Program received signal SIGSEGV, Segmentation fault.
[Switching to Thread -1211124016 (LWP 7616)]
0x08048566 in main () at a.cpp:6
4         for (int i=0;i!=11;++i) tab[i]=0;
(gdb) print i
$1 = 10

ou utiliser la librairie dynamique à l'intérieur d'un débogueur, ddd par exemple. Avant de lancer l'exécution du programme, tapez la commande suivante à l'intérieur du debogueur :

set environment LD_PRELOAD=/usr/lib/libefence.so.0.0

Une fuite mémoire correspond à un espace alloué dynamiquement que l'on a oublié de libéré par la suite. Cela peu amener votre programme a utiliser énormément de mémoire et réduite ainsi ses performaces.

Il existe cependant un moyen très simple de détecter les fuites mémoires potentielles en utilisant valgrind :

valgrind --leak-check=full ./a.exe
==7862== Memcheck, a memory error detector.
==7862== Copyright (C) 2002-2006, and GNU GPL'd, by Julian Seward et al.
==7862== Using LibVEX rev 1658, a library for dynamic binary translation.
==7862== Copyright (C) 2004-2006, and GNU GPL'd, by OpenWorks LLP.
==7862== Using valgrind-3.2.1-Debian, a dynamic binary instrumentation framework.
==7862== Copyright (C) 2000-2006, and GNU GPL'd, by Julian Seward et al.
==7862== For more details, rerun with: -v
==7862== 
==7862== 
==7862== ERROR SUMMARY: 0 errors from 0 contexts (suppressed: 23 from 1)
==7862== malloc/free: in use at exit: 40 bytes in 1 blocks.
==7862== malloc/free: 1 allocs, 0 frees, 40 bytes allocated.
==7862== For counts of detected errors, rerun with: -v
==7862== searching for pointers to 1 not-freed blocks.
==7862== checked 134,252 bytes.
==7862== 
==7862== 40 bytes in 1 blocks are definitely lost in loss record 1 of 1
==7862==    at 0x4021A55: operator new[](unsigned) (vg_replace_malloc.c:195)
==7862==    by 0x8048550: main (a.cpp:4)
==7862== 
==7862== LEAK SUMMARY:
==7862==    definitely lost: 40 bytes in 1 blocks.
==7862==      possibly lost: 0 bytes in 0 blocks.
==7862==    still reachable: 0 bytes in 0 blocks.
==7862==         suppressed: 0 bytes in 0 blocks.
==7862== Reachable blocks (those to which a pointer was found) are not shown.
==7862== To see them, rerun with: --show-reachable=yes

Un grand nombre d'algorithmes ont été étudiés et améliorés de manière à obtenir des temps d'exécution minimaux.

• BLAS (Basic Linear Algebra) est une bibliothèque de sous-programmes liés au calcul vectoriel et matriciel.
• LAPACK (Linear Algebra PACKage) écrit en Fortran 77, concerne la résolution de systèmes d'équations linéaires, recherche de valeurs propres, factorisation de matrices,... LAPACK utilise une partie des sous-programmes de BLAS

Différentes implémentations de BLAS/LAPACK existent et dépendent du processeur utilisé :

• AMD : ACML
• HP : MLIB
• Intel : MKL
• SUN : SPL

Voici un petit exemple de produit matriciel qui utilise le sous-programme SGEMM (Single floating-point GEneral Matrix matrix Multiplication) de la librairie Intel MKL :

int size; // size of square matrix
float *A; // first square matrix
float *B; // second square matrix
float *C; // result

mkl_set_num_threads(2);

// compute C = alpha * A * B + beta * C

cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
	    size, size, size,
	    1.0, /* alpha */
	    A, size,
	    B, size,
	    0.0, /* beta */
	    C, size);

On compile le programme avec les librairies suivantes :

ALL=C icpc -o main.exe obj/*.o -O3 -axT -msse2 
-Wl,--start-group 
/opt/intel/mkl/10.0.011/lib/32/libmkl_core.a
/opt/intel/mkl/10.0.011/lib/32/libmkl_intel.a 
/opt/intel/mkl/10.0.011/lib/32/libmkl_blacs.a 
/opt/intel/mkl/10.0.011/lib/32/libmkl_lapack.a 
/opt/intel/mkl/10.0.011/lib/32/libmkl_intel_thread.a
/opt/intel/mkl/10.0.011/lib/32/libmkl_scalapack.a 
-Wl,--end-group -lguide -lpthread

Voici quelques résultats :

On remarquera que :

• bien que le Core 2 Duo E8400 fonctionne à une fréquence de 3 Ghz, il obtient des temps de calcul plus importants que le Core 2 Quad Q6600 qui ne fonctionne qu'à 2,4 Ghz. Cela est du principalement à la taille du cache L2 qui est de 8 Mo sur le Quad Core alors qu'il n'est que de 6 Mo sur le Dual Core
• la fréquence du FSB ne semble pas influencer les résultats puisque le Q6600 et le Q9300 obtienent sensiblement les mêmes résultats

Le GPU computing est destiné à accélérer les algorithmes massivement parallèles en profitant des nombreuses unités d'exécution présentes au seun d'un GPU. Ces algorithmes doivent être conçus de manière à se décompose en une multitude de petits threads qui seront exécutés en parallèle, par groupes, sur le GPU... Une façon de penser et de programmer totalement différente donc qui a demandé à NVidia de concevoir un langage adapté : C pour CUDA. Celui-ci a été écrit par Ian Buck, qui était à l'origine de Brook GPU, un langage destiné à exploiter le pipeline de rendu 3D, et non directement le coeur de calcul du GPU, ce qui le rendait plus rigide et moins performant.
(Citation: Hardware.fr, Damien Triolet)

L'idée de la programmation GPU repose sur l'utilisation des circuits GPU (Graphical Processing Unit) des cartes graphiques. En effet, ces circuits se révèlent très performants pour réaliser du calcul vectoriel ou matriciel.

Par exemple, dans le cas de la multiplication matricielle, la société GPU-Tech annonce en novembre 2007 des temps 10 foix supérieurs à un CPU pour des flottants en simple précision. Sur un AMD Athlon 64 3800+, on obtient les résultats suivants :

Image du site GPU-Tech

La difficulté de la programmation GPU est de faire cohabiter un CPU avec GPU.

quelques transparents :

• Séminaire IDRIS
• Getting started with CUDA
• Advanced CUDA Training

Le projet Folding@Home (2000) vise à étudier le repliement des protéines. Lorsque des protéines ne se replient pas correctement ou ne sont pas correctement formées, elles peuvent engendrer des maladies (Alzheimer, Encéphalopathie spongiforme bovine (Mad Cow Disease), Parkinson, Myopathie, ...).

Le problème lié au repliement des protéines est qu'il demande énormément de temps de calcul. L'idée du projet Folding@Home a été de demander à des particuliers de télécharger sur leur ordinateur un programme qui s'exécute lorsque la machine est en veille, et effectue des calculs. On dispose ainsi qu'une importante source de calculateurs.

Historique des améliorations du programme de calcul :

• 2000 : Tinker engine
• 2003 : Gromacs, utilisation de code assembleur (à la main) + instructions SSE (amélioration d'un facteur de 20 à 30 par rapport à Tinker)
• 2006 : utilisation de la carte graphique ATI X1900 (nouvelle amélioration d'un facteur de 20 à 30)

CUDA est un outil de développement pour les GPU supportant la technologie CUDA (Compute Unified Device Architecture) disponible sous Linux et Windows. CUDA comporte un certain nombre d'algorithmes optimisés, parmi lesquels on trouve :

• multiplication et transposition de matrices
• transformation de Fourier
• algorithmes de traitement d'images

Pour utiliser CUDA il suffit de disposer d'une carte graphique basée sur le chipset GeForce 8800. Il existe également des cartes graphiques dédiées au calcul comme la C870 dotée de 1,5 Mo de mémoire et d'une puissance de 500 Gflops (Tesla).

Caractéristiques GPU (Source 59Hardware)

Au mois de septembre 2009, NVidia présente sa future architecture Fermi (GT300, GF100), l'ancienne architecture était le Tesla (GT200).

• 3 Milliards de transistors sur une puce gravée en 40 nm
• 16 multiprocesseurs (512 cores organisés par groupes de 32)
• jusqu'à 6 Go de mémoire avec la GDDR5
• mémoire cache de niveau 1 et 2
• support langage C / C++

Le profiling et le tuning de programmes sont deux aspects complémentaires de l'optimisation :

• le profilage consiste à déterminer les hot spots, c'est à dire les endroits du code qui consomment le plus de ressources
• le réglage (tuning) consiste à modifier les hot spots de manière à les rendre plus cool (=frais, donc moins chauds)

Dans la suite de cette section, nous étudions le comportement du programme suivant : prog.c

On utilse les options de compilation classiques avec gprof :

> gcc -o prog.exe prog.c -pg
> ./prog.exe
> gprof prog.exe >prog_gprof.txt

Le résultat en sortie peut être visualisé ici.

Le flat profile est le suivant :

Flat profile:

Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total           
 time   seconds   seconds    calls   s/call   s/call  name    
 99.55     21.94    21.94  2000000     0.00     0.00  proc3
  0.27     22.00     0.06        2     0.03     0.03  fill_array
  0.09     22.02     0.02     2000     0.00     0.01  proc21
  0.09     22.04     0.02     2000     0.00     0.01  proc22
  0.00     22.04     0.00        2     0.00     0.00  allocate_array
  0.00     22.04     0.00        1     0.00    21.98  proc1

On peut donc lire que :

• le sous-programme proc3 a été appelé 2.000.000 de fois et a consommé 99.55% du temps CPU (soit 21.94s sur les 22.04s du temps d'exécution total), c'est donc cette fonction qu'il faudra optimiser
• proc21 a été appelé 2000 fois, ainsi que proc22 : bien que proc21 et proc22 appellent proc3, ici n'est reporté que le temps d'exécution hors appels aux autres sous-programmes.

Le second type d'information reporté par gprof concerne la granularité, c'est à dire le temps total passé dans chaque sous-programme ainsi que le nombre d'appels aux autres sous-programmes :

index % time    self  children    called     name
                                                 
[1]    100.0    0.00   22.04                 main [1]
                0.00   21.98       1/1           proc1 [2]
                0.06    0.00       2/2           fill_array [6]
                0.00    0.00       2/2           allocate_array [7]
-----------------------------------------------
                0.00   21.98       1/1           main [1]
[2]     99.7    0.00   21.98       1         proc1 [2]
                0.02   10.97    2000/2000        proc21 [4]
                0.02   10.97    2000/2000        proc22 [5]
-----------------------------------------------
               10.97    0.00 1000000/2000000     proc21 [4]
               10.97    0.00 1000000/2000000     proc22 [5]
[3]     99.5   21.94    0.00 2000000         proc3 [3]
-----------------------------------------------
                0.02   10.97    2000/2000        proc1 [2]
[4]     49.9    0.02   10.97    2000         proc21 [4]
               10.97    0.00 1000000/2000000     proc3 [3]
-----------------------------------------------
                0.02   10.97    2000/2000        proc1 [2]
[5]     49.9    0.02   10.97    2000         proc22 [5]
               10.97    0.00 1000000/2000000     proc3 [3]
-----------------------------------------------
                0.06    0.00       2/2           main [1]
[6]      0.3    0.06    0.00       2         fill_array [6]
-----------------------------------------------
                0.00    0.00       2/2           main [1]
[7]      0.0    0.00    0.00       2         allocate_array [7]
-----------------------------------------------

On voit donc que :

• [1] la fonction main consomme 100% du temps d'exécution, ce qui est normal puisque c'est elle qui se trouve au sommet de l'arbre d'appel des sous-programmes.
• [2] c'est ensuite au tour de proc1 d'appeler :
  • proc21 : 2000 fois, soit 10.97s
  • proc22 : 2000 fois, soit 10.97s
• [3] proc3 est appelé 2.000.000 de fois :
  • 1.000.000 de fois par proc21
  • 1.000.000 de fois par proc22

VTune est un outil complet de profiling et de tuning qui permet de déterminer les hot spots et d'éclairer sur leur comportement : instruction peu performante, accès mémoire non aligné, nombre important de défauts de cache dus à des données de trop grande taille.

• Call Graph Wizard
• Sampling Wizard (% utilisation par sous-programme)
• Sampling Wizard (proc3, code source)
• Sampling Wizard (proc3, code source et assembleur)

Pour rappel, cette fonction compilée avec des options de compilation différentes donnait des résultats de 8s ou 23s.

  push   ebp
  xor    ecx,ecx
  mov    ebp,esp
  mov    edx,DWORD PTR [ebp+8]
L1:
  movzx  eax,dl
  movsx  eax,BYTE PTR [eax+0x8048bd0]
  add    ecx,eax
  shr    edx,0x8
  test   edx,edx ; sans cette instruction 23s, avec 8s
  jne    L1
  pop    ebp
  mov    eax,ecx
  ret    

Un rapide analyse sous VTune, nous donne les résultats suivants pour les compteurs choisis :

• RAT_STALL.FLAGS : compte le nombre de cycles durant lesquels l'exécution est suspendue en raison d'un accès partiel au registre des FLAGS