IMPLANTATION DES TRANSFORMÉES DE FOURIER RAPIDES (FFT)

Table des matières

1- Objectif de la modification
2- Principe de la modification
3- La modification en pratique
     A- Petite introduction
     B- Modifications de sous-routines déjà existantes
     C- Ajout de nouvelles sous-routines
4- Les résultats de la modification
    A- Résultats des simulations
    B- Temps des simulations
    C- Temps relatif pris par la sous-routine MAC4251
5- Construction de AUX1
6- Calcul au second ordre

1- Objectif de la modification

Il s'agissait d'augmenter la rapidité du programme de simulation Monte Carlo pour être en mesure d'effectuer, dans un temps raisonnable, des simulations pour un réseau de taille 16*16.
 

2- Principe de la modification

Le programme de simulation Monte Carlo Quantique passe la majeure partie de son temps à multiplier des matrices.  La multiplication de matrices carrées N*N pleines (i.e. qui n'ont pas de section nulle a priori) est un ordre de N3 opérations, tandis que la multiplication d'une matrice carrée N*N par une matrice diagonale N est un ordre de N2 opérations.

L'idée de l'implantation des FFT était de tirer profit de ce fait.  En effet, les matrices d'énergie cinétique sont diagonales dans l'espace de Fourier.  Ainsi, en appliquant judicieusement des transformées de Fourier, on peut transformer ces matrices pleines en matrices diagonales.  Bien sûr, les transformées de Fourier elles-mêmes nécessitent un certain nombre d'opérations.  Normalement, la transformation de Fourier d'une matrice carrée N*N pleine est un ordre de N3 opérations: il n'y a donc à présent pas de réduction du nombre d'opérations.  C'est là que les transformées de Fourier rapides (FFT) interviennent.  En effet, ces dernières permettent d'effectuer la même transformation en un ordre de N2log2N opérations.  Ainsi, la multiplication des matrices, qui étaient un ordre de N3 opérations, devient un ordre de N2log2N +N2 opérations.  Le seul problème est que N doit nécessairement être une puissance entière de 2 pour que l'algorithme des FFT fonctionne (i.e. il faut que N=2z, où z=1,2,3, etc.).
 

3- La modification en pratique

    A- Petite introduction

Tout d'abord, il est à noter que pour un réseau de taille N*M, les matrices à multiplier sont de taille NM*NM.  Dans ce qui suit, nous utiliserons les conventions suivantes: T représente la transformée de Fourier directe et T-1 représente la transformée de Fourier inverse.  Nous travaillons ici dans le langage matriciel, de sorte que pour une matrice A quelconque, il n'est pas vrai en général que A*T=T*A.  Ainsi, l'orthogonalité de la transformation de Fourier s'écrit T-1 = Tt , où t fait référence à la transposition.  Évidemment, T* T-1 = I, la matrice identité.  C'est de ce fait que l'on tire profit dans la sous-routine MAC4251.
 

    B- Modifications de sous-routines déjà existantes

a) module HOPPING: des tableaux ont été ajoutés; F qui est la matrice H diagonalisée dans l'espace de Fourier, FM qui est la matrice HM diagonalisée dans l'espace de Fourier et TEMP, qui est un tableau complexe pour les étapes intermédiaires du calcul.

b) sous-routine INITH3: suite au calcul de H, on calcule F:  F=T*H* T-1 ; idem pour FM suite au calcul de HM.  Cela nécessite évidemment plus d'opérations, mais puisque INITH3 n'est appelée qu'une seule fois, le temps requis est négligeable.

c) sous-routine MAC4251: c'est dans cette sous-routine que l'on peut sauver beaucoup d'opérations.
==>Exemple:  il faut effectuer la multiplication exp(H)*G
-Auparavant, on effectuait directement la multiplication par l'appelle d'une sous-routine MAC
-Maintenant, on effectue  exp(H)*G=T-1 *T*exp(H)* T-1 *T*G ;
puisque T*exp(H)* T-1 =exp(F), cela se réduit à  T-1 *exp(F)*T*G, où exp(F) est diagonale.

En tenant compte de toutes les sous-routines nécessaires à ces opérations, les changements effectués dans MAC4251 remplacent un ordre de (NM)3 + (NM)2 opérations par un ordre de 10(NM)2 + (NM)2 log2N + (NM)2 log2M.  Pour N=M=16, cela représente une diminution du nombre d'opérations d'un facteur 14!  Cela suppose que les sous-routines de Numerical Recipes qui font les FFT nécessitent bel et bien le nombre d'opérations spécifié.


     C- Ajout de nouvelles sous-routines

