<-Précédent Retour à l'accueil Contact : etienne"point"sauvage"at"gmail.com
  1. De quoi que c'est-y qu'on cause ?
  2. De quoi avons-nous besoin ?
  3. Est-ce qu'à part frimer devant la machine à café, ça a une utilité ?
  4. Mais pourquoi, comment ?
    1. De la représentation des matrices et de l'arithmétique des pointeurs.
    2. De l'instrument transpositeur.
  5. Des annoncées fonctions
  6. Ta. Dam.

    Assembleur : calculs en virgule flottante

    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é.

  1. De quoi que c'est-y qu'on cause ?
  2. 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.

  3. De quoi avons-nous besoin ?
  4. 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.

  5. Est-ce qu'à part frimer devant la machine à café, ça a une utilité ?
  6. La multiplication matricielle est l'une des briques de base du traitement du signal, par exemple. Donc non, vraiment, ça ne sert à rien.

  7. Mais pourquoi, comment ?
    1. De la représentation des matrices et de l'arithmétique des pointeurs.
    2. 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.

    3. De l'instrument transpositeur.
    4. 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.

  8. Des annoncées fonctions.
  9. 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
    

  10. Ta. Dam.