Previous Next Contents

Appendix A    Remarques sur la validation de programmes

Sans vouloir (ni pouvoir) être exhaustif, voici quelques indications sur les aspects auxquels il faut faire attention quand on fait du calcul scientifique parallèle. Le premier aspect est celui de s'assurer que la coordination des processus est correcte. Pour montrer la difficulté de la chose, en particulier en mémoire partagée, on propose d'étudier quelques algorithmes classiques (voir aussi [Ray84]) d'exclusion mutuelle. A cette occasion, on utilise (sans trop formaliser, on pourra consulter [Win93] pour plus de détails) quelques notions de preuves à la Hoare. Le deuxième aspect est celui de la validité du calcul numérique proprement dit. On n'a pas insisté sur les aspects les plus classiques de stabilité des schémas etc. (voir par exemple [Bam98]). Ici on a fait le choix d'expliquer comment sont représenter les nombres flottants (qui ne sont pas les réèls mathématiques) en mémoire afin d'attirer l'attention sur les dangers de la programmation numérique.

A.1  Validation de la coordination

A.1.1  En mémoire partagée

Algorithmes d'exclusion mutuelle

Essayons de programmer en pseudo-code une première approximation de programme permettant d'assurer l'exclusion mutuelle. Intéressons-nous d'abord à un premier cas dans lequel on a seulement deux processeurs P0 et P1. L'idée est d'assurer que l'on ne puisse jamais être en même temps dans certaine partie du code (appelée section critique) pour P0 et P1.

Le premier code auquel on puisse penser pour Pi est (i=0,1, j=(i+1) mod 2):

VAR turn: 0..1;

LOOP
WAIT UNTIL (turn=i);
   {Section Critique}
turn:=j;
   {Section Non-Critique}
END LOOP
Mais dans ce cas les processus ne peuvent entrer en section critique qu'à tour de rôle, ce qui est très limitatif (en particulier si un des processeurs ou un des processus s'arrête sur une erreur pendant sa section critique...).

Pour résoudre cette difficulté il faut s'arranger pour connaître l'état de chaque processus afin de prendre une décision. Un essai de solution serait:

VAR flag: ARRAY[0..1] OF BOOL;

WAIT UNTIL (flag(j)=FALSE);
flag(i):=TRUE;
   {Section Critique}
flag(i):=FALSE;
Mais hélas maintenant cela n'assure plus l'exclusion mutuelle!

Pour s'en convaincre il suffit de considérer l'exécution suivante:

Pi Pj flag(i) flag(j)
flag(j)=FALSE ?: OUI flag(i)=FALSE ?: OUI FALSE FALSE
flag(i):=TRUE flag(j):=TRUE TRUE TRUE
Section Critique Section Critique - -

On propose alors le code suivant, toujours plus compliqué...

flag(i):=TRUE;
WHILE (flag(j)=TRUE) DO flag(i):=FALSE;
                        WAIT UNTIL NOT(flag(j));
                        flag(i):=TRUE;
END WHILE
   {Section Critique}
flag(i):=FALSE;
Il réalise bien l'exclusion mutuelle mais peut causer un interblocage si les deux processus s'exécutent de façon purement synchrone. A ce moment là, après la première affectation, les deux tests de la boucle WHILE sont vrais, d'où les deux flag sont mis à FALSE simultanément puis à TRUE et ainsi de suite. Les programmes de P0 et de P1 bouclent avant la portion Section Critique.

Une solution est d'instaurer une règle de politesse entre les processeurs. C'est l'algorithme classique décrit dans la section suivante.

L'algorithme de Dekker (1965)

VAR flag: ARRAY[0..1] OF BOOL;
    turn: 0..1;

flag(i):=TRUE;
WHILE (flag(j)=TRUE) DO IF turn=j 
                        THEN flag(i):=FALSE;
                             WAIT UNTIL NOT(turn=j);
                             flag(i):=TRUE;
                        END IF
END WHILE
   {Section Critique}
turn:=j;
flag(i):=FALSE;
Il réalise l'exclusion mutuelle sans blocage mutuel.

L'``algorithme'' de Hyman (1966)

Considérons le code plus simple suivant, proposé par Hyman en 1966:

VAR flag: ARRAY[0..1] OF BOOL;
    turn: 0..1;

flag(i):=TRUE;
WHILE NOT(turn=i) DO
                     WAIT UNTIL NOT(flag(j));
                     turn:=i;
END WHILE
   {Section Critique}
flag(i):=FALSE;
Il est malheureusement faux car, si turn est à 0 et P1 positionne flag(1) à TRUE puis trouve flag(0) à FALSE, alors P0 met flag(0) à TRUE, trouve turn égal à 0 et rentre en section critique. P1 affecte alors 1 à turn et rentre également en section critique!

