Document Actions
Conseils Pratiques
Mise en oeuvre de quelques optimisations algorithmiques
Introduction
Cette page contient quelques exemples de choses à éviter ou à mieux écrire .... et propose une solution.Elle essaye aussi de donner des conseils pour la portabilité des applications.
Pour tout problème ou question, n'hésitez pas à ouvrir un ticket d'assistance.
Type de données pour les flottants
Le type de données pour les nombres réels, appelés aussi flottants, peut réserver des surprises; en effet, pour un même type, on n'a pas forcément la même précision selon la machine utilisée car parfois la signification n'est pas la même.Tout d'abord, partons du principe que l'on souhaite travailler avec environ 15 décimales, soit avec le type DOUBLE PRECISION en général, sauf (liste non exhaustive) pour les machines Cray où ce n'est par la même chose ...
A titre d'exemple, considérons la représentation des données sur un Cray T3E, une SGI Origin 2000 et un IBM p690.
| Cray T3E | ||
| SGI Origin 2000 | ||
| IBM p690 |
Ainsi le même code avec le type DOUBLE PRECISION n'aura pas le même comportement; en effet ce choix a plusieurs conséquences :
- les données prennent deux fois plus de place en mémoire,
- une part importantes des fonctions mathématiques ne supporte pas le type,
- des travaux de re-écriture sont nécessaires.
Ainsi, on peut obtenir ce genre de messages à la compilation :
cf90-1110 f90: WARNING CODE, File = calcul.f, Line = 403, Column = 32
DOUBLE PRECISION is not supported on this platform. REAL will be used.
TL = 1.d0/(1.d0+TT)
Une alternative est de travailler avec un type, dépendant d'un paramètre qui est la seule valeur à modifier pour tout le code.
Par exemple en Fortran 90, les déclarations suivantes permettent de construire le paramètre :
MODULE numerics
INTEGER, PARAMETER :: sp = kind(1.0 ) ! nombre d'octets pour de la simple precision
INTEGER, PARAMETER :: dp = kind(1.0D0) ! nombre d'octets pour de la double precision
INTEGER, PARAMETER :: rp = dp
END MODULE numerics
On peut insérer ces lignes dans un module que l'on employe dans chaque sous programme par l'instruction use.
La déclaration des variables est alors la suivante :
use numerics
INTEGER :: i = 7
REAL(rp) :: epsas = 1E-5_rp
REAL(rp) :: srs6 = 6.0_rp
REAL(rp), DIMENSION(10) :: vxa
L'utilisation de constantes scalaires, d'entiers dans des expressions flottantes s'exprime alors :
vxa(i) = SQRT(24.0_rp) * epsas * srs6 * (2.0_rp * srs6 - 1.0_rp) / REAL(i, rp)
Les valeurs réelles sont toutes calculées selon la précision souhaitée et une seule modification est nécessaire.
De plus, la compilation se fait sans option concernant la précision des flottants, i.e. pas de -r8, -r4 ou autre, la déclaration impose le type.
Mesure du temps dans les programmes
Il existe plusieurs fonctions Fortran qui permettent de mesurer le temps écoulé.- timef :
Exemple :
REAL(8) ELAPSED, timef
ELAPSED = timef()
...
ELAPSED = timef()
Cette function s'applique aux programmes séquentiels. Pour les programmes parallèles comme OpenMP, elle renvoie le temps d'une thread et non pas le temps cumulé. Pour avoir le cumul, il faut employer la function global_timef.
- global_timef :
Exemple :
REAL(8) ELAPSED, global_timef
- mpi_wtime :
Exemple :
REAL(8) start_time, end_time, mpi_wtime, ELAPSED
start_time = mpi_wtime()
...
end_time = mpi_wtime()
ELAPSED = end_time - start_time
- system_clock et cpu_time :
Type de données MPI (ou SHMEM/PVM) pour les flottants
Dans le cas de données flottantes à envoyer ou recevoir d'autre processus, il est nécessaire de fournir le type des données au moyen des constantes MPI_REAL et MPI_DOUBLE_PRECISION. Dans le cas du portage du code d'un calculateur vers un autre, on peut rencontrer les mêmes difficultés que pour les déclarations des flottants. On peut construire une parade analogue à l'aide d'une variable entière ...INTEGER :: ITYPE_REAL
... variable initialisée au début du programme par tous les processus, comme par exemple par :
! Initialize MPI
CALL MPI_INIT(ierr)
!
! Get my process rank
CALL MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ierr)
!
! Get the number of processors
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nbproc, ierr)
!
! Float type
IF ( rp == dp ) THEN
ITYPE_REAL = MPI_DOUBLE_PRECISION
ELSE
ITYPE_REAL = MPI_REAL
END IF
Ainsi la valeur de rp ajuste le type de données
Pour les communications, on a alors :
CALL MPI_SEND ( x(1), 9, ITYPE_REAL, isend, istag, MPI_COMM_WORLD, ierr )
Générateur de nombres pseudo-aléatoires et MPI/OpenMP
Il est facile d'utiliser le générateur pseudo-aléatoire intrinsèque du Fortran 90 (RANDOM_NUMER) dans un code parallèle comme MPI ou OpenMP.Un générateur pseudo-aléatoire est une suite de nombres définie par récurrence à partir d'un ou plusieurs premiers termes, appelés germe(s).
Dans le cas d'une application parallèle, il est impératif que chaque processus (resp. thread) dispose de ses propres germes, différents de ceux de tous les autres processus (resp. threads).
La subroutine RANDOM_SEED permet d'effectuer ces initialisations pour le générateur RANDOM_NUMBER.
Voici la copie écran de la version MPI dans le cas de 4 processus :
GM-% ./a.out -procs 4
my_rank = 2 K = 2 AVANT iseed (1:K) = 1072693248 0
my_rank = 0 K = 2 AVANT iseed (1:K) = 1072693248 0
my_rank = 1 K = 2 AVANT iseed (1:K) = 1072693248 0
my_rank = 3 K = 2 AVANT iseed (1:K) = 1072693248 0
my_rank = 0 K = 2 APRES iseed (1:K) = 0 0
my_rank = 1 K = 2 APRES iseed (1:K) = 214538649 0
my_rank = 2 K = 2 APRES iseed (1:K) = 357564416 0
my_rank = 3 K = 2 APRES iseed (1:K) = 459725676 0
my_rank = 3 X = 0.9829E+00 0.8367E+00 0.4652E-01 0.8633E+00
my_rank = 0 X = 0.1565E-04 0.2631E+00 0.5112E+00 0.9173E+00
my_rank = 1 X = 0.5870E-01 0.4922E+00 0.7565E+00 0.1252E+00
my_rank = 2 X = 0.4312E+00 0.6624E+00 0.2684E+00 0.4868E+00
Le tableau suivant donne accès aux fichiers source des deux versions parallèles et indique la marche à suivre pour les tester :
| Type | Compilation | Exécution | Fichier source |
|---|---|---|---|
| MPI | mpxlf90 -O2 -qfixed=132 -qextname=flush rnd_mpi.f | ./a.out -procs 4 | rnd_mpi.f |
| OpenMP | xlf90_r -qsmp=omp -O2 -qfixed=132 -qextname=flush rnd_omp.f | setenv OMP_NUM_THREADS 4 ; ./a.out | rnd_omp.f |
Bibliothèques scientifiques
Les bibliothèques scientifiques sont des ensembles de sous-programmes testés, validés, optimisés. Utiliser des librairies scientifiques permet de se consacrer uniquement aux nouveaux développements et donc de ne pas devoir réinventer la roue à chaque fois. Elles ont plusieurs qualités qui les rendent très intéressantes :- elles sont (souvent) portables puisqu'indépendantes de l'architecture des machines;
- elles supportent différents types de donnée : réel ou complexe, simple ou double précision;
- elles prennent en compte différents types de stockage : matrice bande, symétrique, hermitienne.
D'autres bibliothèques conçues par les constructeurs présentent des avantages similaires mais ne sont pas portables puisque dépendantes de l'architecture, le cas typique est les transformées de Fourier rapides ou FFT. Elles font toutes la même chose (1D, 2D, 3D, complexe/complexe, complexe/réel ou réel/complexe, etc ...) mais les séquences d'appels changent d'un constructeur à l'autre (tableau de travail réel ou complexe, trigs à stocker de manière particulière, etc ...).
Dans ce cas, on peut procéder de la manière suivante pour essayer de limiter les modifications lors du portage.
On écrit une routine qui sert d'interface entre le programme et la bibliothèque. A cette routine, chacun des sous-programmes transmet les arguments suivants :
- le tableau avec les valeurs en entrée,
- le tableau pour les valeurs en sortie,
- le nombre de points dans chaque direction,
- le sens de la transformation.
Voir les exemples tirés de la page des bibliothèques scientifiques.
Remarque :
ces exemples reposent sur le fait que le tableau d'entrée n'est pas le même que le tableau de sortie (cas général), même si certaines fonctions de la librairie ESSL le supportent.
Enfin, la bibliothèque mathématique standard (fonctions trigonométriques, logarithmiques, etc ...) possède parfois une version optimisée, i.e. plus rapide avec une légère perte de précision. C'est le cas sous AIX avec la librairie MASS (Mathematical Acceleration SubSystem) qui propose deux instances : une instance pour un appel isolé à une fonction et une instance pour un appel sur un vecteur de données à traiter.
L'usage de cette librairie n'entraine aucune modification du code source, il suffit de faire l'édition des liens en ajoutant les deux options -lmass et -lmassv.
Recopiage des données
L'algorithme utilisé peut nécessiter la présence simultanée des valeurs correspondant à deux itérations successives.Dans l'écriture de l'algorithme qui suit, chaque itération demande la recopie du tableau des nouvelles valeurs obtenues dans le tableau des anciennes, (on suppose que les initialisations sont faites).
Exemple :
Version initiale ...
REAL(rp), DIMENSION(n:n) :: F
REAL(rp), DIMENSION(n:n) :: U_new, U_old
REAL(rp) :: h
INTEGER :: i, j
DO WHILE ( MAXVAL( ABS(U_new - Uold) ) > erreur ) U_old = U_new
DO j = 1, n
DO i = 1, nEND DOU_new(i,j) = 0.25_rp *(END DO
U_old(i-1,j) + U_old(i+1,j)
+ U_old(i ,j-1) + U_old(i ,j+1) )
- h * h * F(i,j)
END DO
On peut éviter cette recopie en utilisant un couple d'entiers qui joueront en quelque sorte le role de "pointeurs".
Version proposée ...
REAL(rp), DIMENSION(n:n) :: F
REAL(rp), DIMENSION(n:n,0:1) :: U
REAL(rp) :: h
INTEGER :: i, j
INTEGER :: iold = 0
INTEGER :: inew = 1
INTEGER :: itemp
DO WHILE ( MAXVAL( ABS(U(:,:,inew) - U(:,:,iold) ) ) > erreur )
itemp = iold
iold = inew
inew = itemp
DO j = 1, nEND DODO i = 1, nEND DOU(i,j,inew) = 0.25_rp *(END DOU(i-1,j ,iold) + U(i+1,j ,iold)
+ U(i ,j-1,iold) + U(i ,j+1,iold) )
- h * h * F(i,j)
On remarque que le volume en mémoire reste le même.
On peut aller encore un peu plus loin : dans cette boucle le compilateur n'arrive pas forcément à comprendre que iold et inew sont différents à chaque itération et donc qu'il n'y aura jamais de recouvrement entre ce qui est lu et ce qui est écrit pour un couple (i,j) donné. Ainsi il risque de mettre des "verrous", des recopies en mémoire (cohérence globale) des registres, etc ...
La version suivante, en allongeant un peu le code, permet de gagner encore en performance (10%), ce qui peut être très intéressant si cette boucle est au coeur du code :
DO WHILE ( MAXVAL( ABS(U(:,:,inew) - U(:,:,iold) ) ) > erreur )
itemp = iold
iold = inew
inew = itemp
IF ( inew == 0 ) THENEND DODO j = 1, nELSEDO i = 1, nEND DOU(i,j,0) = 0.25_rp *(END DOU(i-1,j ,1) + U(i+1,j ,1)
+ U(i ,j-1,1) + U(i ,j+1,1) )
- h * h * F(i,j)DO j = 1, nEND IFDO i = 1, nEND DOU(i,j,1) = 0.25_rp *(END DOU(i-1,j ,0) + U(i+1,j ,0)
+ U(i ,j-1,0) + U(i ,j+1,0) )
- h * h * F(i,j)
Remarque :
il est capital que le test IF ... THEN ... ELSE reste en dehors des
deux boucles imbriquées sur i et j.
Calcul de puissances de nombres flottants
Le calcul des puissances peut révéler quelques subtilités ...D'un point de vue mathématique, on a x**y = EXP( y * LN( x ) ). Le calcul d'une exponentielle et d'un logarithme est relativement élevé en terme de coût CPU. Dans le cas d'un exposant entier, il s'agit d'une élévation à la puissance et le calcul ne doit pas être fait avec ces deux fonctions mathématiques; il faut donc faire attention à ce que l'exposant soit bien vu comme une entité de type entier. De plus en ne faisant pas attention on peut facilement perdre du temps CPU .... et même de la précision dans le cas général.
Exemple :
REAL(rp), DIMENSION(n) :: a, b
DO i = 1, n
a(i) = b(i) ** 3.0_rp <-- mauvais car l'exposant est écrit comme un flottantEND DO
a(i) = b(i) ** 3 <-- meilleur car l'opération est faite en deux multiplications
De même lorsque l'on calcule les puissances successives, ....
REAL(rp), DIMENSION(n) :: a, b
REAL(rp) :: x
DO i = 1, n
a(i) = x**i <-- nécessite n*(n-1)/2 multiplicationsEND DO
b(1) = x
DO i = 1, n-1
b(i+1) = b(i) * x <-- nécessite n-1 multiplicationsEND DO
Remarque :
Toutefois d'autres considérations, liées à l'architecture de la machine, peuvent amener à privilégier un algorithme un peu plus cher en terme de nombre d'opérations, s'il en découle un gain significatif en terme de performances (vectorisation, parallélisation). Des comparaisons de performances permettent généralement de déterminer l'algorithme le mieux adapté.
Le phénomène s'aggrave lorsque l'on travaille avec des complexes puisque un complexe est formé de deux réels ...
De même, les puissances successives du i complexe (i^2 = -1) peuvent ne pas être calculées ... mais simplement déduites ... gains significatifs en temps calcul et en précision des résultats.
Evaluation polynômiale
Evaluer un polynôme de degré n peut se faire en O(n) opérations.Cela a aussi pour effet de limiter les erreurs d'arrondi en cumulant les contributions de chaque terme du polynôme. Pour le faire on écrit le polynôme suivant sous la forme suivante :
De manière générale, l'évaluation du polynôme se fait à l'aide de la boucle suivante :
INTEGER :: n, i
REAL(rp), DIMENSION(0:n) :: P
REAL(rp) :: s (valeur du polynôme en x)
s = P(n)
DO i = n-1, 0, -1
s = s * x + P(i)
END DO
Arithmétique entière
Comme les processeurs RISC sont surtout orientés pour le calcul flottant, certaines opérations effectuées sur les nombres entiers ne sont pas optimisées dans les compilateurs : la division est très chère à effectuer.Dans le cas de divisions par puissance de 2, on a un moyen efficace de les effectuer : la fonction ishft.
| Au lieu d'écrire ... | on peut mettre |
| k = i*(i-1)/2 | k = ISHFT( i*(i-1), -1 ) |
| l = j*m/8 | l = ISHFT( j*m, -3 ) |
Remarque :
L'expression i*(i-1)/2 est dangereuse car elle ne donne pas forcément le même résultat selon l'ordre dans lequel les opérations sont effectuées : on sait que i ou i-1 est pair donc le produit est pair ; la division entière par 2 se passe bien.
Dans le cas où i est pair, alors (i-1) est impair et si la division a lieu avant la multiplication, i.e. le calcul de (i-1)/2 avant celui de i*(i-1), le résultat est faux. Cela peut survenir lorsque le niveau d'optimisation est assez haut. Il ne faut pas alors se limiter à un niveau inférieur d'optimisation mais plutôt éliminer ces instructions potentiellement dangereuses.
Le premier argument de la fonction ishft est le nombre à diviser, le
second est l'exposant de la puissance de 2.
On peut noter que cela fonctionne aussi très bien pour la multiplication d'un entier par une puissance de 2 :
Calcul des constantes flottantes
Le calcul des constantes scalaires peut entrainer des pertes de précision significatives. Ainsi il est dommage de travailler avec des variables en double précision si les constantes scalaires sont évaluées en simple précision. On perd en précision et en temps calcul puisque les opérations en double précision coûtent au moins aussi cher que celles en simple précision.
Exemples :
INTEGER, PARAMETER :: rp = kind(1.0D0)
REAL(rp) :: PI_S, a_s, b_s
REAL(rp) :: PI_D, a_d, b_d
PI_S = ACOS( -1.0 ) <-- PI_S = 3.14159274101257324
PI_D = ACOS( -1.0_rp ) <-- PI_D = 3.14159265358979312
a_s = 0.5 * SQRT( 2.0 ) <-- a_s = 0.707106769084930420
a_d = 0.5 * SQRT( 2.0_rp ) <-- a_d = 0.707106781186547573
b_s = 2.0 / 3.0 <-- b_s = 0.666666686534881592
b_d = 2.0_rp / 3.0_rp <-- b_d = 0.666666666666666630
Cela se complique lorsque certaines options de compilation donnent les mêmes résultats dans les deux cas alors que d'autres distinguent les calculs en simple précision des calculs en double. En effet, ce genre de constantes scalaires est évalué lors de la compilation du programme.
Le tableau ci-dessous reprend quelques exemples :
| Séquence de compilation | Résultats identiques ? |
|---|---|
| xlf90_r -qfree=f90 test.f | VRAI |
| xlf90_r -qfree=f90 -O2 test.f | VRAI |
| xlf90_r -qfree=f90 -O3 test.f | VRAI |
| xlf90_r -qfree=f90 -qhot -O3 test.f | FAUX uniquement pour 1/2 * SQRT(1/2) |
| xlf90_r -qfree=f90 -qhot -O3 -qstrict test.f | FAUX uniquement pour 1/2 * SQRT(1/2) |
| xlf90_r -qfree=f90 -qarch=pwr4 test.f | FAUX |
| xlf90_r -qfree=f90 -qarch=pwr4 -qautodbl=dbl4 test.f | VRAI |
Conclusion :
toute seule, sans autre optimisation ou option, l'option -qarch arrive à modifier les résultats !
L'option -qhot ne calcule pas les constantes scalaires de la même manière lorsque la fonction SQRT intervient.
L'option -qautodbl=dbl4 en convertissant toutes les constantes en double précision avant évaluation corrige le problème.
De manière générale, il est conseillé de s'assurer de la précision des calculs dès l'écriture du programme, en imposant explicitement la précision souhaitée, comme dans l'exemple Fortran90 ci-dessus.
Tests arithmétiques flottants selon la précision machine
Bon : IF ( ABS(reel1 - reel2) .LT. machep) THEN ...
Mauvais : IF ( reel1 .EQ. reel2 ) THEN ...
avec machep une valeur de référence pour désigner une quantité négligeable :
machep = 1.0_rp
100 IF (machep + 1.0_rp .GT. 1.0_rp ) THEN
machep = 0.5_rp * machep
GO TO 100
END IF
machep = 10.0_rp * machep
machep est proche de 1E-6 en simple précision et 1E-14 en double précision.
On peut aussi utiliser les fonctions Fortran90 comme EPSILON.
Sauvegarde et reprise
La sauvegarde consiste à générer des fichiers de résultats pour une étude postérieure à l'exécution. On choisira alors une sortie formatée avec des descripteurs fournissant une précision adéquate : 4 ou 5 décimales sont suffisantes pour une représentation graphique classique.Une reprise consiste à générer des fichiers non formatés contenant toutes les informations (variables, paramètres, etc ..) nécessaires à la poursuite d'un calcul sur la même architecture comme s'il n'y avait pas eu d'interruption. Il faut donc faire en sorte que toute l'information présente en mémoire se retrouve dans le fichier sans perte de précision. Prévoir des reprises fréquentes est impératif que ce soit pour pallier à une interruption d'exécution (dépassement du temps cpu alloué dans la requête Loadleveler, par exemple) ou bien lorsque l'on souhaite poursuivre une simulation (convergence insuffisante, par exemple), sans pour autant devoir la relancer depuis l'instant initial.
D'autre part, pour les exploitations graphiques, il existe des formats spécifiques supportés par de nombreux logiciels :
- HDF est un format de données développé par le NCSA (National Center for Supercomputing Applications). La liste des logiciels supportant ce format comprend entre autres MATLAB, PV-Wave et Tecplot.
- NetCDF est un format de données développé par UCAR (University Corporation for Atmospheric Research). La liste des logiciels supportant ce format comprend entre autres AVS, MATLAB et PV-Wave.