SIMD Programming

[UfnCod3r]

:لبخند:
سلام
در این تاپیک قصد دارم تا چیزایی که خودم از SIMD, SSE می دونم رو بطور ساده توضیح بدم بقیه
هم راحت یاد بگیرن و بیفتن روی چرخ :چشمک::شیطان:

اول از همه بهتره اینا رو بخونید:گیج:
http://en.wikipedia.org/wiki/SIMD
http://en.wikipedia.org/wiki/Streaming_SIMD_Extensions
-
-
SIMD مخفف Single instruction, multiple data هست . یعنی یک دستور چند داده
کارش از اسمش معلومه
ینی ما یک دستورالعمل رو مشخص کنیم و اون رو روی چند تا داده اعمال کنیم .
این رو در نظر بگیرد .

float a[4] = {1,1,1,1};
float b[4] = {2,2,2,2};
float c[4];

c[0] = a[0] + b[0];
c[1] = a[1] + b[1];
c[2] = a[2] + b[2];
c[3] = a[3] + b[3];


این کدمون الان حالت برداری داره .
دستورات هر خط یکی هستن فقط داده ها تغییر می کنن.
تو این کد ما 4 تا عدد اعشاری رو با 4 تای دیگه جمع کردیم و ذخیره کردیم .
در مجموع میشه 4 تا عمل جمع و 12 تا عمل خواندن و نوشتن .
خب اگه میشد به این پردازشگر بفهمونیم که از a[0 تا a[3 رو با b[0 تا b[3 جمع کن خیلی خوب میشد .:متفکر:
بله میشه .

به مرور زمان که پردازشگر ها سریع تر شدن ، چند هسته ای شدن ،کش امد و ی سری چیزه دیگه
ی سری دستور عمل ها هم اضافه شدن . درواقع ی هردفعه الحاقی های جدید میاد
یکی از اون ها SSE هست .
با امدن SSE هشت تا ثبات 128 بیتی امد به نام XMM0...XMM7
تو این ثبات ها میشه 4 تا float رو جا داد و ضرب و جمع و دستورای دیگه رو اعمال کرد . چیزای دیگم میشه.
بعد هردفعه نسخه های جدید تر از SSE امد و چیز هایی دیگه مثل اعداد صحیح و
دستورالعمل های جدید اضافه شد .

برای اینکه بفهمید CPU ی شما چه الحاقی هایی رو پشتیبانی می کنه می تونید از برنامه CPUZ استفاده کنید
درحال حاضر میشه گفت اشغال ترین پی سی های امروزی SSE1,2,3 رو پشتیبانی می کنن . پس از این بابت نگران نباشید.

خب برای استفاده از این دستورات ما چند تا راه داریم .

-از اسمبل استفاده کنیم (خریت)
اسمبل نوشتن اصلا گزینه خوبی نیست . خوانایی کد کم میشه. کلی بهمون فشار میاد
برا کامپایلر ها و معماری های مختلف باید کد متفاوت بزنید و کلی دردسر دیگه

- به کامپایلر بگیم خودش کد رو تولید کنه (خریت)
کامپایلر ها اونقدر ها هم زرنگ نیستن .
هرچی کار کنن نمی تونن کد موازی تولید کنن . تو این زمینه نمیشه کاریشون کرد.

راه اخرش اینه که از ی سری توایع از پیش تعریف شده استفاده کنیم.
این توابع در واقع هر کدومشون یک دستور هستن (بجز بعضی ها) که خود کامپایلر بر اساس اونا کد اسمبل رو تولید می کنه .
این طوری نه عذاب اسمبل نوشتن رو می کشیم و هم اینکه برای کامپایلر ها و پلتفرم های مختلف لازم نیست کد رو دست بزنیم

این توابع تو این سرفایل ها قرار دارن. نمی دونم چرا اینطوری نام گزاری کردن:متفکر:


SSE xmmintrin.h
SSE2 emmintrin.h
SSE3 pmmintrin.h
Suppl. SSE3 tmmintrin.h
SSE4.1 smmintrin.h
SSE4.2 nmmintrin.h
AES, PCLMUL wmmintrin.h
AVX immintrin.h


خب حالا همون مثال بالا رو با SSE می نویسیم

float a[4] = {1,1,1,1};
float b[4] = {2,2,2,2};
float c[4];

__m128 xmm0 = _mm_loadu_ps(a);
__m128 xmm1 = _mm_loadu_ps(b);
__m128 xmm2 = _mm_add_ps(xmm0, xmm1);
_mm_storeu_ps(c, xmm2);


الان تو این کد دیگه کلا 4 تا دستوربیشتر انجام نمیشه . MOVUPS, MOVUPS, MOVUPS, ADDPS

در مورد نام تابع ها اگه دقت کنید قلقش میاد دستتون
اینایی که مربوط به SSE هستن همه با _mm_ شروع می شن.
و در اخر هم با ps_ یا ss_ تموم میشن . چیزایی دیگم هست که خودتون باید بفهمید
ps از
packed single گرفته شده
منظور از packed اینه که این عمل روی چهار تا عنصر این ثبات اعمال میشه
single هم منظورش single precision floating point . یا همون اعدادد ممیز شناور تک دقت هست. float خودمون.:لبخندساده::گیج:

ss از single single گرفته شده .
منظور از single اول اینه که این عمل فقط روی 32 بیته اول یا همون عنصر اول اعمال میشه و 3 تای دیگه رو کاری نداره . که درواقع کدمون از حالت موازی خارج میشه .
یکی از مشکلات کدی که کامپایلر تولید می کنه همینه . نمی تونه کد موازی تولید کنه و همش تک تک کار میکنه
single اخرهم به معنای همون single precision floating point هست که گفتم.

بین این ها هم عملکرد تابع رو مشخص میکنه .

پس خیل راحت میشه کار

__m128 _mm_add_ps(__m128 a , __m128 b );

رو فهمید.( ن پ ن بیام بگم add یعنی جمع کردن :عصبانی++:)

اول بهتره m128__ رو توضیح بدم . این قضیش مثل همین توابع هست .
این ی ساختار هست مثل همه ی ساختار های دیگه.
کامپایلر خودش این رو با ثبات های XMM0..XMM7 جایگزین می کنه .
درصورت کمبودن ثبات هم از استک استفاده می کنه .
بهتره این رو فقط داخل توابع استفاده کنید نه به عنوان عضو کلاس یا امثالش.
برای لود کردن داده ها داخل ثبات و برعکس هم باید از
load, loadu , store, storeu استفاده بشه .
u در اینجا به معنی UnAlignment هست
وفتی از _mm_load_ps و _mm_store_ps استفاده میشه Data-align ادرسی که می دیم حتما باید 16 باشه .
ینی ادرسش بایدضریبی از 16 باشه . این حالت خیلی سریع تره.
دلیلش هم برمیگرده به معماری سیستم و انتقال اطلاعات از رم به سی پی یو و برعکس .
این data-alignment هم دنیای هست واس خودش . :لبخند:
ولی با _mm_stoeu_ps, _mm_loadu_ps دیگه کاری به این کارا نداره ولی سرعتش کمتره .

خب حالا می خوایمکد رو تغییر بدیم تا سریع عمل کنه


#define MS_ALIGN16 __declspec(align(16))

MS_ALIGN16 float a[4] = {1,1,1,1};
MS_ALIGN16 float b[4] = {2,2,2,2};
MS_ALIGN16 float c[4];

__m128 xmm0 = _mm_load_ps(a);
__m128 xmm1 = _mm_load_ps(b);
__m128 xmm2 = _mm_add_ps(xmm0, xmm1);
_mm_store_ps(c, xmm2);


با صفتی که به متغیر ها دادیم (__declspec(align(16))) الان اونا در جایی قرار می گیرن که ادرسشون ضریبی از 16 باشه.
البته این برا MSVC کار میکنه برا کامپایلر های دیگه چیز دیگه ای باید بنویسید.
...
تامن بقیش رو می نویسم شما تشکر این ارسال رو بزنید:عصبانی::قلب:

[UfnCod3r]

توضیحات اکثر توابع درنت موجوده . فقط کافیه نام تابع رو تو گوگل بنویسد .

ولی ی سری توابع یکمی پیچیدگی دارن و من اونا رو می گم :لبخند:
این رو در نظر بگیرید

void Float4InvMul_Normal(const float* a, const float* b, float* out)
{
out[0] = sqrt(a[0] * b[3]);
out[1] = sqrt(a[1] * b[2]);
out[2] = sqrt(a[2] * b[1]);
out[3] = sqrt(a[3] * b[0]);
}


اگه بخوایم این رو با SSE بویسم باید یه جوری b رو برعکس کنیم با بقیه مشکلی نیست . این کار رو با تابع _mm_shuffle_ps می توینم انجام بدیم .
این تابع سه تا پرامتر داره . دو تا ثبات و ی عدد 1 بایتی که مشحص می کنه چی رو برداره و کجا بزاره .
این که رو چه حسابی این کارو می کنه برمی گرده به بیت های عدد.:گیج:

این کد رو در نظر بگیرد.

MS_ALIGN16 float a[4] = {10, 20, 30, 40};
_m128 xmm0 = _mm_load_ps(a);
_m128 xmm1 = _mm_shuffle_ps(xmm0, xmm0, 3);

//Result
//XMM0 == [40, 30, 20, 10]
//XMM1 == [10, 10, 10, 40]


من عدد 3 وارد کردم که باینریش میشه

11 00 00 00

از سمت راست دو بیت اول مشخص می کنه که در عنصر اول چی قرار بگیره . دو بیت دوم مشخص می کنه که در عنصر دوم چی قرار بگیره و الی اخر ..
دو تا بیت اول که در مبنای 10 میشه 3 یعنی اینکه عنصر سوم xmm0 رو بزار تو عنصر صفرمxmm1.
دو تا بیت دوم که صفر هستن ینی اینکه عنصر صفرم xmm0 رو بزار تو عنصر اول xmm1
دو تا بیت سوم که صفر هستن ینی اینکه عنصر صفرم xmm0 رو بزار عنصر تو دوم xmm1
دو تا بیت چهارم که صفر هستن ینی اینکه عنصر صفرم xmm0 رو بزار تو عنصر سوم xmm1

درواقع مقدار هر 2 بیت می گه کدوم عنصر رو بردار و جایگاه هر 2 بیت هم میگه کجا بزارش :گیج:

البته شما لازم نیست با بیت ها سر و کله بزنید می تونید از _MM_SHUFFLE استقاده کنید.
به این صورت :

_m128 xmm1 = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(0,0,0,3));


در ضمن اگه دقت کرده باشین این تابع دو تا ثبات میگیره و مقادیر سمت چپ رو از ثبات اولی بر میداره و نصف دیگه رو از دومی . (وارد حاشیه نمی شیم تا این حد بسه :لبخندساده:)

اینم ی نمونه از برعکس کردن

MS_ALIGN16 float a[4] = {10, 20, 30, 40};
_m128 xmm0 = _mm_load_ps(a);
_m128 xmm1 = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(0,1,2,3));
_mm_store_ps(a, xmm1);
//Result =>> a[0] == 40, a[1] == 30, a[2] == 20, a[3] == 10

[UfnCod3r]

خب حالا ی چند تا تابع رو با SSE پیاده سازی می کنیم تا کم کم راه بیفتیم .:شیطان:

تابع ساده Lerp یا همون Linear interpolatio (http://en.wikipedia.org/wiki/Linear_interpolation) رو برای بردار 4 بعدی می نویسم.


void Vec4_Lerp(const float* a, const float* b, float alpha, float* out)
{
out[0] = a[0] + (b[0] - a[0]) * alpha;
out[1] = a[1] + (b[1] - a[1]) * alpha;
out[2] = a[2] + (b[2] - a[2]) * alpha;
out[3] = a[3] + (b[3] - a[3]) * alpha;
}

َ
واس اینکه حرفه ای بازی در بیاریم اسمبل هاش رو هم می زارم :قهقهه:


?Vec4_Lerp@@YAXPBM0MPAM@Z PROC ; Vec4_Lerp, COMDAT

; 82 : {

push ebp
mov ebp, esp

; 83 : out[0] = a[0] + (b[0] - a[0]) * alpha;

mov eax, DWORD PTR _a$[ebp]
mov ecx, DWORD PTR _b$[ebp]
fld DWORD PTR [ecx]
mov edx, DWORD PTR _out$[ebp]
fsub DWORD PTR [eax]
fld DWORD PTR _alpha$[ebp]
fld ST(0)
fmulp ST(2), ST(0)
fld DWORD PTR [eax]
faddp ST(2), ST(0)
fxch ST(1)
fstp DWORD PTR [edx]

; 84 : out[1] = a[1] + (b[1] - a[1]) * alpha;

fld DWORD PTR [ecx+4]
fsub DWORD PTR [eax+4]
fmul ST(0), ST(1)
fadd DWORD PTR [eax+4]
fstp DWORD PTR [edx+4]

; 85 : out[2] = a[2] + (b[2] - a[2]) * alpha;

fld DWORD PTR [ecx+8]
fsub DWORD PTR [eax+8]
fmul ST(0), ST(1)
fadd DWORD PTR [eax+8]
fstp DWORD PTR [edx+8]

; 86 : out[3] = a[3] + (b[3] - a[3]) * alpha;

fld DWORD PTR [ecx+12]
fsub DWORD PTR [eax+12]
fmulp ST(1), ST(0)
fadd DWORD PTR [eax+12]
fstp DWORD PTR [edx+12]

; 87 : }

pop ebp
ret 0
?Vec4_Lerp@@YAXPBM0MPAM@Z ENDP ; Vec4_Lerp





void Vec4_Lerp_SSE(const float* a, const float* b, float alpha, float* out)
{
__m128 xmmA = _mm_load_ps(a);
__m128 xmmB = _mm_load_ps(b);
__m128 xmmAlpha = _mm_set_ps1(alpha);
_mm_store_ps(out, _mm_add_ps(xmmA, _mm_mul_ps(_mm_sub_ps(xmmB, xmmA), xmmAlpha)));
}


ASM

?Vec4_Lerp_SSE@@YAXPBM0MPAM@Z PROC ; Vec4_Lerp_SSE, COMDAT

; 90 : {

push ebp
mov ebp, esp
and esp, -8 ; fffffff8H

; 91 : __m128 xmmA = _mm_load_ps(a);

mov eax, DWORD PTR _a$[ebp]

; 92 : __m128 xmmB = _mm_load_ps(b);

mov edx, DWORD PTR _b$[ebp]
movaps xmm2, XMMWORD PTR [edx]

; 93 : __m128 xmmAlpha = _mm_set_ps1(alpha);
; 94 : _mm_store_ps(out, _mm_add_ps(xmmA, _mm_mul_ps(_mm_sub_ps(xmmB, xmmA), xmmAlpha)));

subps xmm2, XMMWORD PTR [eax]
movss xmm0, DWORD PTR _alpha$[ebp]
mov ecx, DWORD PTR _out$[ebp]
shufps xmm0, xmm0, 0
mulps xmm2, xmm0
addps xmm2, XMMWORD PTR [eax]
movaps XMMWORD PTR [ecx], xmm2

; 95 : }

mov esp, ebp
pop ebp
ret 0
?Vec4_Lerp_SSE@@YAXPBM0MPAM@Z ENDP ; Vec4_Lerp_SSE



این باید دو برابری از حالت معمولیش سریع تر باشه .:شیطان:

حالا تابع بدست اوردن فاصله دو بردار 4 بعدی رو بدست میاریم .
الگوریتمش بدین صورته که اول دو تا بردار رو از همدیگه کم می کنیم و بعد طول اون رو حساب می کنیم . برا محاسبه طول باید در هم ظربشون کنیم و بعد x, y, z, w رو با هم جمع کنیم .


float Vec4_Distance(float* a, float* b)
{
float xx = (a[0]-b[0]) * (a[0]-b[0]);
float yy = (a[1]-b[1]) * (a[1]-b[1]);
float zz = (a[2]-b[2]) * (a[2]-b[2]);
float ww = (a[3]-b[3]) * (a[3]-b[3]);
return sqrt(xx + yy + zz + ww);
}


ASM

?Vec4_Distance@@YAMPAM0@Z PROC ; Vec4_Distance, COMDAT

; 99 : {

push ebp
mov ebp, esp
mov eax, DWORD PTR _a$[ebp]
mov ecx, DWORD PTR _b$[ebp]
fld DWORD PTR [eax+4]
fsub DWORD PTR [ecx+4]
fld DWORD PTR [eax]
fsub DWORD PTR [ecx]
fld DWORD PTR [eax+8]
fsub DWORD PTR [ecx+8]
fld DWORD PTR [eax+12]
fsub DWORD PTR [ecx+12]

; 100 : float xx = (a[0]-b[0]) * (a[0]-b[0]);

fld ST(2)
fmulp ST(3), ST(0)

; 101 : float yy = (a[1]-b[1]) * (a[1]-b[1]);
; 102 : float zz = (a[2]-b[2]) * (a[2]-b[2]);
; 103 : float ww = (a[3]-b[3]) * (a[3]-b[3]);
; 104 : return sqrt(xx + yy + zz + ww);

fxch ST(2)
fstp DWORD PTR tv267[ebp]
fld DWORD PTR tv267[ebp]
fld ST(3)
fmulp ST(4), ST(0)
fxch ST(3)
fstp DWORD PTR tv263[ebp]
fld DWORD PTR tv263[ebp]
faddp ST(3), ST(0)
fmul ST(0), ST(0)
fstp DWORD PTR tv258[ebp]
fld DWORD PTR tv258[ebp]
faddp ST(2), ST(0)
fmul ST(0), ST(0)
fstp DWORD PTR tv254[ebp]
fadd DWORD PTR tv254[ebp]
fstp DWORD PTR tv251[ebp]
fld DWORD PTR tv251[ebp]
call __CIsqrt
fstp DWORD PTR tv244[ebp]
fld DWORD PTR tv244[ebp]

; 105 : }

pop ebp
ret 0







float Vec4_Distance_SSE(float* a, float* b)
{
__m128 xmm0 = _mm_sub_ps(_mm_load_ps(a), _mm_load_ps(b));
xmm0 = _mm_mul_ps(xmm0, xmm0);
__m128 yy = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(3,2,1,1));
__m128 zz = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(3,2,1,2));
__m128 ww = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(3,2,1,3));
return _mm_cvtss_f32(_mm_sqrt_ss(_mm_add_ss(xmm0, _mm_add_ss(yy, _mm_add_ss(zz, ww)))));
}



