| Autor: | Wojciech Muła |
|---|---|
| Dodany: | 13.03.2002 |
| Aktualizacja: | 15.09.2007 |
Treść
Po opis rozkazów odsyłam do dokumentacji procesora. Rozkazy SSE umożliwiają obliczenia na liczbach zmiennoprzecinkowych pojedynczej precyzji (4 bajtowych). Przy korzystaniu z SSE bardzo istotne jest wyrównanie danych — większość rozkazów domyślnie przyjmuje, że dane są na granicy 16 bajtów — gdy tak nie jest generowany jest wyjątek #GP.
Procesory 32-bitowe posiadają osiem, natomiast 64-bitowe szesnaście rejestrów SSE, każdy o rozmiarze 128 bitów (4 x liczba zmiennoprzecinkowa). Ich nazwy w asemblerze to: xmm0, xmm1, ... xmm15.
Uwaga: w komentarzach „?” oznacza nieważny, używany jest do oznaczenia nieistotnych wartości w rejestrach SSE.
Iloczyn skalarny (ang. dot product) dwóch wektorów 3D V1 = [x1, y1, y1] i V2 = [x2, y2, y2] dany jest wzorem:
DP = x1 ⋅ x2 + y1 ⋅ y2 + z1 ⋅ z2
void sse_dot(float vec1[4], float vec2[4], float* result) {
__asm__ volatile (
"movups (%0), %%xmm0 \n" // load vec1: |w1|z1|y1|x1|
"movups (%1), %%xmm1 \n" // load vec2: |w2|z2|y2|x2|
" \n"
"mulps %%xmm1, %%xmm0 \n" // xmm0 := |w1*w2|z1*z2|y1*y2|x1*x2|
"movhlps %%xmm0, %%xmm1 \n" // xmm1 := | . | . |w1*w2|z1*z2|
"addps %%xmm0, %%xmm1 \n" // xmm1 := | . | . | w+y | z+x |
"movaps %%xmm1, %%xmm0 \n" // save xmm1
"shufps $0x01, %%xmm1, %%xmm1 \n" // xmm1 := | . | . | . | w+y |
"addss %%xmm1, %%xmm0 \n" // xmm0[0] := dot product
"movss %%xmm0, (%2) \n"
:
: "r" (vec1), "r" (vec2), "r" (result));
}
Żeby jednak w 100% wykorzystać moc obliczeniową, należy liczyć równocześnie 4 iloczyny skalarne.
void sse_dot4(float* v1x4, float* v2x4, float* results) {
__asm__ volatile (
// load 4 vectors from v1
"movups (%0), %%xmm0 \n"
"movups 16(%0), %%xmm1 \n"
"movups 32(%0), %%xmm2 \n"
"movups 48(%0), %%xmm3 \n"
// load 4 vectors from v2
"movups (%1), %%xmm4 \n"
"movups 16(%1), %%xmm5 \n"
"movups 32(%1), %%xmm6 \n"
"movups 48(%1), %%xmm7 \n"
// perfom parallel multiplications
"mulps %%xmm4, %%xmm0 \n" // xmm0 := |A3|A2|A1|A0|
"mulps %%xmm5, %%xmm1 \n" // xmm1 := |B3|B2|B1|B0|
"mulps %%xmm6, %%xmm2 \n" // xmm2 := |C3|C2|C1|C0|
"mulps %%xmm7, %%xmm3 \n" // xmm3 := |D3|D2|D1|D0|
// (xmm4-xmm7 are free)
// perfom additions
"movaps %%xmm0, %%xmm4 \n"
"unpcklps %%xmm1, %%xmm0 \n" // xmm0 := |B1|A1|B0|A0|
"unpckhps %%xmm1, %%xmm4 \n" // xmm4 := |B3|A3|B2|A2| (xmm1 is free)
"movaps %%xmm2, %%xmm1 \n"
"unpcklps %%xmm3, %%xmm2 \n" // xmm2 := |D1|C1|D0|C0|
"unpckhps %%xmm3, %%xmm1 \n" // xmm1 := |D3|C3|D2|C2| (xmm3 is free)
"addps %%xmm4, %%xmm0 \n" // xmm0 := | B13 | A13 | B02 | A02 | (xmm1 is free)
"addps %%xmm1, %%xmm2 \n" // xmm2 := | D13 | C13 | D02 | C02 | (xmm3 is free)
"movaps %%xmm0, %%xmm1 \n"
"shufps $0b01000100, %%xmm2, %%xmm0 \n" // xmm0 := | D02 | C02 | B02 | A02 |
"shufps $0b11101110, %%xmm2, %%xmm1 \n" // xmm1 := | D13 | C13 | B13 | A13 |
"addps %%xmm1, %%xmm0 \n" // xmm0 := |D0123|C0123|B0123|A0123|
"movups %%xmm0, (%2) \n"
:
: "r" (v1x4), "r" (v2x4), "r" (results)
);
}
Wynikiem iloczynu wektorowego (ang. cross product) dwóch wektorów jest wektor prostopadły (normalny) do płaszczyzny na której znajdują się te wektory. Iloczyn wektorowy liczy się z wyznacznika:
|i j k|
|a b c| = i(b*z-c*y) - j(a*z-c*x) + k(a*y-b*x)
|x y z|
= [b*z - c*y, c*x - a*z, a*y - b*x] =
= [b*z, c*x, a*y] - [c*y, a*z, b*x]
gdzie i, j, k — wersory osi, w uładzie kartezjańskim i = [1, 0, 0], j = [0, 1, 0], oraz k = [0, 0, 1].
void sse_cross_prod(float v1[4], float v2[4], float result[4]) {
__asm__ volatile (
// load vectors index -> 0 1 2 3
"movups (%0), %%xmm0 \n" // xmm0 := | a | b | c | . |
"movups (%1), %%xmm1 \n" // xmm1 := | x | y | z | . |
// clone
"movaps %%xmm0, %%xmm2 \n"
"movaps %%xmm1, %%xmm3 \n"
// shuffle
"shufps $0b00001001, %%xmm0, %%xmm0 \n" // xmm0 := | b | c | a | . |
"shufps $0b00010010, %%xmm1, %%xmm1 \n" // xmm1 := | z | x | y | . |
"shufps $0b00010010, %%xmm2, %%xmm2 \n" // xmm2 := | c | a | b | . |
"shufps $0b00001001, %%xmm3, %%xmm3 \n" // xmm3 := | y | z | x | . |
// multiply
"mulps %%xmm1, %%xmm0 \n" // xmm0 := |b*z|c*x|a*y| . |
"mulps %%xmm3, %%xmm2 \n" // xmm2 := |c*y|a*z|b*x| . |
// and finally subtract
"subps %%xmm2, %%xmm0 \n" // xmm0 := |b*z - c*y|c*x - a*z|a*y - b*x| ... |
"movups %%xmm0, (%2) \n"
:
: "r" (v1), "r" (v2), "r" (result)
);
}
Wynikiem normalizacja wektora n-wymiarowego jest wektor o długości 1 i o takich samych cosinusach kierunkowych jak wektor wejściowy. Wzór jest następujący:
V' = V/|V|
; 11.03.2002
;
; edi -> wektor wejściowy
movups xmm0, [edi] ; załaduj wektor
movaps xmm1, xmm0 ; xmm1 = | 0 | z | y | x |
mulps xmm1, xmm0 ; xmm1 = | 0 | z*z | y*y | x*x |
movaps xmm2, xmm1
movaps xmm3, xmm1
shufps xmm2, xmm1, 0fdh ; xmm2 = | 0 | 0 | 0 | y^2 |
shufps xmm3, xmm1, 0feh ; xmm3 = | 0 | 0 | 0 | z^2 |
addps xmm1, xmm2 ;
addps xmm1, xmm3 ; xmm1 = | 0 | nu | nu | x^2+y^2+z^2 |
sqrtss xmm1 ; xmm1 = | 0 | nu | nu | sqrt(x^2+y^2+z^2) |
; len
shufps xmm1, xmm1, 000h ; xmm1 = | len | len | len | len |
divps xmm0, xmm1 ; xmm0 = | 0 |z/len|y/len|x/len| -- wynik
Obrót punktu o kąt fi wokół środka układu współrzędnych dany jest wzorem:
x' = xcos(fi) − ysin(fi) y' = xsin(fi) + ycos(fi)
; 11.03.2002
segment .data
align 16
sincos:
dd ? ; sinus(fi)
dd ? ; cosinus(fi)
vector:
dd ? ; x
dd ? ; y
negate:
dd 0
dd 0
dd 0x80000000
dd 0
segment .text
movlps xmm0, [vector] ; xmm0 = | ? | ? | y | x |
movlps xmm1, [sincos] ; xmm1 = | ? | ? | c | s |
unpcklps xmm0, xmm0 ; xmm0 = | y | y | x | x |
shufps xmm1, xmm1, 044h ; xmm1 = | c | s | c | s |
xorps xmm1, [negate] ; xmm1 = | c |-s | c | s |
mulps xmm0, xmm1 ; xmm0 = | y*c | -y*s | x*c | x*s |
movaps xmm1, xmm0
shufps xmm1, xmm1, 00eh ; xmm1 = | nu | nu |-y*s | y*c |
addps xmm0, xmm1 ; xmm0 = | nu | nu |xc-ys|xs+yc| -- wynik
Wynikiem mnożenia liczb zespolonych (a + ib) i (c + id) jest liczba zespolona (ac − bd) + i(bc + ad).
; 13-03-02
; edi -> (a+ib) (typ float[2])
; esi -> (c+id) (typ float[2])
; edx -> wynik
segment .data
__negate:
dd 0
dd 0
dd 0
dd 0x80000000
segment .text
movlps xmm0, [edi] ; xmm0 = | nu | nu | b | a |
movlps xmm1, [esi] ; xmm1 = | nu | nu | d | c |
shufps xmm0, xmm0, 044h ; xmm0 = | b | a | b | a |
unpcklps xmm1, xmm1 ; xmm1 = | d | d | c | c |
xorps xmm1, [__negate] ; xmm1 = | -d | d | c | c |
mulps xmm0, xmm1 ; xmm0 = | -bd | ad | bc | ac |
movaps xmm1, xmm0
shufps xmm1, xmm1, 00bh ; xmm1 = | nu | nu | ad | -bd |
addps xmm0, xmm1 ; xmm0 = | nu | nu |ad+bc|ac-bd|
movlps [edx], xmm0
Mnożenie macierzy 4x4 przez wektor 4x1:
[a0, a1, a2, a3] [x] [a0*x + a1*y + a*z + a3] [b0, b1, b2, b3] [y] [b0*x + b1*y + b*z + b3] [ ] * [ ] = [ ] [c0, c1, c2, c3] [z] [c0*x + c1*y + c*z + c3] [d0, d1, d2, d3] [w] [d0*x + d1*y + d*z + d3]
void sse_matvec_mult(float mat[4*4], float vec[4], float result[4]) {
__asm__ volatile (
// load vector
"movups (%1), %%xmm4 \n" // xmm4 := [x, y, w, z] = V
// load matrix
"movups 0x00(%0), %%xmm0 \n" // xmm0 := [a0, a1, a2, a3] = A
"movups 0x10(%0), %%xmm1 \n" // xmm1 := [b0, b1, b2, b3] = B
"movups 0x20(%0), %%xmm2 \n" // xmm2 := [c0, c1, c2, c3] = C
"movups 0x30(%0), %%xmm3 \n" // xmm3 := [d0, d1, d2, d3] = D
// multiply all rows by vector
"mulps %%xmm4, %%xmm0 \n" // xmm0 := A*V
"mulps %%xmm4, %%xmm1 \n" // xmm1 := B*V
"mulps %%xmm4, %%xmm2 \n" // xmm2 := C*V
"mulps %%xmm4, %%xmm3 \n" // xmm3 := D*V
// finally calc dot products
"movaps %%xmm0, %%xmm4 \n" // save A
"unpcklps %%xmm1, %%xmm0 \n" // xmm0 := [a0, b0, a1, b1]
"unpckhps %%xmm1, %%xmm4 \n" // xmm4 := [a2, b2, a3, b3] (xmm1 is free)
" \n"
"movaps %%xmm2, %%xmm1 \n" // save C
"unpcklps %%xmm3, %%xmm2 \n" // xmm2 := [c0, d0, c1, d1]
"unpckhps %%xmm3, %%xmm1 \n" // xmm1 := [c2, d2, c3, d3] (xmm3 is free)
" \n"
"addps %%xmm4, %%xmm0 \n" // xmm0 := [a0+a2, b0+b2, a1+a3, b1+b3] (xmm4 is free)
"addps %%xmm2, %%xmm1 \n" // xmm1 := [c0+c2, bd+d2, c1+c3, d1+d3] (xmm2 is free)
" \n"
"movaps %%xmm0, %%xmm4 \n"
"shufps $0b01000100, %%xmm1, %%xmm0 \n" // xmm0 := [a02, b02, c02, d02 ]
"shufps $0b11101110, %%xmm1, %%xmm4 \n" // xmm4 := [a13, b13, c13, d13 ]
" \n"
"addps %%xmm4, %%xmm0 \n" // xmm0 := [a0123, b0123, c0123, d0123]
"movups %%xmm0, (%2) \n"
"0: \n"
:
: "r" (mat), "r" (vec), "r" (result)
: "memory"
);
}
Mnożenie wektora 1x4 przez macierz 4x4:
[a0, a1, a2, a3] [x*a0 + y*b0 + z*c0 + w*d0]
[b0, b1, b2, b3] [x*a1 + y*b1 + z*c1 + w*d1]
[x, y, z, w] * [ ] = [ ]
[c0, c1, c2, c3] [x*a2 + y*b2 + z*c2 + w*d2]
[d0, d1, d2, d3] [x*a3 + y*b3 + z*c3 + w*d3]
void sse_vecmat_mult(float mat[4*4], float vec[4], float result[4]) {
__asm__ volatile (
// load vector
"movups (%1), %%xmm4 \n" // xmm4 := [x, y, z, w] = V
"movaps %%xmm4, %%xmm5 \n"
"movaps %%xmm4, %%xmm6 \n"
"movaps %%xmm4, %%xmm7 \n"
// load matrix
"movups 0x00(%0), %%xmm0 \n" // xmm0 := [a0, a1, a2, a3] = A
"movups 0x10(%0), %%xmm1 \n" // xmm1 := [b0, b1, b2, b3] = B
"movups 0x20(%0), %%xmm2 \n" // xmm2 := [c0, c1, c2, c3] = C
"movups 0x30(%0), %%xmm3 \n" // xmm3 := [d0, d1, d2, d3] = D
// populate coords in vectors
"shufps $0b00000000, %%xmm4, %%xmm4 \n" // xmm4 := [x, x, x, x]
"shufps $0b01010101, %%xmm5, %%xmm5 \n" // xmm5 := [y, y, y, y]
"shufps $0b10101010, %%xmm6, %%xmm6 \n" // xmm6 := [z, z, z, z]
"shufps $0b11111111, %%xmm7, %%xmm7 \n" // xmm7 := [w, w, w, w]
// multiply all rows by vector
"mulps %%xmm4, %%xmm0 \n" // xmm0 := [x*a0, x*a1, x*a2, x*a3]
"mulps %%xmm5, %%xmm1 \n" // xmm1 := [y*b0, y*b1, y*b2, y*b3]
"mulps %%xmm6, %%xmm2 \n" // xmm2 := [z*c0, z*c1, z*c2, z*c3]
"mulps %%xmm7, %%xmm3 \n" // xmm3 := [w*d0, w*d1, w*d2, w*d3]
// finally calc dot products
"addps %%xmm1, %%xmm0 \n"
"addps %%xmm3, %%xmm2 \n"
"addps %%xmm2, %%xmm0 \n"
// save result
"movups %%xmm0, (%2) \n"
:
: "r" (mat), "r" (vec), "r" (result)
: "memory"
);
}
Transpozycja macierzy to — z technicznego punktu widzenia — zamiana wierszy z wektorami, tj.:
[a0 a1 a2 a3] [a0 b0 c0 d0]
[b0 b1 b2 b3] > [a1 b1 c1 d1]
[d0 d1 d2 d3] > [a2 b2 c2 d2]
[c1 c1 c2 c3] [a3 b3 c3 d3]
M M^T
void sse_transpose(float mat[4*4]) {
__asm__ volatile (
// 0 1 2 3 <- index
"movups 0x00(%0), %%xmm0 \n" // xmm0 := a0 a1 a2 a3
"movups 0x10(%0), %%xmm1 \n" // xmm1 := b0 b1 b2 b3
"movups 0x20(%0), %%xmm2 \n" // xmm2 := c0 c1 c2 c3
"movups 0x30(%0), %%xmm3 \n" // xmm3 := d0 d1 d2 d3
"movaps %%xmm0, %%xmm4 \n"
"unpcklps %%xmm2, %%xmm0 \n" // xmm0 := a0 c0 a1 c1
"unpckhps %%xmm2, %%xmm4 \n" // xmm4 := a2 c2 a3 c3
"movaps %%xmm1, %%xmm2 \n"
"unpcklps %%xmm3, %%xmm1 \n" // xmm1 := b0 d0 b1 d1
"unpckhps %%xmm3, %%xmm2 \n" // xmm2 := b2 d2 b3 d3
"movaps %%xmm0, %%xmm3 \n"
"unpcklps %%xmm1, %%xmm0 \n" // xmm0 := a0 b0 c0 d0
"unpckhps %%xmm1, %%xmm3 \n" // xmm3 := a1 b1 c1 d1
"movaps %%xmm4, %%xmm1 \n"
"unpcklps %%xmm2, %%xmm1 \n" // xmm1 := a2 b2 c2 d2
"unpckhps %%xmm2, %%xmm4 \n" // xmm4 := a3 b3 c3 d3
"movups %%xmm0, 0x00(%0) \n"
"movups %%xmm3, 0x10(%0) \n"
"movups %%xmm1, 0x20(%0) \n"
"movups %%xmm4, 0x30(%0) \n"
: "=r" (mat)
: "0" (mat)
: "memory"
);
}
Poniższa procedura jest szybsza o ok. 2,5 raza od procedury napisanej w języku C i skompilowanej przez GCC z opcjami -O2 -ffast-math.
void sse_matmat_mult(float mat1[4*4], float mat2[4*4], float mat3[4*4]) {
__asm__ volatile (
// load all rows from M2, i.e. Ba, Bb, Bc an Bd
"movups 0x00(%1), %%xmm4 \n" // xmm4 := Ba
"movups 0x10(%1), %%xmm5 \n" // xmm5 := Bb
"movups 0x20(%1), %%xmm6 \n" // xmm6 := Bc
"movups 0x30(%1), %%xmm7 \n" // xmm7 := Bd
#define BLOCK(offset) \
/* load n-th row from matrix M1, for example Aa, and clone it */ \
"movups "#offset"(%0), %%xmm0 \n" /* xmm0 := Aa */ \
"movaps %%xmm0, %%xmm1 \n" \
"movaps %%xmm0, %%xmm2 \n" \
"movaps %%xmm0, %%xmm3 \n" \
/* populate each element */ \
"shufps $0b00000000, %%xmm0, %%xmm0 \n" /* xmm0 := [Aa0, Aa0, Aa0, Aa0] */ \
"shufps $0b01010101, %%xmm1, %%xmm1 \n" /* xmm1 := [Aa1, Aa1, Aa1, Aa1] */ \
"shufps $0b10101010, %%xmm2, %%xmm2 \n" /* xmm2 := [Aa2, Aa2, Aa2, Aa2] */ \
"shufps $0b11111111, %%xmm3, %%xmm3 \n" /* xmm3 := [Aa3, Aa3, Aa3, Aa3] */ \
/* and mul by all M2 rows */ \
"mulps %%xmm4, %%xmm0 \n" /* xmm0 := [Aa0*Ba0, Aa0*Ba1, Aa0*Ba2, Aa0*Ba3] */ \
"mulps %%xmm5, %%xmm1 \n" \
"mulps %%xmm6, %%xmm2 \n" \
"mulps %%xmm7, %%xmm3 \n" \
/* finally calculate n-th row of resultant matrix */ \
"addps %%xmm1, %%xmm0 \n" \
"addps %%xmm3, %%xmm2 \n" \
"addps %%xmm2, %%xmm0 \n" \
"movups %%xmm0, "#offset"(%2) \n"
BLOCK(0x00)
BLOCK(0x10)
BLOCK(0x20)
BLOCK(0x30)
:
: "r" (mat1), "r" (mat2), "r" (mat3)
);
}