Vectorialisation sous Matlab

Sommaire

I. Le logiciel MatLab

MatLab, contraction de Matrix et de Laboratory, est un logiciel optimisé pour la calcul matriciel. Cela a pour conséquence qu'un programme, pouvant être fait soit à l'aide d'une ou plusieurs boucles, soit à l'aide d'un calcul matriciel, s'exécutera généralement plus vite en utilisant un calcul matriciel puisque, dans ce cas, MatLab accède aux données dans un ordre adapté à leur stockage en mémoire.
La gestion des boucles sous MatLab s'est améliorée avec les dernières versions ; cependant, il est souvent plus avantageux de vectorialiser ses calculs de manière à obtenir dans un code une majorité de calculs matriciels à la place des boucles.

II. Configuration matérielle pour les tests

Processeur: 1,6 GHz
RAM: 2 Go
Swap: 4 Go

III. Comment vectorialiser un calcul

Tout d'abord, la plupart des fonctions MatLab peuvent agir sur des vecteurs et non seulement sur des simples nombres ; par exemple, pour calculer les valeurs du logarithme népérien pour les entiers de 1 à 10^6, stockées dans un vecteur ligne v, il suffit de faire:

v=log((1:1e6));

Vous pouvez alors vérifier, à l'aide des quelques lignes suivantes, que cette opération est plus rapide qu'une boucle.

tic;
v=zeros(1,1e6);
for p=1:1e6
    v(p)=log(p);
end
toc;

Le temps d'exécution pour cette boucle est de 3,7s, alors que le temps d'exécution du même calcul écrit de cette manière

tic;v=log((1:1e6));toc;

est de 0,3s.

Soit un facteur 10 entre les deux temps d'exécution.

Ensuite, beaucoup de fonctions MatLab permettent de s'affranchir de calculs à l'aide de boucles. En voici une liste non exhaustive:

• Les opérateurs arithmétiques *, /, ^ peuvent agir sur des tableaux de même taille en effectuant l'opération terme à terme ; il suffit de rajouter un point devant l'opérateur pour préciser qu'on veut réaliser cette opération, c'est-à-dire .*, ./, .^. Par exemple, (1:10).^2 donnera les 10 premiers carrés parfaits.
• repmat - très simple d'utilisation, tapez help repmat ou doc repmat pour avoir une aide quant à cette fonction - permet de répéter un tableau (ce tableau peut être un vecteur ou une matrice) pour en construire un plus grand.
• reshape permet de redimensionner un tableau ; l'utilisation de cette fonction implique de connaître quelque peu comment MatLab "replace" les éléments lors d'un redimensionnement. De nouveau, l'aide MatLab donnée par help reshape ou doc reshape vous sera utile, vous pouvez consulter l'astuce sur les manipulations élémentaires des tableaux pour obtenir des explications différemment formulées.
• Beaucoup d'autres !

IV. Expérience sur un exemple concret

Dans le but d'illustrer l'avantage de la vectorialisation, travaillons sur un exemple concret. Savoir ce que fait le programme que nous allons utiliser n'est pas une nécessité absolue, il a été choisi de manière à montrer que vectorialiser un calcul n'est pas un artifice destiné à améliorer la vitesse d'exécution de simples "cas d'école" mais peut servir dans des cas réels de programmation. Il est donc possible de passer directement aux sous-sections IV.2 et IV.3 sans que la compréhension de l'astuce en soit altérée.

IV.1. Exemple concret choisi

Le programme utilisé pour les tests calcule le champ magnétique généré dans l'espace par une boucle de courant.

Le champ magnétique généré par une telle boucle est donné par la loi de Biot et Savart:

IV.2. Programme non vectorialisé

function B = compute_B_loop(R,I,x,n)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calcule le champ magnétique généré par une boucle de courant
%% R est le rayon de la boucle
%% I est l'intensité du courant
%% x est la matrice des points où on veut calculer le champ B
%% La taille de x est  3 x nombre_de_points
%%     ( x1 ... ...)
%% x = ( y1 ... ...)
%%     ( z1 ... ...)
%% n est le nombre de points de discrétisation pour la boucle de courant 
%% B est une matrice 3 x nombre_de_points
%%     ( Bx1 ... ...)
%% B = ( By1 ... ...)
%%     ( Bz1 ... ...)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu0 = 4*pi*10^(-7);

nx=size(x,2);

B=zeros(3,nx);