a) FFTR: cette sous-routine applique la FFT (directe si ISIGN=1, inverse si ISIGN= -1) à gauche (L=0) ou à droite (L=2) sur le tableau réel DATR (NM,NM) passé en argument.

On commence par construire le tableau complexe TEMP à partir du tableau réel DATR:
TEMP(i,j)=cmplx (DAT(i,j), 0.d0) pour tous les i,j = 1,...,NM.

Ensuite, la sous-routine appelle DROITE (ISIGN) ou GAUCHE (ISIGN).

b) FFTC: cette sous-routine applique la FFT (directe si ISIGN=1, inverse si ISIGN= -1) à gauche (L=0) ou à droite (L=2) sur le tableau complexe TEMP.

La sous-routine appelle DROITE (ISIGN) ou GAUCHE (ISIGN).

c) MULD_TEMP: cette sous-routine multiplie TEMP par la matrice diagonale F (D=1) ou par la matrice diagonale FM (D= -1) par la gauche (L=0) ou par la droite (L=2).

d) DROITE: cette sous-routine applique la FFT (directe si ISIGN=1, inverse si ISIGN= -1) à droite.

On commence par construire le tableau AUX1(N,M) à partir de la i-ème rangée de TEMP.  Voir plus bas l'explication portant sur la construction de AUX1.

Ensuite, la sous-routine appelle FOURROW (AUX1,ISIGN) et FOURCOL (AUX1,ISIGN) qui effectuent respectivement la FFT selon y et selon x.

Puis, on effectue la reconstruction de la i-ème rangée de TEMP à partir de AUX1 (qui est maintenant transformée).

On recommence le même processus pour toutes les rangées de TEMP.

Finalement, si ISIGN= -1, on divise le tableau TEMP par la constante de normalisation NM apparaissant dans la définition de la transformée de Fourier discrète inverse car les sous-routines FOURROW et FOURCOL ne le font pas automatiquement.

e) GAUCHE: cette sous-routine applique la FFT (directe si ISIGN=1, inverse si ISIGN= -1) à gauche.

On commence par construire le tableau AUX1(N,M) à partir de la i-ème colonne de TEMP.  Voir plus bas l'explication portant sur la construction de AUX1.

Ensuite, la sous-routine appelle FOURROW (AUX1,ISIGN) et FOURCOL (AUX1,ISIGN) qui effectuent respectivement la FFT selon y et selon x.

Puis, on effectue la reconstruction de la i-ème colonne de TEMP à partir de AUX1 (qui est maintenant transformée).

On recommence le même processus pour toutes les colonnes de TEMP.

Finalement, si ISIGN= -1, on divise le tableau TEMP par la constante de normalisation NM apparaissant dans la définition de la transformée de Fourier discrète inverse car les sous-routines FOURROW et FOURCOL ne le font pas automatiquement.

f) FOURROW: cette sous-routine provient de Numerical Recipes in Fortran 90; elle remplace chaque rangée (les éléments du tableau à deux dimensions ayant la même première dimension) de la matrice complexe DAT par sa transformée de Fourier (directe si ISIGN=1, inverse si ISIGN= -1).  Il est à noter que tous les éléments de la matrice (tableau) seront NM fois trop grands si c'est la FFT inverse qui a été effectuée.

g) FOURCOL: cette sous-routine provient de Numerical Recipes in Fortran 90; elle remplace chaque colonne (les éléments du tableau à deux dimensions ayant la même deuxième dimension) de la matrice complexe DAT par sa transformée de Fourier (directe si ISIGN=1, inverse si ISIGN= -1).  Il est à noter que tous les éléments de la matrice (tableau) seront NM fois trop grands si c'est la FFT inverse qui a été effectuée.


4- Les résultats de la modification

Il est possible de voir les résultats de certaines simulations faites pour comparer les versions avec et sans FFT.  Il y a trois remarques importantes à faire, qui portent sur les résultats des simulations, le temps des simulations et le temps relatif pris par la sous-routine MAC4251.

    A- Résultats des simulations

En principe, pour un même germe (i.e. pour une même configuration initiale aléatoire du champ Hubbard-Stratonovich), les simulations selon les deux méthodes devraient donner exactement les mêmes résultats.  Bien sûr, on s'attend à une légère différence provenant des imprécisions numériques.  Cependant, on observe dans plusieurs cas que les résultats sont beaucoup plus éloignés que ne le laisseraient prévoir de simples instabilités numériques.  En fait, lorsque le nombre de mesures est élevé, ce phénomène semble se produire de façon presque automatique.