En conclusion, il n'est pas si facile de programmer avec des processus communiquant par mémoire partagée (cf. JAVA par exemple).

L'algorithme de Peterson (1981)

Une vraie amélioration de l'algorithme de Dekker est l'algorithme de Peterson proposé en 1981.

Contentons nous d'abord du cas simple où l'on dispose uniquement de deux processeurs P0 et P1. Le code de Pi est (i=0,1, j=(i+1) mod 2):

VAR flag: ARRAY[0..1] OF BOOL;
    turn: 0..1;

flag(i):=TRUE;
turn:=i;
while (flag(j)=TRUE et non(turn=j)) do
        skip;
od;
   {Critical Section}
flag(i):=FALSE;
On a bien alors exclusion mutuelle sans interblocage et avec équité (pas de famine).

Donnons ici quelques idées de preuve de ces assertions.

Introduisons deux variables auxiliaires after(1) et after(2) (qui servent à indiquer si le contrôle est après turn=i dans Pi en rajoutant les instructions flag(1)=FAL SE et flag(2)=FALSE avant la mise en parallèle des deux processus et after(1)=TRUE, after(2)=TRUE après respectivement turn=1 et turn=2.

En utilisant l'abbréviation I= [turn=1 ou turn=2] on commente le code avec les assertions logiques suivantes:

{ non flag(1) }
flag(1)=TRUE; after(1)=FALSE;
{ flag(1) et non after(1) }
turn=1; after(1)=TRUE;
{ inv: flag(1) et after(1) et I }
while (flag(2) et non(turn=2)) do
   { flag(1) et after(1) et I }
   skip;
od;
{ flag(1) et after(1) et (flag(2) et (non(after(2)) ou turn=2)) }
[ Section Critique CS1]
{ flag(1) }
flag(1)=FALSE;
De même pour le deuxième processus,

{ non flag(2) }
flag(2)=TRUE; after(2)=FALSE;
{ flag(2) et non after(2) }
turn=2; after(2)=TRUE;
{ inv: flag(2) et after(2) et I }
while (flag(1) et non(turn=1)) do
   { flag(2) et after(2) et I }
   skip;
od;
{ flag(2) et after(2) et (flag(1) et (non(after(1)) ou turn=1)) }
[ Section Critique CS2]
{ flag(2) }
flag(2)=FALSE;
Il est relativement aisé de se convaincre ces assertions forment bien un shéma de preuve correct de chacun des processus pris séquentiellement. Seule subsiste la preuve de non-interférence. En fait, seul, 
pre(CS1)=flag(1) et after(1) et (flag(2) et (non after(2) ou (turn=2)))
contient des références à des objets manipulés par l'autre processus, et ce, seulement en les lignes flag(2)=TRUE; after(2)=FALSE; et turn=2; after(2)=TRUE; . Or, 
{ pre(CS1) } flag(2)=TRUE; after(2)=FALSE; { pre(CS1) }
{ pre(CS1) } turn=2; after(2)=TRUE; { pre(CS1) }
et de même par symétrie. De plus, on peut voir que 
pre(CS1) et pre(CS2) implique (turn=1 et turn=2)
Ce qui implique que l'on a toujours 
non(pre(CS1) et pre(CS2))
prouvant ainsi l'exclusion mutuelle.

Généralisation de l'algorithme de Peterson à n processeurs

Le code suivant, codé sur chacun des n processeurs (P0 à Pn-1) réalise bien l'exclusion mutuelle:

VAR flag: ARRAY[0..n-1] OF [-1..n-2]
    turn: ARRAY[0..n-2] OF [0..n-1]

FOR (j=0..n-2) DO
                  flag(i):=j;
                  turn(j):=i;
                  WAIT UNTIL ((all k<>i, flag(k)<j) OR (turn(j)<>i));
END FOR
   { Section Critique }
flag(i):=-1;
En fait c'est la solution pour 2 processeurs itérée n-1 fois. flag(i)=-1 signifie que Pi ne s'engage pas en section critique. turn permet de gèrer les conflits dans un couple de processus. Cet algorithme réalise bien l'exclusion mutuelle sans interblocage et avec équité.

Dans les améliorations possibles, il y a en particulier l'algorithme de Burns (1981). C'est une solution entièment symétrique qui minimise le nombre de variables utilisées ainsi que les valeurs qu'elles peuvent prendre en cours de programme.

A.1.2  Algorithmes distribués

Supposons que chacun des p processeurs Pi (i=0,...,p-1) est muni d'une variable booléenne locale demandeuri indiquant s'il désire rentrer ou non en section critique. Tous les Pi essaient d'exécuter leur section critique SCi si et seulement si demandeuri est vrai. Si demandeuri est faux alors Pi exécute HSCi (``Hors Section Critique''). Seul HSCi peut modifier la variable demandeuri. Il existe alors plusieurs méthodes distribuées permettant d'assurer la propriété de section critique.

Algorithme de LeLann

On organise de façon logique les p processeurs en anneau. On considère un unique jeton tournant sur cet anneau. Le processeur qui possède le jeton est seul autorisé à exécuter sa section critique.

Algorithme de Ricart et Agrawala

Lorsqu'un processeur désire entrer en section critique, il demande la permission à chacun des p-1 autres processeurs. Une fois qu'il a reçu toutes les permissions, il est autorisé à exécuter sa section critique.

Par exemple dans le cas où p vaut 2 et on définit P1 et P2 de la façon suivante:

Mais si les demandes et les accords de permissions sont des messages synchrones on a un problème potentiel d'interblocage. Inversement si les demandes et les accords de permissions sont des messages asynchrones il est possible que les deux processeurs rentrent en section critique au même moment!

Une façon de résoudre le problème est de supposer que l'on a une horloge globale (qui peut être par exemple réalisée à partir d'horloges locales). On estampille les demandes par l'heure de départ du message et par l'identité du demandeur, et on peut alors choisir de demander et d'accorder les permissions par des messages asynchrones.