ASM

?Vec4_Distance_SSE@@YAMPAM0@Z PROC ; Vec4_Distance_SSE, COMDAT

; 108 : {

push ebp
mov ebp, esp
and esp, -8 ; fffffff8H
sub esp, 8

; 109 : __m128 xmm0 = _mm_sub_ps(_mm_load_ps(a), _mm_load_ps(b));

mov eax, DWORD PTR _a$[ebp]
mov ecx, DWORD PTR _b$[ebp]
movaps xmm0, XMMWORD PTR [eax]
subps xmm0, XMMWORD PTR [ecx]

; 110 : xmm0 = _mm_mul_ps(xmm0, xmm0);

mulps xmm0, xmm0

; 111 : __m128 yy = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(3,2,1,1));
; 112 : __m128 zz = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(3,2,1,2));
; 113 : __m128 ww = _mm_shuffle_ps(xmm0, xmm0, _MM_SHUFFLE(3,2,1,3));

movaps xmm1, xmm0
shufps xmm1, xmm0, 231 ; 000000e7H
movaps xmm2, xmm0
shufps xmm2, xmm0, 230 ; 000000e6H

; 114 : return _mm_cvtss_f32(_mm_sqrt_ss(_mm_add_ss(xmm0, _mm_add_ss(yy, _mm_add_ss(zz, ww)))));

