| <-Précédent | Retour à l'accueil | Contact : etienne"at"les-sauvages.fr |
Afin de pouvoir faire des calculs plus compliqués que les simples opérations arithmétiques sur des vecteurs, et surtout pour travailler avec les matrices, nous avons besoin d'une opération qui s'appelle le produit matriciel. Pour vous la faire courte, il s'agit de calculer l'ensemble des produits scalaires possibles entre les lignes d'une matrice A et les colonnes d'une matrice B.
Simple, me direz-vous. Reprenons le code du produit scalaire. Hélas il faudrait transposer la matrice B pour transformer ses colonnes en lignes, et ensuite on perdrait l'avantage de la contiguité de la mémoire.
Néanmoins vous avez tout à fait raison. En C++ des familles, l'opération s'écrit :
void crossProduct(int m, int n, int p, double* A, double* B, double* C)
{
for (auto i = 0; i < m; i++)
for (auto j = 0; j < p; j++)
{
C[i * p + j] = 0;
for (auto k = 0; k < n; k++)
{
C[i * p + j] += A[i * n + k] * B[k * p + j];
}
}
}
Hélas en assembleur, la procédure va être plus complexe.
Je rappelle qu'on travaille avec MASM de Visual Studio sous Windows, sur architecture 64 bits, avec un processeur supportant AVX2 et FMA. On se relaxe, on ouvre ses chakras.
Oui. On va s'ouvrir le champ des possibles en calcul matriciel, à nous les décompositions PLU, les vecteurs propres et autres compositions matricielles. L'éclate totale.
Regardez bien le code C montré plus haut. La fonction prend 6 arguments. Or, Microsoft a décidé que le nombre d'arguments qu'on peut passer à une fonction par les registres est de 4. Pas un de plus, moins si l'on veut. Or nous, nous en avons besoin de 6. Enfer et damnation, que va-t-il advenir de nos arguments supplémentaires ?
Hé bien, oyez bien braves gens, ils vont être passés sur la pile.
La pile ? Késaco ? La pile, mes bons amis, est un espace de stockage de toutes sortes de valeurs, orienté dans le sens des adresses décroissantes. À l'envers, quoi. Plus on met de choses sur la pile, et plus son adresse sera petite. Un registre général lui est consacré : RSP (SP comme Stack Pointer). Et jouer avec la pile, c'est un sport très couru des compilateurs et développeurs assembleurs.
Ceci étant dit, qu'advient-il de nos paramètres à nous ? Accrochez-vous c'est tordu : les 4 premiers sont passés par les registres RCX, RDX, R8 et R9, et les autres sont placés sur la pile. Mais attention, ça pourrait être trop simple. Avant de mettre les paramètres supplémentaires sur la pile, la fonction appelante doit réserver sur la pile un espace dans lequel on puisse stocker les 4 premiers paramètres en considérant qu'il s'agit de registres : on appelle ça le shadow store. On trouvera donc notre paramètre 5 à RSP + ((4 + 1) * 8) + 8 = 48 (on commence par stocker RBP), et le sixième à RSP + 56.
On se repose, pour les calculs, sur cette instruction : VFMADD231PD. Elle travaille notamment sur des registres YMM, et effectue la multiplication de ses opérandes 2 et 3, qu'elle ajoute et stocke dans le premier opérande.
PUBLIC crossProduct_asm
crossProduct_asm PROC
push rbp
mov rbp, rsp
push r12
push r13
push r14
mov [rbp+16], ecx ;m
; mov [rbp+24], edx ;n
mov [rbp+32], r8d ;p
mov [rbp+40], r9 ;&A
mov r10, [rbp+48] ;&B
mov r14, r10 ;r14 <- &B
mov r11, [rbp+56] ;&C
mov r12, r11 ;r12<-&C
mov rax, r8 ;rax <- p
shr rax, 2 ; p/4->nombre de paquets de colonnes de B à traiter
shl r8, 3 ;Taille en octets d'une ligne de B : 8 * p
mov r13, rdx ;r13 <- n
shr rdx, 2 ; n/4->nombre de paquets de colonnes de A à traiter
shl r13, 3 ;Taille en octets d'une ligne de A : 8 * n
boucle:
vmovapd ymm0, YMMWORD PTR [r11] ;On charge le contenu de C (précédemment calculé)
vmovapd ymm2, YMMWORD PTR [r9] ;On charge le contenu de A
vpermpd ymm1, ymm2, 00000000b ;On prend le premier nombre dans A
vfmadd231pd ymm0, ymm1, YMMWORD PTR [r10] ;On multiplie et additionne
add r10, r8 ;On prend la ligne suivante de B
vpermpd ymm1, ymm2, 01010101b ;On prend le second nombre dans A
vfmadd231pd ymm0, ymm1, YMMWORD PTR [r10] ;On multiplie et additionne
add r10, r8 ;On prend la ligne suivante de B
vpermpd ymm1, ymm2, 10101010b ;On prend le troisième nombre dans A
vfmadd231pd ymm0, ymm1, YMMWORD PTR [r10] ;On multiplie et additionne
add r10, r8 ;On prend la ligne suivante de B
vpermpd ymm1, ymm2, 11111111b ;On prend le quatrième nombre dans A
vfmadd231pd ymm0, ymm1, YMMWORD PTR [r10] ;On multiplie et additionne
vmovntpd YMMWORD PTR[r11], ymm0 ;On stocke C
add r11, r8 ;On prend la ligne suivante de C
add r9, r13 ;On prend la ligne suivante dans A
mov r10, r14 ;On réinitialise la position dans B
loop boucle ;On boucle sur m, le nombre de lignes de A
;reinit r9
mov r9, [rbp+40] ;On remet &A à sa valeur initiale
mov r10, r14
add r10, 32 ; on ajoute 4 colonnes à B
mov r14, r10 ; On stocke ces nouvelles 4 colonnes
mov r11, r12 ;On remet &C à sa valeur initiale
add r11, 32 ;On lui ajoute 4 colonnes
mov r12, r11 ;Qu'on stocke
mov ecx, [rbp+16] ;On réinitialise ecx à m
dec eax ; on décrémmente p/4
cmp eax, 0
ja boucle ;On boucle sur p, le nombre de colonnes de B
add r9, 32 ;On ajoute 4 colonnes à A
mov [rbp+40], r9
mov r10, [rbp+48] ;On réinitialise B
add r10, r8
add r10, r8
add r10, r8
add r10, r8
mov r14, r10
mov [rbp+48], r10 ;On réinitialise B
mov r11, [rbp+56] ;On réinitialise C
mov r12, r11
mov eax, [rbp+32] ;On réinitialise eax à p/4
shr eax, 2
dec rdx ;On passe aux colonnes de A suivantes.
cmp rdx, 0
ja boucle
pop r14
pop r13
pop r12
pop rbp
ret
crossProduct_asm ENDP
Attention, ce code ne fonctionne que sur des matrices alignées sur 32 octets, alors que malloc aligne sur du 16 octets. Il faut également que n et p soient des nombres divisibles par 4.
Le déterminant est un élément important en mathématiques. Il permet de savoir si des équations ont des solutions, ce genre de choses.
La méthode de calcul est longue et fastidieuse. Si longue et fastidieuse que pour les matrices de taille supérieure à, mettons, 10, on ne l'applique pas. On utilise d'autres techniques pour avoir un déterminant plus simple à calculer. Néanmoins, nous allons la coder comme exercice, car nous allons devoir faire appel à malloc() et free(), mais aussi à la récursivité.
La récursivité, c'est simplement quand une fonction s'appelle elle-même. La factorielle en est l'exemple typique : factorielle de 5 égale 5 fois factorielle de 4. En pseudo-code : fact(i) = i * fact(i-1). On met une condition d'arrêt, bien sûr. Ici, c'est si i == 1, fact(i) = 1;
PUBLIC lmatrice_asm
lmatrice_asm PROC
push rsi
push rdi
lea rsi, [r8 + rcx * 8]
lea r8, [r9 + rdx * 8]
dec rcx
mov r10, rcx
mov rax, rcx
mul rcx
mov rcx, rax
mov rdi, r9
boucle:
cmp rdi, r8
jne ecriture
add rsi, 8
lea r8, [r8 + r10 * 8]
ecriture:
movsq
loop boucle
pop rdi
pop rsi
ret
lmatrice_asm ENDP
PUBLIC det_asm
det_asm PROC
;double det_asm(int n, double mat[])
;{
; if (n == 1)
; {
; return mat[0];
; }
; double resultat = 0.;
; int k = n - 1;
; double signe;
; signe = 1.;
; double* lmat = new double[k * k];
; for (int i = 0; i < n; i++)
; {
; lmatrice_asm(n, i, mat, lmat);
; resultat = resultat + signe * mat[i] * det_asm(k, lmat);
; signe = -signe;
; }
; delete(lmat);
; return resultat;
;}
cmp rcx, 1
jne continue
vmovsd xmm0, QWORD PTR [rdx]
ret
continue:
mov QWORD PTR [rsp+24], rbp ; sauvegarde de rbp dans la stack shadow
push rdi
push r14
push r15
sub rsp, 80 ; 00000050H ; On laisse 80 octets libres sur la pile.
lea r15d, DWORD PTR [rcx-1] ; r15d <- k
movaps XMMWORD PTR [rsp+64], xmm5 ; On stocke xmm5 en haut de la pile
movsd xmm5, QWORD PTR __real@3ff0000000000000
mov edi, ecx ; edi <- n
mov eax, r15d ; eax <- k
mov rbp, rdx ; rdx <- l
movaps XMMWORD PTR [rsp+48], xmm7 ; On stocke xmm7 plus bas dans la pile
xorps xmm7, xmm7 ; xmm7 <- 0
lea r8, QWORD PTR [rax*8] ; r8 <- k * 8
mul r8 ; rax <- k * k * 8
mov rcx, rax ; rcx <- rax
call malloc
mov r14, rax ; r14 <- lmat
mov QWORD PTR [rsp+112], rbx; on sauve RBX car il n'est pas volatile
xor ebx, ebx ;RAZ de rbx
mov QWORD PTR [rsp+120], rsi
mov rsi, rbp
movaps XMMWORD PTR [rsp+32], xmm4; on sauvegarde xmm4
movsd xmm4, QWORD PTR __xmm@8000000000000000
; nop
boucle:
mov r9, r14 ; lmat
mov r8, rbp ; mat
mov edx, ebx ; l
mov ecx, edi ; k
call lmatrice_asm
mov rdx, r14 ;lmat
mov ecx, r15d ;k
call det_asm ; det_asm
movaps xmm1, xmm5 ;On met le signe dans xmm1
inc ebx
mulsd xmm1, QWORD PTR [rsi]; Qu'on multiplie par mat[i]
xorps xmm5, xmm4
add rsi, 8
; vfmadd231pd xmm7, xmm0, xmm1
mulsd xmm0, xmm1
addsd xmm7, xmm0
cmp ebx, edi
jl SHORT boucle
movaps xmm4, XMMWORD PTR [rsp+32]
mov rsi, QWORD PTR [rsp+120]
mov rcx, r14
call free
movsd xmm0, xmm7
mov rbp, QWORD PTR [rsp+128]
mov rbx, QWORD PTR [rsp+112]
movaps xmm5, XMMWORD PTR [rsp+64]
movaps xmm7, XMMWORD PTR [rsp+48]
add rsp, 80 ; 00000050H On détruit la pile
pop r15
pop r14
pop rdi
ret
det_asm ENDP
J'obtiens, principalement pour de belles matrices, un temps d'exécution qui peut être divisé par deux pour le produit matriciel. La différence est moins flagrante pour le déterminant.