لقد مرت سنوات عديدة منذ أن تعرفت على التعليمات MMX و SSE و AVX اللاحقة على معالجات Intel. في وقت من الأوقات ، بدوا وكأنهم نوع من السحر على خلفية مجمّع x86 ، الذي كان شيئًا دنيويًا منذ فترة طويلة. لقد ربطوني كثيرًا قبل بضع سنوات كان لدي فكرة أن أكتب برنامج عرض البرامج الخاص بي للعبة واحدة مشهورة. الشيء الذي وعد بهذه التعليمات وعدني بذلك. في مرحلة ما ، فكرت حتى في كتابتها. لكن اتضح أن كتابة النص أكثر تعقيدًا من التعليمات البرمجية.
في ذلك الوقت ، كنت أرغب في تجنب المشاكل مع دعم المعالجات المختلفة. أردت أن أتمكن من التحقق من العارض الخاص بي على أقصى مبلغ متاح. لا يزال لدي أصدقاء مع معالجات AMD القديمة ، وكان سقفهم SSE3. لذلك ، في ذلك الوقت قررت أن أقتصر على SSE3 كحد أقصى. لذلك كانت هناك مكتبة رياضية متجهة ، أقل قليلاً من تنفيذها بالكامل على SSE ، مع إدراج نادر قبل SSE3. ومع ذلك ، في مرحلة ما ، تساءلت عن الحد الأقصى للأداء الذي يمكنني الحصول عليه من المعالج لعدد من العمليات الحسابية الحسابية. إحدى هذه العمليات هي مضاعفة 4 مصفوفات في 4.
في الواقع ، قررت أن أقوم بهذا العمل أكثر من أجل الترفيه. لقد كتبت بالفعل وكنت أستخدم ضرب المصفوفات لعرض برنامجي على SSE ويبدو أنه كافٍ بالنسبة لي. ولكن بعد ذلك قررت أن أرى عدد التدابير التي يمكنني الضغط عليها من حيث المبدأ من خلال ضرب مصفوفتين عشريتين 4 × 4. في SSE الحالي ، هذه 16 دورة ساعة. صحيح أن التحول الأخير إلى
IACA 3 بدأ في إظهار 19 ، حيث بدأت في كتابة 1 * لبعض التعليمات بدلاً من 0 *. يبدو أنه في وقت سابق كان مجرد خلل في المحلل.
لفترة وجيزة حول المرافق المستخدمة
لتحليل التعليمات البرمجية ، استخدمت أداة
Intel Code Code Analyzer الشهيرة. للتحليل ، استخدم بنية Haswell (HSW) ، كحد أدنى مع دعم AVX2. كتابة التعليمات البرمجية هي أيضًا مريحة جدًا للاستخدام:
دليل إنتل إنترسكس ودليل تحسين إنتل .
للتجميع ، أستخدم مجتمع MSVS 2017 من وحدة التحكم. أكتب الرمز في الإصدار مع الجوهر. تكتب مرة واحدة ، وعادة ما تعمل على الفور على منصات مختلفة. بالإضافة إلى ذلك ، لا يدعم برنامج التحويل البرمجي x64 VC ++ أداة التجميع المضمنة ، ولكني أريد أن يعمل تحت x64 أيضًا.
نظرًا لأن هذه المقالة تتجاوز إلى حد ما مستوى المبتدئين في برمجة SIMD ، فلن أصف التسجيلات أو التعليمات أو رسم (أو حلاقة) صور جميلة ومحاولة تعلم البرمجة باستخدام تعليمات SIMD. موقع إنتل مليء بالوثائق الممتازة والواضحة والمفصلة.
أردت أن أجعل كل شيء أسهل ... ولكن اتضح كما هو الحال دائمًا
هذا هو المكان الذي تبدأ فيه اللحظة ، مما يعقد الكثير من التنفيذ والمقال. لذلك ، سوف أتطرق إلى ذلك قليلاً. ليس من المثير للاهتمام بالنسبة لي أن أكتب ضرب المصفوفة مع تخطيط صف قياسي من العناصر. من احتاج لذلك ، ودرسوا في الجامعات أو بمفردهم. هدفنا هو الإنتاجية. أولاً ، لقد تحولت إلى تخطيط العمود منذ فترة طويلة. يستند عارض البرامج الخاص بي إلى OpenGL API ، وبالتالي ، لتجنب عمليات النقل غير الضرورية ، بدأت في تخزين العناصر في أعمدة. هذا مهم أيضًا لأن ضرب المصفوفات ليس بالغ الأهمية. ضرب المصفوفات 2-5-10 جيدا. هذا كل ما في الأمر. ثم نضرب المصفوفة النهائية بآلاف أو ملايين القمم. وهذه العملية أكثر أهمية. يمكنك ، بالطبع ، تبديل كل مرة. ولكن لماذا ، إذا كان من الممكن تجنب ذلك.
ولكن نعود إلى المصفوفات حصريا. حددنا التخزين في الأعمدة. ومع ذلك ، يمكن أن يكون أكثر تعقيدًا. من الأنسب بالنسبة لي تخزين العناصر الأقدم للمتجهات وصفوف المصفوفة في سجلات SIMD بحيث يكون
x في أعلى تعويم (الفهرس 3) ، و
w في القاصر (الفهرس 0). هنا ، على ما يبدو ، سيتعين علينا التراجع مرة أخرى عن السبب.
الشيء هو أنه في عارض البرامج في متجه ، يجب عليك معالجة المكون
w في كثير من الأحيان (يتم تخزين
1 / z هناك) ، وهو ملائم جدًا للقيام بذلك من خلال إصدار
_ss من العملية (العمليات حصريًا مع المكون في العوامة السفلية
لسجل xmm ) ، دون لمس
x ، y ، ض . لذلك ، في سجل SSE ، يتم تخزين المتجه في ترتيب مفهومة
x ، y ، z ، w ، وذاكرة في الاتجاه المعاكس
w ، z ، y ، x .
علاوة على ذلك ، يتم تنفيذ جميع خيارات الضرب أيضًا من خلال الوظائف الفردية. يتم ذلك لأنني أستخدمها لاستبدال الخيار المطلوب اعتمادًا على نوع التعليمات المدعومة.
وصف جيد هنا.نقوم بتطبيق الوظائف الأساسية
الضرب مع الحلقات ، ترتيب الصف
خيار تخطيط خط العناصرfor (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { r[i][j] = 0.f; for (int k = 0; k < 4; ++k) { r[i][j] += m[i][k] * n[k][j]; } } }
كل شيء بسيط وواضح هنا. لكل عنصر نقوم به 4 ضربات و 3 إضافات. في المجموع ، هذه 64 مضاعفات و 48 إضافات. وهذا دون الأخذ بعين الاعتبار قراءة عناصر السجل.
كل شيء حزين ، باختصار. بالنسبة لهذا الخيار ، بالنسبة للدورة الداخلية ، أصدرت IACA:
3.65 دورات ساعة لتجميع x86 و 2.97 ساعات لتجميع x64 . لا تسأل لماذا الأرقام الكسرية. لا اعرف. IACA 2.1 لم يعاني من هذا. على أي حال ، يجب ضرب هذه الأرقام بحوالي 4 * 4 * 4 = 64. حتى إذا كنت تأخذ x64 ، تكون النتيجة حوالي 192 مقياسًا. من الواضح أن هذا تقدير تقريبي. لا أرى فائدة من تقييم الأداء بشكل أكثر دقة لهذا الخيار.
تنفيذ الحلقات ، أمر العمود
المصفوفة المنقولة ، إعادة ترتيب مؤشرات الصف والعمود for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { r[j][i] = 0.f; for (int k = 0; k < 4; ++k) { r[j][i] += m[k][i] * n[j][k]; } } }
دورة الضرب ، تخزين موجه SIMD
يضاف تخزين الخطوط بالترتيب العكسي في الذاكرة for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { r[j][3-i] = 0.f; for (int k = 0; k < 4; ++k) { r[j][3-i] += m[k][3-i] * n[j][3-k]; } } }
يبسط هذا التنفيذ إلى حد ما فهم ما يحدث في الداخل ، ولكن من الواضح أنه غير كافٍ.
دروس المساعد
من أجل راحة فهم وكتابة المرجع وكود التصحيح ، من الملائم تنفيذ فئتين إضافيتين. لا شيء أكثر من ذلك ، كل شيء هو فقط من أجل الفهم. ألاحظ أن تنفيذ فئات المتجهات والمصفوفة الكاملة هو سؤال صعب منفصل ، ولا يتم تضمينه في موضوع هذه المقالة.
فئات المتجهات والمصفوفة struct alignas(sizeof(__m128)) vec4 { union { struct { float w, z, y, x; }; __m128 fmm; float arr[4]; }; vec4() {} vec4(float a, float b, float c, float d) : w(d), z(c), y(b), x(a) {} static bool equ(float const a, float const b, float t = .00001f) { return fabs(ab) < t; } bool operator == (vec4 const& v) const { return equ(x, vx) && equ(y, vy) && equ(z, vz) && equ(w, vw); } }; struct alignas(sizeof(__m256)) mtx4 {
الوظيفة المرجعية للاختبارات
نظرًا لأن الترتيب المقبول للعناصر في المصفوفة يعقد الفهم كثيرًا ، فلن نتضايق أيضًا من الوظيفة المرجعية القابلة
للفهم ، والتي ستظهر في عمليات التنفيذ المستقبلية أن كل شيء يعمل بشكل صحيح. سنقوم بمقارنة النتائج اللاحقة معها.
لإنشائه ، ما عليك سوى أخذ الدورة وتوسيعها void mul_mtx4_mtx4_unroll(__m128* const _r, __m128 const* const _m, __m128 const* const _n) { mtx4 const& m = **reinterpret_cast<mtx4 const* const*>(&_m); mtx4 const& n = **reinterpret_cast<mtx4 const* const*>(&_n); mtx4& r = **reinterpret_cast<mtx4* const*>(&_r); r._00 = m._00*n._00 + m._01*n._10 + m._02*n._20 + m._03*n._30; r._01 = m._00*n._01 + m._01*n._11 + m._02*n._21 + m._03*n._31; r._02 = m._00*n._02 + m._01*n._12 + m._02*n._22 + m._03*n._32; r._03 = m._00*n._03 + m._01*n._13 + m._02*n._23 + m._03*n._33; r._10 = m._10*n._00 + m._11*n._10 + m._12*n._20 + m._13*n._30; r._11 = m._10*n._01 + m._11*n._11 + m._12*n._21 + m._13*n._31; r._12 = m._10*n._02 + m._11*n._12 + m._12*n._22 + m._13*n._32; r._13 = m._10*n._03 + m._11*n._13 + m._12*n._23 + m._13*n._33; r._20 = m._20*n._00 + m._21*n._10 + m._22*n._20 + m._23*n._30; r._21 = m._20*n._01 + m._21*n._11 + m._22*n._21 + m._23*n._31; r._22 = m._20*n._02 + m._21*n._12 + m._22*n._22 + m._23*n._32; r._23 = m._20*n._03 + m._21*n._13 + m._22*n._23 + m._23*n._33; r._30 = m._30*n._00 + m._31*n._10 + m._32*n._20 + m._33*n._30; r._31 = m._30*n._01 + m._31*n._11 + m._32*n._21 + m._33*n._31; r._32 = m._30*n._02 + m._31*n._12 + m._32*n._22 + m._33*n._32; r._33 = m._30*n._03 + m._31*n._13 + m._32*n._23 + m._33*n._33; }
تم رسم الخوارزمية الكلاسيكية بوضوح هنا ، من الصعب ارتكاب خطأ (ولكن يمكنك :-)). على ذلك ، أصدرت IACA:
تدابير x86 - 69.95 ، x64 - 64 تدابير . هنا حوالي 64 دورة وسننظر في تسريع هذه العملية في المستقبل.
تنفيذ SSE
خوارزمية SSE الكلاسيكية
لماذا الكلاسيكية؟ لأنه منذ فترة طويلة في تنفيذ
FVec كجزء من MSVS. بادئ ذي بدء ، سنكتب كيف نقدم عناصر المصفوفة في سجلات SSE. يبدو بالفعل أبسط هنا. مجرد مصفوفة منقولة.
// 00, 10, 20, 30 // m[0] - SIMD / 01, 11, 21, 31 // m[1] 02, 12, 22, 32 // m[2] 03, 13, 23, 33 // m[3]
نأخذ رمز
إلغاء التسجيل من المتغير أعلاه. بطريقة ما هو غير ودي لـ SSE. تتكون أول مجموعة من الصفوف من نتائج عمود المصفوفة الناتجة:
r._00 ، r._01 ، r._02 ، r._03 . لدينا هذا العمود ، لكننا بحاجة إلى صف. نعم ، و
m ،
n تبدو غير ملائمة للحسابات. لذلك ، نعيد ترتيب خطوط الخوارزمية بحيث تكون النتيجة
r حكيمة.
// , r[0] r00 = m00*n00 + m01*n10 + m02*n20 + m03*n30; r10 = m10*n00 + m11*n10 + m12*n20 + m13*n30; r20 = m20*n00 + m21*n10 + m22*n20 + m23*n30; r30 = m30*n00 + m31*n10 + m32*n20 + m33*n30; // , r[1] r01 = m00*n01 + m01*n11 + m02*n21 + m03*n31; r11 = m10*n01 + m11*n11 + m12*n21 + m13*n31; r21 = m20*n01 + m21*n11 + m22*n21 + m23*n31; r31 = m30*n01 + m31*n11 + m32*n21 + m33*n31; // , r[2] r02 = m00*n02 + m01*n12 + m02*n22 + m03*n32; r12 = m10*n02 + m11*n12 + m12*n22 + m13*n32; r22 = m20*n02 + m21*n12 + m22*n22 + m23*n32; r32 = m30*n02 + m31*n12 + m32*n22 + m33*n32; // , r[3] r03 = m00*n03 + m01*n13 + m02*n23 + m03*n33; r13 = m10*n03 + m11*n13 + m12*n23 + m13*n33; r23 = m20*n03 + m21*n13 + m22*n23 + m23*n33; r33 = m30*n03 + m31*n13 + m32*n23 + m33*n33;
لكن هذا بالفعل أفضل بكثير. ما الذي نراه في الواقع؟ وفقًا لأعمدة الخوارزمية في كل مجموعة ، لدينا صفوف المصفوفة
m المعنية:
م [0] = {00،10،20،30} ، م [1] = {01،11،21،31} ، م [2] = {02،12،22،32} ، م [3] = {03،13،23،33} ،
التي يتم ضربها في نفس عنصر المصفوفة
n . على سبيل المثال ، بالنسبة للمجموعة الأولى فهي:
n._00، n._10، n._20، n._30 . وعناصر المصفوفة
n لكل مجموعة من صفوف الخوارزمية تكمن مرة أخرى في صف واحد من المصفوفة.
ثم كل شيء بسيط: نأخذ ببساطة صفوف المصفوفة
m حسب الفهرس ، أما بالنسبة للعناصر
n ، فنحن نأخذ صفها ونقوم
بتبديلها عبر جميع عناصر السجل الأربعة من خلال تعليمات
التبديل ، للضرب في صف المصفوفة
m في السجل. على سبيل المثال ، بالنسبة للعنصر
n._00 (تذكر أن تحوله في السجل له فهرس 3) ، فسيكون:
_mm_shuffle_ps (n [0]، n [0]، _MM_SHUFFLE (3،3،3،3))
في شكل مبسط ، تبدو الخوارزمية كما يلي:
// n[0]={00,10,20,30} r[0] = m[0] * n00 + m[1] * n10 + m[2] * n20 + m[3] * n30; // n[1]={01,11,21,31} r[1] = m[0] * n01 + m[1] * n11 + m[2] * n21 + m[3] * n31; // n[2]={02,12,22,32} r[2] = m[0] * n02 + m[1] * n12 + m[2] * n22 + m[3] * n32; // n[3]={03,13,23,33} r[3] = m[0] * n03 + m[1] * n13 + m[2] * n23 + m[3] * n33;
تطبيق SSE الأساسي void mul_mtx4_mtx4_sse_v1(__m128* const r, __m128 const* const m, __m128 const* const n) { r[0] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[0], n[0], _MM_SHUFFLE(0,0,0,0))))); r[1] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[1], n[1], _MM_SHUFFLE(0,0,0,0))))); r[2] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[2], n[2], _MM_SHUFFLE(0,0,0,0))))); r[3] = _mm_add_ps( _mm_add_ps( _mm_mul_ps(m[0], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(3,3,3,3))), _mm_mul_ps(m[1], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(2,2,2,2)))), _mm_add_ps( _mm_mul_ps(m[2], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(1,1,1,1))), _mm_mul_ps(m[3], _mm_shuffle_ps(n[3], n[3], _MM_SHUFFLE(0,0,0,0))))); }
الآن نغير العناصر
n في الخوارزمية إلى
المراوغة المقابلة ، الضرب في
_mm_mul_ps ، المجموع بمقدار
_mm_add_ps ،
وبذلك تكون قد انتهيت. يعمل. لكن الكود يبدو أسوأ بكثير من الخوارزمية نفسها. إلى هذا الرمز ، أصدرت IACA:
x86 - 18.89 ، x64 - 16 دورة . هذا هو 4 مرات أسرع من سابقتها. في SSE تسجيل تعويم 4th. علاقة خطية تقريبا.
تزيين تنفيذ SSE
لا يزال ، في القانون يبدو فظيعا. سنحاول تحسين هذا عن طريق كتابة القليل من السكر النحوي.
يمكن للمترجم أن يدمج هذه الوظائف بشكل مثالي (على الرغم من عدم وجود __forceinline بأي شكل من الأشكال).
لذا فإن الكود يتحول ... void mul_mtx4_mtx4_sse_v2(__m128* const r, __m128 const* const m, __m128 const* const n) { r[0] = m[0]*shuf<3>(n[0]) + m[1]*shuf<2>(n[0]) + m[2]*shuf<1>(n[0]) + m[3]*shuf<0>(n[0]); r[1] = m[0]*shuf<3>(n[1]) + m[1]*shuf<2>(n[1]) + m[2]*shuf<1>(n[1]) + m[3]*shuf<0>(n[1]); r[2] = m[0]*shuf<3>(n[2]) + m[1]*shuf<2>(n[2]) + m[2]*shuf<1>(n[2]) + m[3]*shuf<0>(n[2]); r[3] = m[0]*shuf<3>(n[3]) + m[1]*shuf<2>(n[3]) + m[2]*shuf<1>(n[3]) + m[3]*shuf<0>(n[3]); }
وهي بالفعل أفضل بكثير وأكثر قابلية للقراءة. إلى ذلك ، أنتجت IACA النتيجة المتوقعة تقريبًا:
x86 - 19 (ولماذا لا تكون كسرية؟) ، X64 - 16 . في الواقع ، لم يتغير الأداء ، لكن الكود أكثر جمالا ومفهوما.
مساهمة ضئيلة في التحسين في المستقبل
دعونا نقدم تحسينًا آخر على مستوى الوظيفة التي ظهرت مؤخرًا في الإصدار الحديدي. عملية
الإضافة المتعددة (fma) .
fma (a، b، c) = a * b + c .
تنفيذ إضافة متعددة __m128 mad(__m128 const a, __m128 const b, __m128 const c) { return _mm_add_ps(_mm_mul_ps(a, b), c); }
لماذا هذا ضروري؟ بادئ ذي بدء ، للتحسين في المستقبل. على سبيل المثال ، يمكنك ببساطة استبدال
mad fma في الرمز النهائي من خلال نفس وحدات الماكرو التي تريدها. لكننا سنضع الأساس للتحسين الآن:
البديل مع إضافة متعددة void mul_mtx4_mtx4_sse_v3(__m128* const r, __m128 const* const m, __m128 const* const n) { r[0] = mad(m[0], shuf<3>(n[0]), m[1]*shuf<2>(n[0])) + mad(m[2], shuf<1>(n[0]), m[3]*shuf<0>(n[0])); r[1] = mad(m[0], shuf<3>(n[1]), m[1]*shuf<2>(n[1]) + mad(m[2], shuf<1>(n[1]), m[3]*shuf<0>(n[1])); r[2] = mad(m[0], shuf<3>(n[2]), m[1]*shuf<2>(n[2])) + mad(m[2], shuf<1>(n[2]), m[3]*shuf<0>(n[2])); r[3] = mad(m[0], shuf<3>(n[3]), m[1]*shuf<2>(n[3])) + mad(m[2], shuf<1>(n[3]), m[3]*shuf<0>(n[3])); }
IACA:
x86 - 18.89 ، x64 - 16 . كسري مرة أخرى. مع ذلك ، ينتج IACA أحيانًا نتائج غريبة. لم يتغير الرمز كثيرا. ربما أسوأ قليلا. لكن التحسين يتطلب أحيانًا مثل هذه التضحيات.
نحن نمر إلى الحفظ من خلال _mm_stream
توصي أدلة التحسين المختلفة مرة أخرى بعدم سحب ذاكرة التخزين المؤقت لعمليات الحفظ الجماعي. عادة ما يكون هذا مبررًا عندما تقوم بمعالجة القمم التي تبلغ آلاف أو أكثر. ولكن بالنسبة للمصفوفات ، ربما لا يكون هذا مهمًا. ومع ذلك ، سأضيفه على أي حال.
خيار حفظ الدفق void mul_mtx4_mtx4_sse_v4(__m128* const r, __m128 const* const m, __m128 const* const n) { _mm_stream_ps(&r[0].m128_f32[0], mad(m[0], shuf<3>(n[0]), m[1]*shuf<2>(n[0])) + mad(m[2], shuf<1>(n[0]), m[3]*shuf<0>(n[0]))); _mm_stream_ps(&r[1].m128_f32[0], mad(m[0], shuf<3>(n[1]), m[1]*shuf<2>(n[1])) + mad(m[2], shuf<1>(n[1]), m[3]*shuf<0>(n[1]))); _mm_stream_ps(&r[2].m128_f32[0], mad(m[0], shuf<3>(n[2]), m[1]*shuf<2>(n[2])) + mad(m[2], shuf<1>(n[2]), m[3]*shuf<0>(n[2]))); _mm_stream_ps(&r[3].m128_f32[0], mad(m[0], shuf<3>(n[3]), m[1]*shuf<2>(n[3])) + mad(m[2], shuf<1>(n[3]), m[3]*shuf<0>(n[3]))); }
لم يتغير شيء في الوقت المناسب هنا ، من الكلمة على الإطلاق. ولكن ، وفقًا للتوصيات ، لا نلمس الآن ذاكرة التخزين المؤقت مرة أخرى.
تنفيذ AVX
خيار AVX الأساسي
بعد ذلك ، ننتقل إلى المرحلة التالية من التحسين. يتم تضمين العوامة الرابعة في سجل SSE ، وفي AVX يبلغ بالفعل 8. أي ، هناك فرصة نظرية لتقليل عدد العمليات التي يتم إجراؤها ، وزيادة الإنتاجية إن لم يكن بمقدار النصف ، ثم 1.5 مرة على الأقل. لكن شيئًا يخبرني أنه لن يكون كل شيء بهذه البساطة مع الانتقال إلى AVX. هل يمكننا الحصول على البيانات اللازمة من السجلات المزدوجة؟
دعونا نحاول معرفة ذلك. مرة أخرى ، نكتب خوارزمية الضرب المستخدمة أعلاه. لا يمكنك القيام بذلك ، ولكن من الأسهل التعامل مع الرمز عندما يكون كل شيء قريبًا ولا يتعين عليك التمرير لأعلى نصف الصفحة.
// : 00, 10, 20, 30, 01, 11, 21, 31, 02, 12, 22, 32, 03, 13, 23, 33 // SSE: r0 = m0*n00 + m1*n10 + m2*n20 + m3*n30 r1 = m0*n01 + m1*n11 + m2*n21 + m3*n31 r2 = m0*n02 + m1*n12 + m2*n22 + m3*n32 r3 = m0*n03 + m1*n13 + m2*n23 + m3*n33
عند الإخراج ، نتوقع الحصول على النتيجة في
ymm = {r0: r1} و
ymm = {r2: r3} . إذا كان إصدار خوارزمياتنا في إصدار SSE معممًا للأعمدة ، فسنحتاج الآن إلى تعميمها على الصفوف. لذا فإن التصرف كما في حالة خيار SSE لن يعمل.
إذا أخذنا المصفوفة
m في التسجيلات
ymm ، نحصل على
ymm = {m0: m1} و
ymm = {m2: m3} على التوالي. في السابق ، كان لدينا أعمدة مصفوفة فقط في السجل ، والآن الأعمدة والصفوف.
إذا حاولت التصرف كما كان من قبل ، فيجب عليك ضرب
ymm = {m0: m1} في التسجيل
ymm = {n00، n00، n00، n00}: {n10، n10، n10، n10} . نظرًا لأن
n00 و
n01 موجودان في نفس الصف من المصفوفة
n ، بناءً على المجموعة المتاحة من تعليمات AVX ، فإن
تناثرها بواسطة
ymm سيكون مكلفًا. يعمل كل من
الخلط والتبديل بشكل منفصل لكل من أربع
عوامات (ارتفاع ومنخفض
xmm ) داخل سجلات
ymm .
إذا أخذنا
ymm من المصفوفة
n ، نحصل على كلا العنصرين
n00 و
n10 بأعلى 2
xmm داخل سجل
ymm .
{n00، n10، n20، n30}: {n01، n11، n21، n31} . عادة ، يكون مؤشر التعليمات الموجودة من 0 إلى 3. ويتناول العناوين العائمة داخل تسجيل
xmm واحد فقط من اثنين من تسجيل
ymm الداخلي. لا يمكن نقل
n10 من
xmm الأكبر سناً إلى الأصغر سناً. ثم يجب تكرار هذا التركيز عدة مرات. لا يمكننا تحمل هذه الخسارة في التدابير. من الضروري الخروج بشيء آخر.
اعتدنا على تعميم الأعمدة ، ولكن الآن الصفوف. لذلك ، سنحاول الذهاب بطريقة مختلفة قليلاً. نحتاج إلى الحصول على النتيجة في
{r0: r1} . هذا يعني أنه يجب تحسين الخوارزمية ليس في خطوط منفصلة من الخوارزمية ، ولكن في خطين في وقت واحد. وهنا ، ما كان ناقصًا في
خلط ورق اللعب سيكون إضافة لنا. ننظر إلى ما سيكون لدينا في تسجيلات
ymm عندما نعتبر المصفوفة
n .
n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31} n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33}
نعم ، نلاحظ أنه في أجزاء
xmm مختلفة من تسجيل
ymm لدينا عناصر
00 و
01 . يمكن مضاعفة الحالة بالتسجيل من خلال الأمر
permute في
{_00 ، _00 ، _00 ، _00}: {_ 01 ، _01 ، _01 ، _01} ، مع الإشارة إلى فهرس واحد فقط لكل من أجزاء
xmm . هذا هو بالضبط ما نحتاجه. في الواقع ، يتم استخدام المعاملات أيضًا في خطوط مختلفة. الآن فقط في تسجيل
ymm المقابل للضرب سيكون من الضروري الاحتفاظ بـ
{m0: m0} ، أي الصف الأول المكرر من المصفوفة
m .
لذا ، نرسم الخوارزمية بمزيد من التفاصيل. نقرأ الصفوف المزدوجة للمصفوفة
m في تسجيلات
ymm :
mm[0] = {m0:m0} mm[1] = {m1:m1} mm[2] = {m2:m2} mm[3] = {m3:m3}
ثم سنحسب الضرب على النحو التالي:
r0r1 = mm[0] * {n00,n00,n00,n00:n01,n01,n01,n01} + // permute<3,3,3,3>(n0n1) mm[1] * {n10,n10,n10,n10:n11,n11,n11,n11} + // permute<2,2,2,2>(n0n1) mm[2] * {n20,n20,n20,n20:n21,n21,n21,n21} + // permute<1,1,1,1>(n0n1) mm[3] * {n30,n30,n30,n30:n31,n31,n31,n31} // permute<0,0,0,0>(n0n1) r2r3 = mm[0] * {n02,n02,n02,n02:n03,n03,n03,n03} + // permute<3,3,3,3>(n2n3) mm[1] * {n12,n12,n12,n12:n13,n13,n13,n13} + // permute<2,2,2,2>(n2n3) mm[2] * {n22,n22,n22,n22:n23,n23,n23,n23} + // permute<1,1,1,1>(n2n3) mm[3] * {n32,n32,n32,n32:n33,n33,n33,n33} // permute<0,0,0,0>(n2n3)
نعيد كتابة أكثر وضوحا:
r0r1 = mm[0]*n0n1<3,3,3,3>+mm[1]*n0n1<2,2,2,2>+mm[2]*n0n1<1,1,1,1>+mm[3]*n0n1<0,0,0,0> r2r3 = mm[0]*n2n3<3,3,3,3>+mm[1]*n2n3<2,2,2,2>+mm[2]*n2n3<1,1,1,1>+mm[3]*n2n3<0,0,0,0>
أو بشكل مبسط:
r0r1 = mm[0]*n0n1<3> + mm[1]*n0n1<2> + mm[2]*n0n1<1> + mm[3]*n0n1<0> r2r3 = mm[0]*n2n3<3> + mm[1]*n2n3<2> + mm[2]*n2n3<1> + mm[3]*n2n3<0>
يبدو أن كل شيء واضح.
يبقى فقط لكتابة التنفيذ void mul_mtx4_mtx4_avx_v1(__m128* const r, __m128 const* const m, __m128 const* const n) { __m256 mm0 = _mm256_set_m128(m[0], m[0]); __m256 mm1 = _mm256_set_m128(m[1], m[1]); __m256 mm2 = _mm256_set_m128(m[2], m[2]); __m256 mm3 = _mm256_set_m128(m[3], m[3]); __m256 n0n1 = _mm256_load_ps(&n[0].m128_f32[0]); __m256 y1 = _mm256_permute_ps(n0n1, 0xFF);
فيما يلي الأرقام المثيرة للاهتمام من IACA:
x86 - 12.53 ، x64 - 12 . على الرغم ، بالطبع ، أردت أفضل. غاب عن شيء.
AVX الأمثل بالإضافة إلى السكر النحوي
يبدو أنه في الكود أعلاه ، لم يتم استخدام AVX بكامل طاقته. نجد أنه بدلاً من تعيين سطرين متطابقين في تسجيل
ymm ، يمكننا استخدام
البث ، والذي يمكن أن يملأ تسجيل
ymm بقيم
xmm متطابقة. أيضا ، على طول الطريق ، إضافة بعض "السكر النحوي" لوظائف AVX.
تنفيذ AVX المحسن __m256 operator + (__m256 const a, __m256 const b) { return _mm256_add_ps(a, b); } __m256 operator - (__m256 const a, __m256 const b) { return _mm256_sub_ps(a, b); } __m256 operator * (__m256 const a, __m256 const b) { return _mm256_mul_ps(a, b); } __m256 operator / (__m256 const a, __m256 const b) { return _mm256_div_ps(a, b); } template <int i> __m256 perm(__m256 const v) { return _mm256_permute_ps(v, _MM_SHUFFLE(i, i, i, i)); } template <int a, int b, int c, int d> __m256 perm(__m256 const v) { return _mm256_permute_ps(v, _MM_SHUFFLE(a, b, c, d)); } template <int i, int j> __m256 perm(__m256 const v) { return _mm256_permutevar_ps(v, _mm256_set_epi32(i, i, i, i, j, j, j, j)); } template <int a, int b, int c, int d, int e, int f, int g, int h> __m256 perm(__m256 const v) { return _mm256_permutevar_ps(v, _mm256_set_epi32(a, b, c, d, e, f, g, h)); } __m256 mad(__m256 const a, __m256 const b, __m256 const c) { return _mm256_add_ps(_mm256_mul_ps(a, b), c); } void mul_mtx4_mtx4_avx_v2(__m128* const r, __m128 const* const m, __m128 const* const n) { __m256 const mm[] { _mm256_broadcast_ps(m+0), _mm256_broadcast_ps(m+1), _mm256_broadcast_ps(m+2), _mm256_broadcast_ps(m+3) }; __m256 const n0n1 = _mm256_load_ps(&n[0].m128_f32[0]); _mm256_stream_ps(&r[0].m128_f32[0], mad(perm<3>(n0n1), mm[0], perm<2>(n0n1)*mm[1])+ mad(perm<1>(n0n1), mm[2], perm<0>(n0n1)*mm[3])); __m256 const n2n3 = _mm256_load_ps(&n[2].m128_f32[0]); _mm256_stream_ps(&r[2].m128_f32[0], mad(perm<3>(n2n3), mm[0], perm<2>(n2n3)*mm[1])+ mad(perm<1>(n2n3), mm[2], perm<0>(n2n3)*mm[3])); }
وهنا تكون النتائج أكثر إثارة للاهتمام. تنتج IACA أرقامًا:
x86 - 10 ، x64 - 8.58 ، والتي تبدو أفضل بكثير ، ولكن لا تزال غير مرتين.
خيار AVX + FMA (نهائي)
لنقم بمحاولة أخرى. الآن سيكون من المنطقي استدعاء مجموعة تعليمات FMA مرة أخرى ، حيث تمت إضافتها إلى المعالجات بعد AVX. فقط قم بتغيير
مول فردي
+ أضف لعملية واحدة. على الرغم من أننا ما زلنا نستخدم تعليمات الضرب لإعطاء المترجم المزيد من الفرص للتحسين ، والمعالج للتنفيذ المتوازي للمضاعفات. عادةً ما ألقي نظرة على الرمز الذي تم إنشاؤه في أداة التجميع للتأكد من الخيار الأفضل.
في هذه الحالة ، نحتاج إلى حساب
a * b + c * d + e * f + g * h . يمكنك القيام بهذا الجبين:
fma (a، b، fma (c، d، fma (e، f، g * h))) . ولكن ، كما نرى ، من المستحيل إجراء عملية واحدة هنا دون إكمال العملية السابقة. وهذا يعني أننا لن نكون قادرين على استخدام القدرة على القيام بالمضاعفات المزدوجة ، كما يسمح لنا خط أنابيب SIMD بذلك. إذا قمنا بتحويل الحسابات
fma قليلاً (a، b، c * d) + fma (e، f، g * h) ، فسوف نرى أنه يمكننا موازاة الحسابات. أولاً قم بعمل
مضاعفتين مستقلتين ، ثم عمليتي
fma مستقلتين.
تنفيذ AVX + FMA __m256 fma(__m256 const a, __m256 const b, __m256 const c) { return _mm256_fmadd_ps(a, b, c); } void mul_mtx4_mtx4_avx_fma(__m128* const r, __m128 const* const m, __m128 const* const n) { __m256 const mm[]{ _mm256_broadcast_ps(m + 0), _mm256_broadcast_ps(m + 1), _mm256_broadcast_ps(m + 2), _mm256_broadcast_ps(m + 3) }; __m256 const n0n1 = _mm256_load_ps(&n[0].m128_f32[0]); _mm256_stream_ps(&r[0].m128_f32[0], fma(perm<3>(n0n1), mm[0], perm<2>(n0n1)*mm[1])+ fma(perm<1>(n0n1), mm[2], perm<0>(n0n1)*mm[3])); __m256 const n2n3 = _mm256_load_ps(&n[2].m128_f32[0]); _mm256_stream_ps(&r[2].m128_f32[0], fma(perm<3>(n2n3), mm[0], perm<2>(n2n3)*mm[1])+ fma(perm<1>(n2n3), mm[2], perm<0>(n2n3)*mm[3])); }
IACA:
x86 - 9.21 ، x64 - 8 . الآن جيد جدًا. ربما سيقول شخص ما ما يمكن القيام به بشكل أفضل ، لكني لا أعرف كيف.
المعايير
ألاحظ على الفور أن هذه الأرقام لا ينبغي اعتبارها الحقيقة المطلقة. حتى مع اختبار ثابت ، يسبحون في حدود معينة.
بل وأكثر من ذلك يتصرفون بشكل مختلف على منصات مختلفة. مع أي تحسين ، خذ قياسات خاصة لحالتك.جدول المحتويات
- الوظيفة: اسم الوظيفة. تعمل النهاية على s - مع دفق ، حركة عادية مختلفة (بدون دفق). أضيفت للوضوح ، لأن هذا مهم بما فيه الكفاية.
- دورات IACA: عدد علامات التجزئة لكل وظيفة محسوبة بواسطة IACA
- الدورات المقاسة: عدد المقاييس المقاسة (أقل هو أكثر)
- تسريع IACA: عدد القياسات في خط صفري / عدد القياسات في خط
- التسارع المقاس: عدد القياسات في خط الصفر / عدد القياسات في الخط (كلما كان ذلك أفضل)
بالنسبة إلى loop_m ، تم ضرب القراد من المقالة في 64. أي أن هذه قيمة تقريبية جدًا. في الواقع ، اتضح ذلك.i3-3770:
i7-8700 ك:
كود الاختبار في المصدر. إذا كانت هناك اقتراحات معقولة حول كيفية تحسينها ، فاكتب في التعليقات.مكافأة من عالم الخيال
في الواقع ، إنه من عالم الخيال لأنه إذا رأيت معالجات تدعم AVX512 ، فربما في الصور. ومع ذلك ، حاولت تنفيذ الخوارزمية. هنا لن أشرح أي شيء ، تشبيهًا كاملاً مع AVX + FMA. الخوارزمية هي نفسها ، عمليات أقل فقط.كما يقولون ، سأترك الأمر هنا __m512 operator + (__m512 const a, __m512 const b) { return _mm512_add_ps(a, b); } __m512 operator - (__m512 const a, __m512 const b) { return _mm512_sub_ps(a, b); } __m512 operator * (__m512 const a, __m512 const b) { return _mm512_mul_ps(a, b); } __m512 operator / (__m512 const a, __m512 const b) { return _mm512_div_ps(a, b); } template <int i> __m512 perm(__m512 const v) { return _mm512_permute_ps(v, _MM_SHUFFLE(i, i, i, i)); } template <int a, int b, int c, int d> __m512 perm(__m512 const v) { return _mm512_permute_ps(v, _MM_SHUFFLE(a, b, c, d)); } __m512 fma(__m512 const a, __m512 const b, __m512 const c) { return _mm512_fmadd_ps(a, b, c); } void mul_mtx4_mtx4_avx512(__m128* const r, __m128 const* const m, __m128 const* const _n) { __m512 const mm[]{ _mm512_broadcast_f32x4(m[0]), _mm512_broadcast_f32x4(m[1]), _mm512_broadcast_f32x4(m[2]), _mm512_broadcast_f32x4(m[3]) }; __m512 const n = _mm512_load_ps(&_n[0].m128_f32[0]); _mm512_stream_ps(&r[0].m128_f32[0], fma(perm<3>(n), mm[0], perm<2>(n)*mm[1])+ fma(perm<1>(n), mm[2], perm<0>(n)*mm[3])); }
الأرقام رائعة: x86 - 4.79 ، x64 - 5.42 (IACA مع بنية SKX). هذا على الرغم من حقيقة أن الخوارزمية لديها 64 مضاعفات و 48 إضافات.كود PS من المقالة
هذه تجربتي الأولى في كتابة مقال. شكرا لكم جميعا على تعليقاتكم. يساعدون في جعل الشفرة والمقال أفضل.