| <-Précédent | Retour à l'accueil | Contact : etienne"at"les-sauvages.fr | Suivant-> |
Nous allons nous attaquer à l'une des fonctions phares du calcul matriciel, celle pour laquelle on a fait le moins de progrès depuis son invention : la multiplication matricielle. Je crois que c'est du niveau des classes propédeutiques, ce qui ce traduit depuis mon jargon suranné en maths sup' - maths spé.
Une matrice est un tableau de nombres qui ont un sens ensemble. Les vrais matheux définissent ça, je crois, comme une application affine, mais je leur laisse l'usufruit, la pleine propriété, la rente et la jouissance de la vraie définition. Moi, je sais qu'une matrice, c'est un tableau de nombres qui vont ensemble, comme les valeurs des pixels d'une image, l'historique des cours du marché A de la bourse, le tableau des notes d'une classe, voire l'évolution de mes mensurations en fonction du temps. Oui, la vieillesse, c'est la victoire de la gravité.
Bon, on a depuis quelques siècles appris à travailler avec ces matrices, et on en a défini les opérations usuelles. On en vient donc à la multiplication. La multiplication matricielle est assez spéciale : elle est assez simple dans son esprit, mais terriblement longue et ennuyeuse à calculer. Pour multiplier deux matrices carrées de taille n, on fait, en appliquant la définition de la multiplication, n x n x n multiplications, et quasiment autant de sommes. Pour n = 2, ça fait déjà 8 multiplications à se farcir. Pour n = 3, on passe directement à 27. Une imagette de 100 par 100, ça fait 1.000.000 (un million) de multiplications. Sérieusement, les gars : un million de multiplications pour une ch'tiote image de 100 pixels, vous vous imaginez les faire, à la main ?
Non, moi non plus. Le mieux qu'on ait fait niveau algorithmique, c'est enlever une multiplication sur 8 au prix d'une brassée d'additions en plus. Et l'usage de la multiplication matricielle est, comment dire... courant. Par exemple, une scène 3D utilise la multiplication matricielle quelques (de 1 à ...plus que ça) fois par image : les rotations d'object se font par multiplication matricielle, ainsi que les changements de point de vue. Autant dire que Counter Battle Royale : la liche 4.2 dark age golden edition, elle vous en sort 6000 par seconde.
Récapitulons : un éditeur de texte (j'utilise Geany), YASM, gcc, un système POSIX (j'utilise Debian), une machine 64 bits et BLAS.
BLAS se prononce "sous-programmes basiques d'algègre linéaire", Basic Linear Algebra Subprograms. Sous-programme veut dire fonction, basique signifie qu'on n'y trouvera pas la transformation submatricelle de Zanvisky à coefficients transcendants, et algèbre linéaire qu'on ne devrait pas trop trouver de sinus et exponentielles en échange de la possibilité de travailler avec des vecteurs et des matrices. Nous allons l'utiliser parce que la méthode naïve en C de calcul matriciel va, je vous vends la mèche, se faire éclater de chez éclater par le code assembleur. Pour redonner un peu de tenue aux langages de haut niveau, nous allons nous comparer à BLAS, référence en terme de rapidité de calcul. BLAS, il y en a plein de versions j'ai choisi CBLAS parce qu'on va l'appeler en C.
La multiplication matricielle est l'une des briques de base du traitement du signal, par exemple. Donc non, vraiment, ça ne sert à rien.
Vous savez pourquoi Trinity est plate comme une limande ? Parce que dans la Matrice, on n'est qu'en une seule dimension.
Je ne sais pas si vous avez remarqué, mais la mémoire est accessible via un simple index, genre mov rax, [rdx]. C'est un monde monodimensionnel. Nos matrices sont bidimensionnelles, une ordonnée et une abscisse. On doit en conséquence "linéariser" la matrice. On peut les représenter de plusieurs façons, mais la plus courante est celle-ci :
La matrice A(n, m) est représentée par un tableau de nombres de dimension n x m. L'élément A(i, j) est stocké à l'index (i * m + j) du tableau. Une autre convention possible est de stocker l'élément à l'index (i + j * n).
On aura toujours une opération à faire sur les index pour accéder à un élément arbitraire d'une matrice : l'index doit en plus être multiplié par la taille de l'élément. On parle de ceci comme étant l'arithmétique des pointeurs.
Outre les techniques dont on commence à avoir l'habitude, l'astuce va consister à tripoter l'une des matrices avant de faire les calculs.
L'algorithme naïf (celui qu'on va coder) indique que le résultat à la ligne i et la colonne j est le produit scalaire du vecteur ligne i du multiplicateur et du vecteur colonne j du multiplicande. L'important ici est "vecteur ligne" puis "vecteur colonne". On a deux vecteurs à utiliser, l'un dans un sens et l'autre dans l'autre sens. On ne pourra remplir un registre vectorisé que valeur par valeur, parce que pour l'un des vecteurs, les valeurs ne sont pas contigues en mémoire. On parcourt A et B selon j : on a A [i * m + j], valeurs contigues, et B [j * n + k], valeurs non contigues. Et accéder à des valeurs non contigues en mémoire, c'est long.
Pour ce faire, nous allons transposer l'une des matrices. Transposer, ça veut dire faire l'opération qui transformera les lignes en colonnes et les colonnes en lignes. La transposée de la matrice B se note B'. On a B'[i, j] = B[j, i]
Voici la fonction
transpose: ;On va dire : RDI = a', RSI = a, RDX = nombre de lignes de a, RCX = nombre de colonnes de a mov r9, rcx mov r10, rdi mov r8, rdx shl rdx, 3 .boucle: lodsq mov [rdi], rax add rdi, rdx loop .boucle mov rcx, r9 add r10, 8 mov rdi, r10 dec r8 jnz .boucle ret
Donc dans RDI, on a l'adresse de la matrice transposée, sous prétexte que le D de RDI est là pour Destination. Cette fonction ne fait pas l'allocation mémoire, donc il faut faire attention à la validité de son pointeur. Dans RSI, avec un S comme Source, on a la matrice à transposer. Dans RDX, le nombre de lignes et dans RCX le nombre de colonnes.
Donc la première chose que l'on fait, c'est de sauvegarder ces valeurs dans des registres qui ne craignent pas. C'est inutile pour RSI cependant : on ne va parcourir l'entrée qu'une fois. C'est une fonction en O(n).
Ensuite on multiplie RDX par 8 pour avoir le nombre d'octets dans une colonne. Oui, une colonne a bien nombre de lignes éléments.
On boucle sur une ligne : on lit un double dans RAX à partir de RSI et on le stocke dans RDI. Ensuite, on fait pointer RDI sur la ligne suivante de la transposée.
On recharge le compteur de boucle, on remet RDI sur la ligne suivante, et on reboucle.
J'ai dit que la multiplication matricielle était constituée de produits scalaires. Le produit scalaire de deux vecteurs A et B est la somme des multiplications A[i] * B[i]. Zou :
dotProduct: ;Ordre des paramètres : RDI, RSI, RDX, RCX, R8 et R9 ;On va dire : RDI = c, RSI = a, RDX = b, rcx = n xorpd xmm1, xmm1 xor rax, rax test rcx, rcx jz .retour shr rcx, 1 jnc .somme jz .impair inc rax test rdx, 0xF jnz .bNonAligne test rsi, 0xF jnz .aNonAligne .somme: movapd xmm0, [rsi] mulpd xmm0, [rdx] addpd xmm1, xmm0 add rsi, 16 ;On passe à la paire de double suivante add rdx, 16 loop .somme ;Et on boucle sur RCX .general: test rax, 1 jz .enregistrement .impair: movsd xmm0, [rsi] mulsd xmm0, [rdx] addsd xmm1, xmm0 add rsi, 8 ;On passe au double suivant add rdx, 8 .enregistrement: movhlps xmm0, xmm1 addsd xmm1, xmm0 movsd [rdi], xmm1;On stocke le résultat dans la matrice : on vient de calculer un terme. .retour: ret .bNonAligne: test rsi, 0xF jnz .aBNonAlignes .sommeBNonAligne: movupd xmm0, [rdx] mulpd xmm0, [rsi] addpd xmm1, xmm0 add rsi, 16 ;On passe à la paire de double suivante add rdx, 16 loop .sommeBNonAligne ;Et on boucle sur RCX jmp .general .aBNonAlignes: .sommeNonAligne: movupd xmm0, [rsi] movupd xmm2, [rdx] mulpd xmm0, xmm2 addpd xmm1, xmm0 add rsi, 16 ;On passe à la paire de double suivante add rdx, 16 loop .sommeNonAligne ;Et on boucle sur RCX jmp .general .aNonAligne: .sommeANonAligne: movupd xmm0, [rsi] mulpd xmm0, [rdx] addpd xmm1, xmm0 add rsi, 16 ;On passe à la paire de double suivante add rdx, 16 loop .sommeANonAligne ;Et on boucle sur RCX jmp .general
AAAAAAAAAAAAAAhhhhhhhh !! Tout ça pour un bête produit scalaire ? Arf, oui. Non mais sérieux, t'as mangé un clown ? A quoi ça sert tout ça ?
Je conçois votre étonnement, ami lecteur. Vous souvenez-vous de l'alignement mémoire ? Il se trouve que nous ne pouvons plus faire l'économie, dans cette fonction, de la détection et du traitement des cas où les vecteurs ne sont pas alignés. Afin d'utiliser les instructions les plus efficaces, nous avons quatre cas :
On teste l'alignement ainsi : si l'un des 4 derniers bits de RSI (resp. RDX) est à 1, alors A (resp. B) n'est pas aligné. La subtilité se base sur le fait qu'on peut mettre nue adresse mémoire alignée en paramètre de mulpd.
Du coup, il ne reste plus qu'à mettre tout ça ensemble : on transpose, on boucle pour les produits scalaires :
crossProduct: ;Ordre des paramètres : RDI, RSI, RDX, RCX, R8 et R9 ;On va dire : RDI = c, RSI = a, RDX = b, rcx = n, r8 = M, R9 = P push rbx push rbp push r12 push r13 push r14 push r15 mov rbx, rdi ;sauvegarde de c mov rbp, rsi ;sauvegarde de a mov r12, rdx ;sauvegarde de b mov r13, rcx ;sauvegarde de n mov r14, r8 ;sauvegarde de m mov r15, r9 ;sauvegarde de p ;Appel de malloc pour avoir une matrice transposée mov rax, rcx shl rax, 3 mul r9 mov rdi, rax call malloc ;Calcul de la transposition mov rdi, rax mov rsi, r12 ;sauvegarde de RDX : b mov r11, rdi mov rdx, r13 mov rcx, r15 call transpose mov rdi, rbx mov r8, r14 ;sauvegarde de R8 : m mov r9, r15 ;sauvegarde de R9 : p shl r9, 3 .boucle: mov rdx, r11 ;sauvegarde de RDX : b mov r14, r15 ;sauvegarde de R9 : p .boucleLignes: mov rsi, rbp ; restauration de a mov rcx, r13 call dotProduct add rdi, 8 dec r14 jnz .boucleLignes mov rdx, r13 shl rdx, 3 add rbp, rdx dec r8 jnz .boucle .retour: mov rdi, r11 call free pop r15 pop r14 pop r13 pop r12 pop rbp pop rbx ret