addss xmm2, xmm1
movaps xmm1, xmm0
shufps xmm1, xmm0, 229 ; 000000e5H
addss xmm1, xmm2
addss xmm0, xmm1
sqrtss xmm0, xmm0
movss DWORD PTR tv180[esp+8], xmm0
fld DWORD PTR tv180[esp+8]

; 115 : }

mov esp, ebp
pop ebp
ret 0





این یکمی پیچیدگیش بیشتره .
a, b رو بارگزاری کردیم و از همدیگه کم کردیم و ضرب در خودش کردیم .
حالا باید چهار تا عددی که تو ثبات XMM هستن رو با هم جمع کنیم . برای این کار دستوری وجود نداره پس باید ی فکر دیگه کنیم . (البته در SSE4.1 ی چی هست که بعدا می گم :چشمک: )
برای اینکار ما می توینم y, z, w رو ببریم تو خونه ی اول ثبات و بعد با _mm_add_ss جمعشون کنیم .
درسته که این قسمتش دو تا suffle استفاده کردیم و و از حالت برداری هم خارج شدیم ولی خب داده های ما تو چند تادستور قبل قشنگ ضربو تقسیم شدن و درحال حاضر داخل ثبات هم هستن
به علاوه اینکه جذری که SSE داره خیلی سریعه و درکل تو این کد هم باز سرعت بیشتری داریم :شیطان:
تابع _mm_cvtss_f32 هم نام گزاریش این طوری الهام گرفته شده .

ctv == convert
ss == single single
f32 == float

ینی عنصر اول یه ثبات XMM رو به یک float تبدیل می کنه . که البته تبدیلی وجود نداره . چون float رو به float تبدیل نمی کنن . فقط یه MOVSS انجام میده .

[dearsahar]

خیلی عالی بود