A.2  Validation du calcul numérique

A.2.1  De l'imprécision en général

Un ordinateur calcule essentiellement faux: les entiers machines ne sont pas les entiers mathématiques (ils sont bornés), et les flottants ne sont qu'une représentation approchée finie des nombres réels. Par exemple tous les nombres transcendants (p, e etc.) sont approximés, de même que les rationnels ayant une expansion binaire infinie: attention,
1
10
=0.00011001100110011001100···
et est donc tronqué en machine (voir la section suivante pour plus de détails). Les opérations élémentaires ne vérifient pas les lois algébriques usuelles (associativité en particulier) parce-que des arrondis se glissent entre chaque opération. Enfin, le problème de l'imprécision ne touche pas seulement les quelques dernières décimales comme on pourrait le croire. Par exemple f(a,b)=333.75b6+a2(11a2b2-b6-121b4-2)+5.5b8+a/ 2b. En toutes précisions sur SUN: f(77617,33096)=1.172603..., en fait f(77617,33096)=-0.827396.... Le calcul numérique mal compris peut amener à des catastrophes...

A.2.2  La norme IEEE754

C'est la norme de calcul (1985) sur les nombres flottants, à laquelle vous aurez peut-être accès en C (en utilisant des librairies, telles ieee754.h sur SUN). Ce qu'il faut retenir c'est que toutes les représentations actuelles sur ordinateur convergent vers un tel standard. Une bonne partie du hardware (sauf quelques Cray) suivent la norme, mais les compilateurs ne sont pas encore à niveau (Microsoft...), en particulier pour la gestion des exceptions arithmétiques (``underflow'' etc.).

Représentation en mémoire
La norme définit les nombres standards, r = s*n*2k+1-N, avec s=-1,1, 1-2K < k < 2K, -2N < n < 2N normalisés pour que k soit tel que r=s*2k(1+f) avec f < 1, en plusieurs versions, Elle définit les nombres dénormalisés (pour gérer les ``underflow'' de façon graduelle), r=s*n*2k+1-N=s*2k(0+f) avec k=2-2K et 0 < n < 2N-1 c.a.d. 0<f<1.

Elle définit également +¥ et -¥ (remarquer que leurs inverses, +0 et -0 sont là également).

Enfin, elle définit les NaN ``Not a Number'' signés ou pas (par exemple 0 * ¥).

Quelques exemples
Pour un flottant simple, le maximum normalisé est de 3.402 82347*1038, le minimum normalisé positif est de 1.17549435*10-38, le maximum dénormalisé positif est de 1.17549421*10-38 et le minimum dénormalisé positif est de 1.40129846*10-45.

Autour de 1 l'erreur maximale (``unit in the last place'') est de 2-23 pour un flottant simple c'est à dire environ 1.19200928955*10-7.