L'explication de cet écart est la suivante.  Lors de l'exécution du programme, les tests pour déterminer si un spin-flip sur un site sera effectué ou non sont fait un très grand nombre de fois (cliquez ici pour savoir ce qu'est un spin-flip).  En fait, ce test est effectué (mes+warm)*tran*sites fois, où mes est le nombre de mesures, warm le nombre de réchauffements, tran le nombre de tranches de temps imaginaire et sites le nombre de sites.

Pour savoir si le spin-flip est accepté, on compare une probabilité calculée à partir de la fonction de Green à un nombre aléatoire.  Les probabilités calculées selon les deux méthodes devraient en principe être identiques, mais il y a de légère différences numériques.  Il peut donc arriver que le nombre aléatoire ALEA se situe entre la probabilité P(a) calculée avec une méthode et la probabilité P(b) calculée avec l'autre méthode.  Si c'est le cas, le changement sera accepté dans un cas et refusé dans l'autre.  À partir de ce moment, les configuration du champ Hubbard-Stratonovich ne seront plus les mêmes dans les deux cas et il sera alors normal d'obtenir des résultats qui diffèrent.  Étant donnée la très légère différence initiale entre P(a) et P(b), il peut sembler très improbable que ALEA se situe entre les deux.  Cependant, il faut rappeler que le test de spin-flip est effectué un très grand nombre de fois.

Des tests ont été fait pour vérifier que cette explication est en effet la bonne.  On comprend ainsi pourquoi la divergence entre les deux méthodes semble presque certaine lorsqu'un grand nombre de mesures sont effectuées.
 

    B- Temps des simulations

Pour ce qui est des 8*8, les tests que nous avons effectué montre que la méthode avec les FFT est légèrement plus lente que la méthode sans FFT.  Cela est anormal, car dans ce cas, la sous-routine MAC4251 devraient demander 4 fois moins d'opérations avec les FFT que sans.  Cependant, il faut dire que le nombre d'opérations n'est pas uniquement ce qui détermine la vitesse d'un programme.  Il y a des effets d'accès à la mémoire, de compilateur, de rapidité d'exécution de certaines opérations comparativement à d'autres qui font qu'il est à toute fin pratique impossible de prévoir d'avance les effets d'une telle modification sur le temps de calcul.  En faisant des simulations sur d'autre noeuds, nous avons obtenu des résultats où cette fois la méthode avec FFT était notablement plus rapide que l'autre, même pour un 8*8.

Malheureusement, pour les 16*16, nous ne pouvions utiliser ces autres noeuds plus avantageux pour les FFT car il n'y avait qu'un seul noeud (elix01) qui comportait suffisamment de mémoire pour pouvoir effectuer ces simulations.  Pour cette taille, la méthode avec les FFT était en moyenne 2.2 fois plus rapide.  Cela n'est pas suffisant pour pouvoir envisager de faire d'importantes simulations.

Pour pouvoir augmenter la rapidité de la méthode avec les FFT, il restait donc deux pistes de travail.  Premièrement, trouver le noeud sur lequel les FFT sont les plus rapides.  Deuxièmement, trouver la sous-routine la plus rapide possible pour effectuer les FFT.

Cependant, suite à la réalisation d'un profil d'une simulation pour un réseau 16*16, nous avons réalisé que le temps relatif pris par la sous-routine MAC4251 ne nous permettrait pas d'améliorer suffisamment la vitesse des simulations pour faire des réseaux 16*16.
 

    C- Temps relatif pris par la sous-routine MAC4251

Comme le temps pris par la sous-routine INITH3 est négligeable et que les autres modifications de sous-routine ne concernent que MAC4251, c'est dans l'exécution de cette dernière que l'on peut sauver du temps.  On suppose (avec beaucoup de bon sens) que les sous-routines qui n'ont pas été modifiées prendront à peu près le même temps (il peut y avoir des effets d'accès à la mémoire étant donné la déclaration de nouveaux tableaux dans le module HOPPING, mais ces effets ne devraient pas être très importants).  Si par exemple MAC4251 représente 75% du temps total de calcul, il sera impossible en théorie d'améliorer la rapidité du programme par un facteur supérieur à 4.

Nous avons donc réalisé des profils pour des réseaux 4*4, 8*8 et 16*16 pour voir l'évolution du temps relatif pris par la sous-routine MAC4251.  Pour un 4*4, MAC4251 représentait environ 40% du temps de calcul total, tandis que pour un 8*8, c'était environ 50% du temps.  Nous espérions donc que pour un 16*16, le temps relatif soit encore plus important.  Malheureusement, il est retombé à 41.2% !  On ne peut donc espérer augmenter la vitesse par un facteur supérieur à 1.7!

En fait, comme mentionné dans le temps des simulations, le programme avec FFT était en moyenne 2.2 fois plus rapide.  Cependant, le profil nous a permis de réaliser qu'il serait vain d'espérer pouvoir accélérer autant les programmes que nous l'aurions voulu.  Cela tient également au fait que l'améliration principale, qui aurait porté sur les sous-routines qui effectuent les FFT, aurait été négligeable.  En effet, le profil nous informe que les sous-routines FOURROW et FOURCOL ne représentent au total que 8.5% du temps de calcul.  Même avec des sous-routines infiniment plus rapides, le changement serait minime.

Il peut sembler curieux que le programme avec FFT était en moyenne 2.2 fois plus rapide alors qu'on ne peut théoriquement pas espérer en augmenter la vitesse par un facteur supérieur à 1.7.  Cette anomalie provient du fait que le profil est nécessairement généré sans option d'optimisation, alors que les test ont été effectués avec des programmes optimisés.  Un autre test effectué à partir de programmes non optimisés donne plutôt une accélération d'un facteur 1.7, exactement la valeur maximale.  Bien sûr, cette valeur est une limite que nous ne devrions pas atteindre, mais il faut se rappeler trois choses.  Premièrement, ce facteur provient d'une estimation du nombre d'opérations et non pas d'un calcul exact de ce nombre.  Deuxièmement, deux opérations différentes ne prennent pas nécessairement le même temps.  Autrement dit, ce n'est pas parce qu'il y a deux fois plus d'opérations qu'automatiquement le temps est exactement le double.  Troisièmement, il y a toujours des effets internes (utilisation de la "cache", etc.) dont on peut difficilement tenir compte pour faire une comparaison précise des temps de calcul.
 

5- Construction de AUX1

Si on multiplie une matrice A par la droite avec une matrice B (i.e. A*B), ce seront les rangées de A et les colonnes de B qui interviendront.  Voilà pourquoi dans la sous-routine DROITE, on construit AUX1 à partir des rangées de TEMP (et dans la sous-routine GAUCHE, à partir des colonnes de TEMP).

Dans une rangée (colonne) de TEMP, l'information relative aux NM sites du réseau est contenue.  Ce sont les coordonnées selon x qui varient en premier.
Ainsi, les N premiers éléments de la rangée (colonne) correspondent aux sites (x,y) suivants:
(1,1); (2,1); (3,1); ...; (N,1)
Les N éléments suivants (i.e. les éléments N+1 à 2N) correspondent aux sites (x,y) suivants:
(1,2); (2,2); (3,2); ...; (N,2)
Et ainsi de suite jusqu'aux N derniers éléments correspondant aux sites (x,y) suivants:
(1,M); (2,M); (3,M); ...; (N,M)

Le tableau AUX1 (N,M) est formé comme suit:
-la première colonne contient les N premiers éléments de la i-ème rangée (colonne) de TEMP
-la deuxième colonne contient les N éléments suivants de la i-ème rangée (colonne) de TEMP
-et ainsi de suite jusqu'à la Mième colonne, qui contient les N derniers éléments de la i-ème rangée (colonne) de TEMP

Dans un même colonne de AUX1, les éléments correspondent donc tous au même site y.  De même, dans une même rangée de AUX1, les éléments correspondent tous au même site x.  C'est pourquoi on dit que FOURROW fait la transformée de Fourier selon y (qui varie dans une même rangée) et que FOURCOL fait la transformée de Fourier selon x (qui varie dans une même colonne).
 

6- Calcul au second ordre

Suite à l'automatisation du calcul au second ordre par Steve Allen, nous avons décidé d'implanter également les FFT pour ce dernier.  Les modifications supplémentaires du programme suite à cet ajout sont les suivantes:
 

a) module HOPPING: des tableaux ont été ajoutés; F2M qui est la matrice H2M diagonalisée dans l'espace de Fourier et F2P qui est la matrice H2P diagonalisée dans l'espace de Fourier.

b) sous-routine INITH3: si il y a lieu d'avoir calcul au second ordre (voir automatisation du calcul au second ordre), voici les ajouts.  Suite au calcul de H2P, on calcule F2P:
F2P=T*H2P* T-1 ; idem pour F2M suite au calcul de H2M.

c) sous-routine MULD_TEMP: puisque l'on peut maintenant multiplier TEMP par F2M ou par F2P, il existe deux options de plus à MULD_TEMP.  Si D= -2, on multipliera par F2M; si D=2, on multipliera par F2P.

d) sous-routine MACGD: c'est dans cette sous-routine que l'on peut sauver beaucoup d'opérations.  Le principe est exactement le même que pour MAC4251, sauf qu'on utilise les matrices de deuxième ordre (F2M et F2P) au lieu des matrices de premier ordre (F et FM).