for p=1:nx
    for q=0:n-1
        normxmy3=((x(1,p)-R*cos(2*pi/n*q))^2+(x(2,p)-R*sin(2*pi/n*q))^2+x(3,p)^2)^(3/2);
        B(1,p)=B(1,p)+cos(2*pi/n*q)*x(3,p)/normxmy3;
        B(2,p)=B(2,p)+sin(2*pi/n*q)*x(3,p)/normxmy3;
        B(3,p)=B(3,p)...
            +(-sin(2*pi/n*q)*(x(2,p)-R*sin(2*pi/n*q))-...
            cos(2*pi/n*q)*(x(1,p)-R*cos(2*pi/n*q)))/normxmy3;
    end
end


B = mu0*R*I/(2*n)*B;
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

IV.3. Programme vectorialisé

function B = compute_B(R,I,x,n)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calcule le champ magnétique généré par une boucle de courant
%% R est le rayon de la boucle
%% I est l'intensité du courant
%% x est la matrice des points où on veut calculer le champ B
%% La taille de x est  3 x nombre_de_points
%%     ( x1 ... ...)
%% x = ( y1 ... ...)
%%     ( z1 ... ...)
%% n est le nombre de points de discrétisation pour la boucle de courant 
%% B est une matrice 3 x nombre_de_points
%%     ( Bx1 ... ...)
%% B = ( By1 ... ...)
%%     ( Bz1 ... ...)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu0 = 4*pi*10^(-7);

nx=size(x,2);