Opérations
La norme impose que +, -, *, /, Ö soient calculés avec une imprécision ne pouvant dépasser la moitié de la distance entre deux flottants du type choisi autour du résultat exact (s'il n'y a pas d'``overflow''). Quand on a le choix entre deux arrondis possibles, on choisit celui qui a une mantisse paire.

On a également accès à des modes d'arrondi,

A.2.3  Que faut-il faire?

Dans le cadre du calcul scientifique se préoccuper des problèmes de précision de calcul est fondamental. Il faut se poser des questions tant au niveau de schéma purement mathématique de calcul qu'au niveau informatique.

D'un point de vue mathématique, si on a un problème mal conditionné (par exemple le calcul de l'inverse d'une matrice avec un tout petit déterminant), on aura toutes les raisons de craindre le pire numériquement.

D'un point de vue informatique, on peut avoir certaines catastrophes dues à une mauvaise compréhension de l'ordre de grandeur des quantités en jeu.

Dans tous les cas, après une nécessaire réflexion mathématique sur le schéma, il faudra penser à valider numériquement le programme écrit. On pourra par exemple, dans le cas de l'équation des ondes, comparer une solution explicite dans un cas simple à la résolution approchée (par le schéma discret choisi) de cette même solution. L'analyse mathématique du shéma doit donner une indication de l'écart entre les deux méthodes.

A.2.4  Un exemple

On considère le calcul itératif de la puissance nième du nombre d'or F=Ö5-1/ 2 par deux méthodes différentes,

Pour la méthode (1):

#include <stdio.h>
#include <math.h>

main()
{
  float x,y,z;
  int i;
  x=1;
  y=(sqrt(5)-1)/2;
  for (i=1;i<=20;i++)
    {
      z=x;
      x=y;
      y=z-y;
      printf("phi^%d=%f\n",i,x);
      }
}
Résultats:

phi^1=0.618034
phi^2=0.381966
phi^3=0.236068
phi^4=0.145898
phi^5=0.090170
phi^6=0.055728
phi^7=0.034442
phi^8=0.021286
phi^9=0.013156
phi^10=0.008130
phi^11=0.005026
phi^12=0.003103
phi^13=0.001923
phi^14=0.001180
phi^15=0.000743
phi^16=0.000437
phi^17=0.000306
phi^18=0.000131
phi^19=0.000176
phi^20=-0.000045
Visiblement fortement erroné...

En fait, à chaque tour de boucle les erreurs sur x, y et z sont données par la transformation affine,

æ
ç
ç
è
dx'
dy'
dz'
ö
÷
÷
ø
=
æ
ç
ç
è
0 1 0
1 -1 0
1 0 0
ö
÷
÷
ø		 æ
ç
ç
è
dx
dy
dz
ö
÷
÷
ø		 + æ
ç
ç
è
0
a
0
ö
÷
÷
ø

Avec au premier tour de boucle dy=e=2-23. On voit facilement que a £ 2-24 (en fait est de l'ordre de 2-23-i). D'autre part la transformation linéaire sur les erreurs a deux valeurs propres non nulles r1=-1+Ö 5/ 2 et r0=-1-Ö 5/ 2 < -1 !. D'où une erreur à la vingtième itération de l'ordre de e/ r1-r0 r120 + | r121 | -1/ r1-1 a (environ 4*10-4+2.4*10-3) supérieur à r020 (environ 6.6*10-5).

.5cm

Maintenant, au tour du programme suivant la méthode (2), d'ailleurs pas nécessairement plus lent:

#include <stdio.h>
#include <math.h>

main()
{
  float t,p;
  int i;
  t=1;
  p=(sqrt(5)-1)/2;
  for (i=1;i<=20;i++)
    {
      t=t*p;
      printf("phi^%d=%f\n",i,t);
      }
}
Résultats:

phi^1=0.618034
phi^2=0.381966
phi^3=0.236068
phi^4=0.145898
phi^5=0.090170
phi^6=0.055728
phi^7=0.034442
phi^8=0.021286
phi^9=0.013156
phi^10=0.008131
phi^11=0.005025
phi^12=0.003106
phi^13=0.001919
phi^14=0.001186
phi^15=0.000733
phi^16=0.000453
phi^17=0.000280
phi^18=0.000173
phi^19=0.000107
phi^20=0.000066
Il est facile de voir que l'erreur à la vingtième itération est de l'ordre de [ r0-2-23,r0+2-23 ]20 c'est à dire environ [ 6.61067063328*10-5,6.61072163724*10-5 ].

Dans le domaine du calcul scientifique, on retrouvera de tels systèmes dynamiques (pour des schémas explicites). Heureusement, la plupart du temps, la condition de stabilité des schémas assurera une non-dégradation des erreurs. Pour les schémas implicites, on aura à résoudre des systèmes linéaires. Là, selon les méthodes utilisées, on peut avoir des choix algorithmiques meilleurs que d'autres en termes de précision de calcul. Par exemple, on choisira plutôt le pivot le plus grand en norme dans la méthode d'élimination de Gauss.

Previous Next Contents