| <-Précédent | Retour à l'accueil | Contact : etienne"at"les-sauvages.fr |
Pour tous ceux qui sont passés sur les bancs d'une école d'ingénieur ou de Math Sup', l'une des premières choses que l'on voit dans le cours d'algèbre linéaire est la Décomposition PLU. P comme Permutation, L comme Lower (basse) et U comme Upper(haute). Il s'agit de décomposer la matrice en 3 : une matrice de permutation, qui ne fait que permuter des lignes, une matrice triangulaire vers le bas et une triangulaire vers le haut. Je vous conseille l'article Wikipedia
Je rappelle à mes petits camarades que la matrice qu'on se propose de décomposer et que nous appellerons A dans la suite de ce chapitre, je rappelle donc que cette matrice se doit d'être carrée.
Wikipédia étant mon ami, il me fournit un code C. Code que nous allons légèrement modifier pour obtenir ceci :
/* INPUT: A - array of pointers to rows of a square matrix having dimension N
* Tol - small tolerance number to detect failure when the matrix is near degenerate
* OUTPUT: Matrix A is changed, it contains a copy of both matrices L-E and U as A=(L-E)+U such that P*A=L*U.
* The permutation matrix is not stored as a matrix, but in an integer vector P of size N+1
* containing column indexes where the permutation matrix has "1". The last element P[N]=S+N,
* where S is the number of row exchanges needed for determinant computation, det(P)=(-1)^S
*/
int LUPDecompose(double* A, int N, double Tol, int* P)
{
for (auto i = 0; i <= N; i++)
P[i] = i; //Unit permutation matrix, P[N] initialized with N
for (auto i = 0; i < N; i++) {
double maxA = 0.0;
auto imax = i;
for (auto k = i; k < N; k++)
{
auto absA = fabs(A[k * N + i]);
if (absA > maxA) {
maxA = absA;
imax = k;
}
}
if (maxA < Tol) return 0; //failure, matrix is degenerate
if (imax != i) {
//pivoting P
auto j = P[i];
P[i] = P[imax];
P[imax] = j;
//pivoting rows of A
for (auto k = 0; k < N; k++)
{
auto buf = A[i*N + k];
A[i*N + k] = A[imax*N + k];
A[imax*N + k] = buf;
}
//counting pivots starting from N (for determinant)
P[N]++;
}
for (auto j = i + 1; j < N; j++) {
A[j * N + i] /= A[i * N + i]; //On divise par l'élément diagonal
for (auto k = i + 1; k < N; k++)
A[j * N + k] -= A[j * N + i] * A[i * N + k];//On enlève
}
}
return 1; //decomposition done
}
}
Essentiellement, la modification porte sur le passage d'une matrice en double pointeur à une matrice à pointeur simple, qui nous permettra d'allouer la mémoire en une seule opération, et de libérer cette même mémoire en, tiens comme c'est curieux, une seule opération aussi. Cela nous assure également que la matrice sera tockée en un seul bloc mémoire, ce qui nous permettra de faire de l'arithmétique des pointeurs.
Quelques commentaires sur ce code :
On peut diviser une itération sur i en deux parties : d'abord, la recherche de la meilleure permutation, et ensuite, le calcul des éléments de L et U. La recherche de la meilleure permutation s'effectue aussi en deux parties : d'abord, la recherche du maximum de la colonne i à partir de l'index i, et ensuite, la permutation. Nous allons dans un premier temps mettre ces parties dans des fonctions. En assembleur MASM, cela donne :
PUBLIC fabs_max_colonne_asm fabs_max_colonne_asm PROC ;rcx <- n ;rdx <- i ;r8 <- A xorps xmm0, xmm0 movsd xmm3, QWORD PTR __xmm@7fffffffffffffff7fffffffffffffff mov rax, rdx ;RAX <- i imul eax, ecx ;RAX <- i * n shl rax, 3 ;RAX <- i * n * sizeof(double) lea rax, [rax + rdx * 8] lea r8, [r8 + rax] ;RAX <- A[i, i] mov r9, rcx shl r9, 3 sub rcx, rdx mov r10, rcx boucle: movsd xmm1, QWORD PTR[r8] xorps xmm2, xmm2 subsd xmm2, xmm1 ; movsd xmm2, xmm1 ; andpd xmm2, xmm3 maxsd xmm1, xmm2 comisd xmm0, xmm1 jb nouveauMax ; jg ne fonctionne pas finBoucle: add r8, r9 loop boucle ret nouveauMax: movapd xmm0, xmm1 mov rax, r10 sub rax, rcx add rax, rdx jmp finBoucle fabs_max_colonne_asm ENDP ;Permute two lines of a matrix PUBLIC permute_lines_asm permute_lines_asm PROC ;rcx <- i ;rdx <- j ;r8 <- n ;r9 <- A push rsi push rdi mov rsi, rcx ; rsi<- i * n* 8 imul rsi, r8 shl rsi, 3 add rsi, r9 mov rdi, rdx ; rdi<- j * n* 8 imul rdi, r8 shl rdi, 3 add rdi, r9 boucle: lodsq xchg rsi, rdi mov r10, rax lodsq sub rdi, 8 stosq xchg rsi, rdi sub rdi, 8 mov rax, r10 stosq xchg rsi, rdi dec r8 jnz boucle pop rdi pop rsi ret permute_lines_asm ENDP
Quelques commentaires : tout d'abord, la fonction fabs_max_colonne_asm retourne deux valeurs. Pour l'appeler en C, la valeur de retour sera déterminée par ce qu'on donnera en signature de la fonction : int fabs_max_colonne_asm retournera l'index du maximum, tandis que double fabs_max_colonne_asm retournera la valeur du minimum.
Ensuite, la fonction valeur absolue d'un double n'existe pas. On utilise le fait que la valeur absolue d'un nombre est le maximum entre ce nombre et son opposé. C'est la méthode utilisée ici.
PUBLIC PLUDecompose_asm
PLUDecompose_asm PROC
; for (auto i = 0; i < N; i++) {
; for (auto j = i + 1; j < N; j++) {
; A[j * N + i] /= A[i * N + i]; //On divise par l'élément diagonal
; for (auto k = i + 1; k < N; k++)
; A[j * N + k] -= A[j * N + i] * A[i * N + k];//On enlève le pivot * l'élément courant
; }
;RCX <- n, RDX <- A
;r8 <- k
;xmm0 <- A[j * N + i]
;r9 <- i
;r10 <- j
xor r9, r9
push r12
push r13
push r14
push r15
mov r13, rcx ;R13 <- n
mov r12, rdx ;R12 <- A
boucle_i:
mov r11, r9 ;R11 <- i
mov rcx, r13
mov rdx, r9 ;rdx <- i
mov r8, r12 ;R8 <- A
call fabs_max_colonne_asm
mov rcx, rax ;RCX <- i[max]
mov rdx, r11 ;RDX <- i
mov r8, r13 ;R8 <- n
mov r9, r12 ;R9 <- A
call permute_lines_asm
mov rdx, r12 ;RDX <- A
mov r9, r11 ;R9 <- i
mov rcx, r13 ;RCX <- n
mov r10, r9 ;R10 <- i
inc r10 ;R10 <- i + 1, ie. R10 <- j
cmp r10, rcx
jae Fin
mov r14, r9 ;R14 <- i
imul r14d, ecx ;R14 <- i * n
shl r14, 3 ;R14 <- i * n * sizeof(double)
boucle_j:
mov r15, r10 ;R15 <- j
imul r15d, ecx ;R15 <- j * n
shl r15, 3 ;R15 <- j * n * sizeof(double)
lea rax, [r15 + r9 * 8] ;RAX <- (j * n + i) * sizeof(double)
add rax, rdx ;RAX <- A[j * n + i]
movsd xmm0, QWORD PTR [rax];XMM0 <- *A[j * n + i]
lea r8, [r14 + r9 * 8] ;R8 <- (i * n + i) * sizeof(double)
add r8, rdx ;R8 <- A[i * n + i]
movsd xmm1, QWORD PTR [r8];XMM1 <- *A[i * n + i]
divsd xmm0, xmm1 ;XMM0 <- A[j * n + i]/A[i * n + i]
movsd QWORD PTR [rax], xmm0; A[j * N + i] /= A[i * N + i]
mov r8, r9 ;R8 <- i
inc r8 ;R8 <- i + 1, ie k
movsd xmm2, xmm0
boucle_k:
lea rax, [r14 + r8 * 8] ;RAX <- (i * n + k) * sizeof(double)
add rax, rdx ;RAX <- A[i * n + k]
movsd xmm1, QWORD PTR [rax];XMM1 <- *A[i * n + k]
mulsd xmm0, xmm1 ;XMM0 <- A[j * N + i] * A[i * N + k]
lea rax, [r15 + r8*8] ;RAX <- (j * n + k) * sizeof(double)
add rax, rdx ;RAX <- A[j * n + k]
movsd xmm1, QWORD PTR [rax];XMM1 <- *A[j * n + k]
subsd xmm1, xmm0 ;XMM1 <- A[j * n + k] - A[j * N + i] * A[i * N + k]
movsd QWORD PTR [rax], xmm1;A[j * N + k] -= A[j * N + i] * A[i * N + k]
movsd xmm0, xmm2
inc r8
cmp r8, rcx
jb boucle_k
inc r10
cmp r10, rcx
jb boucle_j
inc r9
cmp r9, r13
jb boucle_i
Fin:
pop r15
pop r14
pop r13
pop r12
ret 0
PLUDecompose_asm ENDP
Ben déjà, ça marche. La matrice A est modifiée pour contenir à la fois L et U, et la matrice P, essentiellement creuse, est ramenée à un vecteur contenant les index permutés. C'est à ma connaissance la première implémentation de la décomposition PLU en assembleur. Ensuite, il nous faut optimiser cette fonction, parce qu'elle est pour l'instant plus lente que la version générée par le compilateur.
Qu'optimiser ? Hé bien, principalement la boucle la plus interne, celle qui va être appelée le plus grand nombre de fois, la boucle sur k. Il faut d'abord enlever toute instruction superflue. Le MOVSD XMM0, XMM2 est en trop.
Ensuite, le MULSD peut prendre comme source une adresse mémoire. Donc supprimons le MOVSD précédent.
Pour continuer, on remarque qu'on utilise le même registre rax pour stocker deux valeurs différentes. Cela empêche le processeur de charger en même temps les 2 valeurs, créant une chaîne de dépendance. Donc, on va plutôt utiliser les registres RSI et RDI pour stocker ces deux valeurs.
Enfin, nous remarquons que cette boucle peut être vectorisée. Nous allons en créer une version qui s'exécute s'il reste au moins 4 itérations à faire, et une autre pour moins de 3 opérations.
On peut également intégrer le code des deux fonctions appelées dans la boucle i pour éviter les CALL quelque peu consommateurs de ressources.
On obtient ceci :
PUBLIC PLUDecompose_asm
PLUDecompose_asm PROC
; for (auto i = 0; i < N; i++) {
; for (auto j = i + 1; j < N; j++) {
; A[j * N + i] /= A[i * N + i]; //On divise par l'élément diagonal
; for (auto k = i + 1; k < N; k++)
; A[j * N + k] -= A[j * N + i] * A[i * N + k];//On enlève le pivot * l'élément courant
; }
;RCX <- n, RDX <- A
;r8 <- k, r9<- P
;xmm0 <- A[j * N + i]
;rbx <- i
;r10 <- j
push r12
push r13
push r14
push r15
push rbx
push rsi
push rdi
mov r13, rcx ;R13 <- n
mov r12, rdx ;R12 <- A
; for (auto i = 0; i <= N; i++)
; P[i] = i; //Unit permutation matrix, P[N] initialized with N
boucle_P:
lea rax, [r9 + rcx * 4]
mov DWORD PTR[rax], ecx
loop boucle_P
; lea rax, [r9 + rcx * 4]
mov DWORD PTR[r9], 0
mov rcx, r13
xor rbx, rbx
; movsd xmm3, QWORD PTR __xmm@7fffffffffffffff7fffffffffffffff
boucle_i:
mov r11, rbx ;R11 <- i
mov rcx, r13
mov rdx, r11 ;rdx <- i
mov r8, r12 ;R8 <- A
;rcx <- n
;rdx <- i
;r8 <- A
xorps xmm0, xmm0
mov rax, r11 ;RAX <- i
imul eax, ecx ;RAX <- i * n
shl rax, 3 ;RAX <- i * n * sizeof(double)
lea rax, [rax + rdx * 8]
lea r8, [r8 + rax] ;RAX <- A[i, i]
mov rbx, rcx
shl rbx, 3
sub rcx, rdx
mov r10, rcx
boucle_n:
movsd xmm1, QWORD PTR[r8]
xorps xmm2, xmm2
subsd xmm2, xmm1
; movsd xmm2, xmm1
; andpd xmm2, xmm3
maxsd xmm1, xmm2
comisd xmm0, xmm1
jb nouveauMax ; jg ne fonctionne pas
finBoucle:
add r8, rbx
loop boucle_n
cmp rax, r11
je post_perm
;Permutation des lignes
mov ecx, eax ;RCX <- i[max]
mov r8d, r13d ;R8 <- n
;rcx <- i
;r11 <- j
;r8 <- n
;r12 <- A
;r9 <- P
; //pivoting P
; auto j = P[i];
; P[i] = P[imax];
; P[imax] = j;
lea r14, [r9 + r8 * 4] ;r14 <- P[N]
mov r15d, DWORD PTR[r14] ;r15 <- *P[N]
inc r15d ;P[N]++
mov DWORD PTR[r14], r15d ;stocke
mov esi, ecx ; rsi<- i * n* 8
imul esi, r8d
shl esi, 3
add rsi, r12
mov edi, r11d ; rdi<- j * n* 8
imul edi, r8d
shl edi, 3
add rdi, r12
boucle:
lodsq ;Charge dans RAX
xchg rsi, rdi
mov r10, rax
lodsq
sub rdi, 8
stosq
xchg rsi, rdi
sub rdi, 8
mov rax, r10
stosq
xchg rsi, rdi
dec r8
jnz boucle
lea rdi, [r9 + rcx * 4] ;rdi <- P[i]
lea rsi, [r9 + r11 * 4] ;rsi <- P[j]
mov r14d, DWORD PTR[rsi] ;r14 <- *P[j]
mov r15d, DWORD PTR[rdi] ;r15 <- *P[i]
mov DWORD PTR[rsi], r15d ;P[j] <- P[i]
mov DWORD PTR[rdi], r14d ;P[i] <- P[j]
post_perm:
mov ebx, r11d ;RBX <- i
mov ecx, r13d ;RCX <- n
mov r10d, r11d ;R10 <- i
inc r10d ;R10 <- i + 1, ie. R10 <- j
cmp r10d, ecx
jae Fin
mov r14d, r11d ;R14 <- i
imul r14d, ecx ;R14 <- i * n
shl r14d, 3 ;R14 <- i * n * sizeof(double)
boucle_j:
mov r15d, r10d ;R15 <- j
imul r15d, ecx ;R15 <- j * n
shl r15d, 3 ;R15 <- j * n * sizeof(double)
lea eax, [r15d + r11d * 8] ;RAX <- (j * n + i) * sizeof(double)
add rax, r12 ;RAX <- A[j * n + i]
movsd xmm0, QWORD PTR [rax] ;XMM0 <- *A[j * n + i]
lea r8d, [r14d + r11d * 8] ;R8 <- (i * n + i) * sizeof(double)
add r8, r12 ;R8 <- A[i * n + i]
movsd xmm1, QWORD PTR [r8] ;XMM1 <- *A[i * n + i]
divsd xmm0, xmm1 ;XMM0 <- A[j * n + i]/A[i * n + i]
movsd QWORD PTR [rax], xmm0 ; A[j * N + i] /= A[i * N + i]
mov r8d, r11d ;R8 <- i
inc r8d ;R8 <- i + 1, ie k
vpermpd ymm0, ymm0, 0
boucle_k:
mov eax, ecx
sub eax, r8d
cmp eax, 3
jbe boucle_k_pair
lea edi, [r15d + r8d*8] ;RAX <- (j * n + k) * sizeof(double)
add rdi, r12 ;RAX <- A[j * n + k]
vmovupd ymm2, YMMWORD PTR [rdi] ;YMM2 <- *A[j * n + k]
lea esi, [r14d + r8d * 8] ;RAX <- (i * n + k) * sizeof(double)
add rsi, r12 ;RAX <- A[i * n + k]
vmulpd ymm1, ymm0, YMMWORD PTR [rsi] ;YMM1 <- A[j * N + i] * A[i * N + k]
vsubpd ymm2, ymm2, ymm1 ;YMM2 <- A[j * n + k] - A[j * N + i] * A[i * N + k]
; vfnmadd231pd ymm2, ymm0, YMMWORD PTR [rsi]
vmovupd YMMWORD PTR [rdi], ymm2 ;A[j * N + k] -= A[j * N + i] * A[i * N + k]
add r8d, 4
fin_boucle_k:
cmp r8d, ecx
jb boucle_k
inc r10d
cmp r10d, ecx
jb boucle_j
inc ebx
cmp ebx, r13d
jb boucle_i
Fin:
pop rdi
pop rsi
pop rbx
pop r15
pop r14
pop r13
pop r12
ret 0
nouveauMax:
movapd xmm0, xmm1
mov eax, r10d
sub eax, ecx
add eax, edx
jmp finBoucle
boucle_k_impair:
lea esi, [r14d + r8d * 8] ;RAX <- (i * n + k) * sizeof(double)
add rsi, r12 ;RAX <- A[i * n + k]
vmulsd xmm1, xmm0, QWORD PTR [rsi] ;XMM1 <- A[j * N + i] * A[i * N + k]
lea edi, [r15d + r8d * 8] ;RAX <- (j * n + k) * sizeof(double)
add rdi, r12 ;RAX <- A[j * n + k]
movsd xmm2, QWORD PTR [rdi] ;XMM2 <- *A[j * n + k]
subsd xmm2, xmm1 ;XMM2 <- A[j * n + k] - A[j * N + i] * A[i * N + k]
movsd QWORD PTR [rdi], xmm2 ;A[j * N + k] -= A[j * N + i] * A[i * N + k]
inc r8d
jmp fin_boucle_k
boucle_k_pair:
cmp eax, 1
jbe boucle_k_impair
lea esi, [r14d + r8d * 8] ;RAX <- (i * n + k) * sizeof(double)
add rsi, r12 ;RAX <- A[i * n + k]
vmulpd xmm1, xmm0, XMMWORD PTR [rsi];YMM1 <- A[j * N + i] * A[i * N + k]
lea edi, [r15d + r8d*8] ;RAX <- (j * n + k) * sizeof(double)
add rdi, r12 ;RAX <- A[j * n + k]
movupd xmm2, XMMWORD PTR [rdi] ;YMM2 <- *A[j * n + k]
subpd xmm2, xmm1 ;YMM2 <- A[j * n + k] - A[j * N + i] * A[i * N + k]
vmovupd XMMWORD PTR [rdi], xmm2 ;A[j * N + k] -= A[j * N + i] * A[i * N + k]
add r8d, 2
jmp fin_boucle_k
PLUDecompose_asm ENDP
Et comme