%% x est la matrice des points de l'espace, on la répète
x=repmat(reshape(x.',1,nx,3),[n 1 1]);

%% y est la matrice des points de la boucle
y=R*[cos(2*pi/n*(0:n-1));sin(2*pi/n*(0:n-1));zeros(1,n)];

%% On répète la matrice y
y=repmat(reshape(y.',n,1,3),[1 nx 1]);

%% On peut alors calculer matriciellement tous les vecteurs SM (i.e. yx)
xmy=x-y; clear x y;

%% On calcule tous les vecteurs dl de la formule (on multipliera par R à la fin)
dy=[-sin(2*pi/n*(0:n-1));cos(2*pi/n*(0:n-1));zeros(1,n)];

%% On répète cette matrice
dy=repmat(reshape(dy.',n,1,3),[1 nx 1]);

%% On peut alors calculer matriciellement tous les vecteurs dl^SM de la formule (i.e. dy^yx)
cdyxmy=cross(dy,xmy,3); clear dy;

%% On calcule au passage toutes les normes de SM (yx)
norm_xmy_3=repmat(dot(xmy,xmy,3).^(3/2),[1 1 3]); clear xmy;

%% Il ne reste plus qu'à assembler B
B=cdyxmy./norm_xmy_3; clear cdyxmy norm_xmy_3;
B=sum(B);
B = mu0*R*I/(2*n)*B;
B = squeeze(B);
B = B.';
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

IV.4. Temps d'exécution pour des calculs importants

Exécutons maintenant ces deux programmes pour les points de [-5,5]x[-5,5]x[-5,5] sur une grille régulière 50x50x50 avec 50 points de discrétisation sur la boucle de courant.

L'exécution du programme non vectorialisé

[X Y Z]=meshgrid(linspace(-5,5,50));
x=[reshape(X,1,50^3);reshape(Y,1,50^3);reshape(Z,1,50^3)];
clear X Y Z;
tic;B_loop=compute_B_loop(1,1,x,50);toc;
clear B_loop;

nous donne un temps d'exécution de 20,8s.

L'exécution du programme vectorialisé

[X Y Z]=meshgrid(linspace(-5,5,50));
x=[reshape(X,1,50^3);reshape(Y,1,50^3);reshape(Z,1,50^3)];
clear X Y Z;
tic;B=compute_B(1,1,x,50);toc;
clear B;

nous donne un temps d'exécution de 14,4s.

Le programme non vectorialisé prend donc près de 1,5 fois plus de temps que le calcul vectorialisé.

Imaginez donc le résultat sur un programme qui tourne en quelques heures. C'est toujours plus plaisant d'avoir les résultats après "seulement" 3h que d'attendre 4h30. D'autant plus si les résultats sont faux et qu'il faut relancer le code après correction... ;-)

V. Inconvénients de la vectorialisation

• Comme nous pouvons le voir avec le programme compute_B, même commenté, un calcul vectorialisé est difficilement lisible pour une personne qui souhaiterait l'utiliser, le modifier, le compléter ou encore le traduire dans un autre langage (C, Fortran, etc...) pour l'inclure dans un code plus gros.
• Le deuxième inconvénient est bien évidemment l'espace mémoire. En effet, un programme vectorialisé a parfois besoin de stocker des matrices plus grosses que celles nécessaires pour stocker le résultat final, c'est en particulier vrai pour notre exemple. Tentez par exemple de lancer le calcul suivant (on prend maintenant une grille régulière 100x100x100).

[X Y Z]=meshgrid(linspace(-5,5,100));
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];
clear X Y Z;
tic;B=compute_B(1,1,x,50);toc;
clear B;

La place mémoire est insuffisante pour l'exécution du programme.

Alors que le programme non vectorialisé

[X Y Z]=meshgrid(linspace(-5,5,100));
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];
clear X Y Z;
tic;B_loop=compute_B_loop(1,1,x,50);toc;
clear B_loop;

s'exécute sans problème.

VI. Solutions aux précédents problèmes

• Pour le problème de lisibilité, il va de soi que commenter au maximum son code vectorialisé (même un code non vectorialisé bien sûr !) peut être une aide précieuse pour un éventuel futur utilisateur.
• Pallier au second problème est un peu plus compliqué. Il s'agit généralement de garder la vectorialisation pour faire passer les calculs "en blocs", mais de segmenter pour faire passer des blocs de taille adaptée pour ne pas saturer l'espace mémoire. Voici donc comment procéder pour notre exemple. Ici le plus gros nombre de points (celui qui posera donc problème) est celui qui correspond à la discrétisation du volume. Nous allons donc segmenter le bloc de points envoyé pour le calcul du champ magnétique.

function B = compute_B(R,I,x,n)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calcule le champ magnétique généré par une boucle de courant
%% R est le rayon de la boucle
%% I est l'intensité du courant
%% x est la matrice des points où on veut calculer le champ B
%% La taille de x est  3 x nombre_de_points
%%     ( x1 ... ...)
%% x = ( y1 ... ...)
%%     ( z1 ... ...)
%% n est le nombre de points de discrétisation pour la boucle de courant 
%% B est une matrice 3 x nombre_de_points
%%     ( Bx1 ... ...)
%% B = ( By1 ... ...)
%%     ( Bz1 ... ...)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu0 = 4*pi*10^(-7);

nx=size(x,2);
%% La taille l d'un bloc sera de 32768 points
l = 32768;
%% Les calculs suivants sont identiques au code précédent vectorialisé
%% sauf qu'on travaille par blocs
seg=(1:l:nx);
if seg(end)~=nx, seg(end+1)=nx; end
seg(end)=seg(end)+1;
nseg=size(seg,2)-1;
display([num2str(nseg), ' blocs']);
B = zeros(nx,3);
for pq=1:nseg
    display(['Bloc ',num2str(pq)]);
    npq=seg(pq+1)-seg(pq);
    xpq=repmat(reshape(x(:,seg(pq):seg(pq+1)-1).',1,npq,3),[n 1 1]);

    y=R*[cos(2*pi/n*(0:n-1));sin(2*pi/n*(0:n-1));zeros(1,n)];
    y=repmat(reshape(y.',n,1,3),[1 npq 1]);
    xpqmy=xpq-y; clear xpq y;
    dy=[-sin(2*pi/n*(0:n-1));cos(2*pi/n*(0:n-1));zeros(1,n)];
    dy=repmat(reshape(dy.',n,1,3),[1 npq 1]);
    cdyxpqmy=cross(dy,xpqmy,3); clear dy;
    norm_xpqmy_3=repmat(dot(xpqmy,xpqmy,3).^(3/2),[1 1 3]); clear xmy;

    B(seg(pq):seg(pq+1)-1,:)=squeeze(sum(cdyxpqmy./norm_xpqmy_3)); clear cdyxpqmy norm_xmy_3;
end
B = mu0*R*I/(2*n)*B;
B = B.';
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Lançons maintenant le calcul précédent qui ne s'exécutait pas faute d'insuffisance d'espace mémoire.

[X Y Z]=meshgrid(linspace(-5,5,100));
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];
clear X Y Z;
tic;B=compute_B(1,1,x,50);toc;
clear B;

Il lui faut environ 124s pour s'exécuter.


Le programme non vectorialisé

[X Y Z]=meshgrid(linspace(-5,5,100));
x=[reshape(X,1,100^3);reshape(Y,1,100^3);reshape(Z,1,100^3)];
clear X Y Z;
tic;B_loop=compute_B_loop(1,1,x,50);toc;
clear B_loop;

s'exécute quant à lui en 169s.

Ici encore, il faut au programme non vectorialisé environ 1,5 fois plus de temps qu'au programme "semi-vectorialisé" pour s'exécuter.