أصناف من SIMD

أثناء تطوير meshoptimizer ، يطرح السؤال غالبًا: "هل يمكن لهذه الخوارزمية استخدام SIMD؟"

المكتبة موجهة نحو الأداء ، لكن SIMD لا توفر دائمًا مزايا كبيرة للسرعة. لسوء الحظ ، يمكن أن تجعل SIMD الشفرة أقل قدرة على الحركة وأقل قابلية للصيانة. لذلك ، في كل حالة ، من الضروري البحث عن حل وسط. عندما يكون الأداء قصوى ، يجب عليك تطوير وصيانة تطبيقات SIMD منفصلة لمجموعات تعليمات SSE و NEON. في حالات أخرى ، تحتاج إلى فهم تأثير استخدام SIMD. سنحاول اليوم تسريع شبكة المبسطة ، وهي خوارزمية جديدة تمت إضافتها مؤخرًا إلى المكتبة ، باستخدام مجموعات التعليمات SSEn / AVXn.



بالنسبة لمعيارنا ، نقوم بتبسيط نموذج بوذا التايلاندي من 6 ملايين مثلث إلى 0.1٪ من هذا الرقم. نستخدم برنامج التحويل البرمجي Microsoft Visual Studio 2019 لبنية x64 الهدف. يمكن للخوارزمية العددية إجراء مثل هذا الترشيد في حوالي 210 مللي ثانية في مؤشر ترابط Intel Core i7-8700K واحد (عند حوالي 4.4 جيجاهرتز). هذا يتوافق مع حوالي 28.5 مليون مثلث في الثانية. ربما هذا يكفي في الممارسة العملية ، لكنني كنت فضولي لاستكشاف أقصى قدرات الجهاز.

في بعض الحالات ، يمكن مزامنة الإجراء بتقسيم الشبكة إلى أجزاء ، ولكن من الضروري إجراء تحليل إضافي للتوصيل من أجل الحفاظ على الحدود ، لذلك سنقوم في الوقت الحالي بتقييد أنفسنا بتحسينات SIMD البحتة.

سبعة قياسات


لفهم إمكانات التحسين ، سنقوم بإجراء التوصيف باستخدام Intel VTune. قم بتشغيل الإجراء 100 مرة للتأكد من وجود بيانات كافية.



قمت هنا بتشغيل وضع microarchitecture لإصلاح وقت تنفيذ كل وظيفة ، وكذلك العثور على الاختناقات. نرى أن الترشيد يتم تنفيذه باستخدام مجموعة من الوظائف ، كل منها يتطلب عددًا معينًا من الدورات. يتم فرز قائمة الوظائف حسب الوقت. ها هم من أجل التنفيذ لجعل الخوارزمية أسهل للفهم:

  • rescalePositions مواضع كل القمم في مكعب واحد للتحضير computeVertexIds باستخدام computeVertexIds
  • computeVertexIds يحسب معرّف كمي 30 بت لكل رأس على شبكة موحدة بحجم معين ، حيث يتم قياس كل محور على شبكة (حجم الشبكة 10 بت ، وبالتالي فإن المعرف هو 30)
  • countTriangles بحساب العدد التقريبي للمثلثات التي سيقوم المبتكر بإنشائها بحجم شبكة معين ، مع افتراض اتحاد كل القمم في خلية شبكة واحدة
  • fillVertexCells تملأ جدولًا fillVertexCells كل الرؤوس للخلايا المقابلة ؛ جميع القمم ذات المعرف نفسه تتوافق مع خلية واحدة
  • fillCellQuadrics تملأ هيكل Quadric (مصفوفة متماثلة Quadric ) لكل خلية لتعكس المعلومات الإجمالية عن الهندسة المقابلة
  • fillCellRemap بحساب مؤشر قمة الرأس لكل خلية ، واختيار واحد من القمم في هذه الخلية ، ويقلل من التشوه الهندسي
  • يعرض filterTriangles المجموعة النهائية من المثلثات وفقًا لجداول vertex-cell-vertex التي تم إنشاؤها مسبقًا ؛ يمكن أن ينتج عن الترجمة الساذجة معدل مثلثات مكررة بنسبة 5٪ تقريبًا ، وبالتالي تقوم الدالة بتصفية التكرارات.

يتم تنفيذ computeVertexIds و countTriangles عدة مرات: تحدد الخوارزمية حجم الشبكة لدمج القمم ، وإجراء بحث ثنائي معجل لتحقيق العدد المستهدف من المثلثات (في الحالة 6000 لدينا) وحساب عدد المثلثات التي ستنشئها كل حجم شبكة في كل تكرار. يتم إطلاق وظائف أخرى مرة واحدة. في ملفنا ، تعطي خمس ممرات بحث حجم شبكة مستهدف قدره 40 3 .

تفيد VTune بشكل مفيد أن الوظيفة الأكثر كثافة للموارد هي تلك التي تقوم بحساب الترميزات: فهي تستغرق ما يقرب من نصف إجمالي وقت التنفيذ البالغ 21 ثانية. هذا هو الهدف الأول لتحسين SIMD.

SIMD قطعة قطعة


دعنا نلقي نظرة على الكود المصدري لـ fillCellQuadrics لفهم ما يحسبه بالضبط:

 static void fillCellQuadrics(Quadric* cell_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* vertex_cells) { for (size_t i = 0; i < index_count; i += 3) { unsigned int i0 = indices[i + 0]; unsigned int i1 = indices[i + 1]; unsigned int i2 = indices[i + 2]; unsigned int c0 = vertex_cells[i0]; unsigned int c1 = vertex_cells[i1]; unsigned int c2 = vertex_cells[i2]; bool single_cell = (c0 == c1) & (c0 == c2); float weight = single_cell ? 3.f : 1.f; Quadric Q; quadricFromTriangle(Q, vertex_positions[i0], vertex_positions[i1], vertex_positions[i2], weight); if (single_cell) { quadricAdd(cell_quadrics[c0], Q); } else { quadricAdd(cell_quadrics[c0], Q); quadricAdd(cell_quadrics[c1], Q); quadricAdd(cell_quadrics[c2], Q); } } } 

تتكرر الوظيفة على كل المثلثات ، وتحسب المربعي لكل منها ، وتضيفها إلى مربعات كل خلية. المربعي - مصفوفة متماثلة 4 × 4 ، معروضة كـ 10 أرقام فاصلة عائمة:

 struct Quadric { float a00; float a10, a11; float a20, a21, a22; float b0, b1, b2, c; }; 

يتطلب حساب المربعي حل معادلة المستوى للمثلث ، وإنشاء مصفوفة من الدرجة الثانية ووزنها باستخدام منطقة المثلث:

 static void quadricFromPlane(Quadric& Q, float a, float b, float c, float d) { Q.a00 = a * a; Q.a10 = b * a; Q.a11 = b * b; Q.a20 = c * a; Q.a21 = c * b; Q.a22 = c * c; Q.b0 = d * a; Q.b1 = d * b; Q.b2 = d * c; Qc = d * d; } static void quadricFromTriangle(Quadric& Q, const Vector3& p0, const Vector3& p1, const Vector3& p2, float weight) { Vector3 p10 = {p1.x - p0.x, p1.y - p0.y, p1.z - p0.z}; Vector3 p20 = {p2.x - p0.x, p2.y - p0.y, p2.z - p0.z}; Vector3 normal = { p10.y * p20.z - p10.z * p20.y, p10.z * p20.x - p10.x * p20.z, p10.x * p20.y - p10.y * p20.x }; float area = normalize(normal); float distance = normal.x*p0.x + normal.y*p0.y + normal.z*p0.z; quadricFromPlane(Q, normal.x, normal.y, normal.z, -distance); quadricMul(Q, area * weight); } 

يبدو أن هناك الكثير من عمليات الفاصلة العائمة ، بحيث يمكن موازنتها باستخدام SIMD. أولاً ، نحن نمثل كل متجه باعتباره متجه SIMD عريض النطاق ، Quadric أيضًا هيكل Quadric إلى 12 رقمًا عائمًا بدلاً من 10 بحيث يتناسب تمامًا مع ثلاث سجلات SIMD (زيادة الحجم لا يؤثر على الأداء) وتغيير ترتيب الحقول بحيث الحسابات أصبح quadricFromPlane أكثر اتساقا:

 struct Quadric { float a00, a11, a22; float pad0; float a10, a21, a20; float pad1; float b0, b1, b2, c; }; 

هنا ، لا تتسق بعض العمليات الحسابية ، ولا سيما المنتج القياسي ، مع الإصدارات السابقة من SSE. لحسن الحظ ، ظهرت تعليمات لمنتج عددي في SSE4.1 ، وهو أمر مفيد للغاية بالنسبة لنا.

 static void fillCellQuadrics(Quadric* cell_quadrics, const unsigned int* indices, size_t index_count, const Vector3* vertex_positions, const unsigned int* vertex_cells) { const int yzx = _MM_SHUFFLE(3, 0, 2, 1); const int zxy = _MM_SHUFFLE(3, 1, 0, 2); const int dp_xyz = 0x7f; for (size_t i = 0; i < index_count; i += 3) { unsigned int i0 = indices[i + 0]; unsigned int i1 = indices[i + 1]; unsigned int i2 = indices[i + 2]; unsigned int c0 = vertex_cells[i0]; unsigned int c1 = vertex_cells[i1]; unsigned int c2 = vertex_cells[i2]; bool single_cell = (c0 == c1) & (c0 == c2); __m128 p0 = _mm_loadu_ps(&vertex_positions[i0].x); __m128 p1 = _mm_loadu_ps(&vertex_positions[i1].x); __m128 p2 = _mm_loadu_ps(&vertex_positions[i2].x); __m128 p10 = _mm_sub_ps(p1, p0); __m128 p20 = _mm_sub_ps(p2, p0); __m128 normal = _mm_sub_ps( _mm_mul_ps( _mm_shuffle_ps(p10, p10, yzx), _mm_shuffle_ps(p20, p20, zxy)), _mm_mul_ps( _mm_shuffle_ps(p10, p10, zxy), _mm_shuffle_ps(p20, p20, yzx))); __m128 areasq = _mm_dp_ps(normal, normal, dp_xyz); // SSE4.1 __m128 area = _mm_sqrt_ps(areasq); // masks the result of the division when area==0 // scalar version does this in normalize() normal = _mm_and_ps( _mm_div_ps(normal, area), _mm_cmpneq_ps(area, _mm_setzero_ps())); __m128 distance = _mm_dp_ps(normal, p0, dp_xyz); // SSE4.1 __m128 negdistance = _mm_sub_ps(_mm_setzero_ps(), distance); __m128 normalnegdist = _mm_blend_ps(normal, negdistance, 8); __m128 Qx = _mm_mul_ps(normal, normal); __m128 Qy = _mm_mul_ps( _mm_shuffle_ps(normal, normal, _MM_SHUFFLE(3, 2, 2, 1)), _mm_shuffle_ps(normal, normal, _MM_SHUFFLE(3, 0, 1, 0))); __m128 Qz = _mm_mul_ps(negdistance, normalnegdist); if (single_cell) { area = _mm_mul_ps(area, _mm_set1_ps(3.f)); Qx = _mm_mul_ps(Qx, area); Qy = _mm_mul_ps(Qy, area); Qz = _mm_mul_ps(Qz, area); Quadric& q0 = cell_quadrics[c0]; __m128 q0x = _mm_loadu_ps(&q0.a00); __m128 q0y = _mm_loadu_ps(&q0.a10); __m128 q0z = _mm_loadu_ps(&q0.b0); _mm_storeu_ps(&q0.a00, _mm_add_ps(q0x, Qx)); _mm_storeu_ps(&q0.a10, _mm_add_ps(q0y, Qy)); _mm_storeu_ps(&q0.b0, _mm_add_ps(q0z, Qz)); } else { // omitted for brevity, repeats the if() body // three times for c0/c1/c2 } } } 

لا يوجد شيء مثير للاهتمام بشكل خاص في هذا الرمز ؛ نحن نستخدم بكثرة تعليمات التحميل / المتجر غير المحاذاة. على الرغم من أنه يمكن محاذاة مدخلات Vector3 ، يبدو أنه لا توجد عقوبة ملحوظة للقراءات غير المحاذاة. يرجى ملاحظة أنه في النصف الأول من ناقلات الوظائف لا يتم استخدامها ، وهو أمر جيد - تحتوي ناقلاتنا على ثلاثة مكونات ، وفي بعض الحالات واحد فقط (راجع حساب المساحات / المساحة / المسافة) ، بينما يقوم المعالج بإجراء 4 عمليات على التوازي. في أي حال ، دعونا نرى كيف ساعد التوازي.



يتم تشغيل مائة بداية من fillCellQuadrics في 5.3 ثوانٍ بدلاً من 9.8 ثانية ، مما يوفر حوالي 45 مللي ثانية في كل عملية - ليس سيئًا ، ولكن ليس مؤثرًا جدًا. في العديد من الإرشادات ، نستخدم ثلاثة مكونات بدلاً من أربعة ، بالإضافة إلى الضرب الدقيق ، مما يعطي تأخيرًا كبيرًا إلى حد ما. إذا سبق لك أن كتبت تعليمات SIMD ، فأنت تعلم كيفية القيام بالمنتج القياسي بشكل صحيح.

للقيام بذلك ، يجب عليك القيام بأربعة متجهات في وقت واحد. بدلاً من تخزين متجه واحد كامل في سجل SIMD واحد ، نستخدم ثلاث سجلات - في أحدنا نخزن أربعة مكونات لـ x ، والآخر سنخزن ، وفي z الثالث. هناك حاجة إلى أربعة متجهات للعمل في آن واحد: هذا يعني أننا سنقوم بمعالجة أربعة مثلثات في وقت واحد.

لدينا العديد من المصفوفات مع فهرسة ديناميكية. عادةً ما يساعد على نقل البيانات إلى صفائف محضرة من مكونات x / y / z (أو بالأحرى ، عادةً ما يتم استخدام سجلات SIMD الصغيرة ، على سبيل المثال ، float x[8], y[8], z[8] ، لكل من الرؤوس الثمانية في الإدخال البيانات: يطلق عليها AoSoA (صفائف بنيات الصفيف) وتعطي توازنًا جيدًا بين مكان ذاكرة التخزين المؤقت وسهولة التحميل في سجلات SIMD) ، لكن الفهرسة الديناميكية هنا لن تعمل بشكل جيد للغاية ، لذا قم بتحميل البيانات لأربعة مثلثات كالمعتاد ، وقم بنقل المتجهات باستخدام طريقة ملائمة الماكرو _MM_TRANSPOSE .

من الناحية النظرية ، تحتاج إلى حساب كل مكون من أربعة محارف منتهية في سجل SIMD الخاص به (على سبيل المثال ، سيكون لدينا __m128 Q_a00 مع أربعة مكونات من a00 المحددة). في هذه الحالة ، تتلاءم العمليات على الكوادارات تمامًا مع تعليمات SIMD عريضة النطاق ، ويبطئ التحويل فعليًا من الكود - وبالتالي ، فإننا نغير المتجهات الأولية فقط ، ثم نعيد معادلات المستوى إلى الوراء ونعمل نفس الكود الذي استخدمناه لحساب الكواداركيات ، لكن نكررها اربع مرات فيما يلي الكود الذي يحسب بعد ذلك معادلات المستوى (يتم حذف الأجزاء المتبقية للإيجاز):

 unsigned int i00 = indices[(i + 0) * 3 + 0]; unsigned int i01 = indices[(i + 0) * 3 + 1]; unsigned int i02 = indices[(i + 0) * 3 + 2]; unsigned int i10 = indices[(i + 1) * 3 + 0]; unsigned int i11 = indices[(i + 1) * 3 + 1]; unsigned int i12 = indices[(i + 1) * 3 + 2]; unsigned int i20 = indices[(i + 2) * 3 + 0]; unsigned int i21 = indices[(i + 2) * 3 + 1]; unsigned int i22 = indices[(i + 2) * 3 + 2]; unsigned int i30 = indices[(i + 3) * 3 + 0]; unsigned int i31 = indices[(i + 3) * 3 + 1]; unsigned int i32 = indices[(i + 3) * 3 + 2]; // load first vertex of each triangle and transpose into vectors with components (pw0 isn't used later) __m128 px0 = _mm_loadu_ps(&vertex_positions[i00].x); __m128 py0 = _mm_loadu_ps(&vertex_positions[i10].x); __m128 pz0 = _mm_loadu_ps(&vertex_positions[i20].x); __m128 pw0 = _mm_loadu_ps(&vertex_positions[i30].x); _MM_TRANSPOSE4_PS(px0, py0, pz0, pw0); // load second vertex of each triangle and transpose into vectors with components __m128 px1 = _mm_loadu_ps(&vertex_positions[i01].x); __m128 py1 = _mm_loadu_ps(&vertex_positions[i11].x); __m128 pz1 = _mm_loadu_ps(&vertex_positions[i21].x); __m128 pw1 = _mm_loadu_ps(&vertex_positions[i31].x); _MM_TRANSPOSE4_PS(px1, py1, pz1, pw1); // load third vertex of each triangle and transpose into vectors with components __m128 px2 = _mm_loadu_ps(&vertex_positions[i02].x); __m128 py2 = _mm_loadu_ps(&vertex_positions[i12].x); __m128 pz2 = _mm_loadu_ps(&vertex_positions[i22].x); __m128 pw2 = _mm_loadu_ps(&vertex_positions[i32].x); _MM_TRANSPOSE4_PS(px2, py2, pz2, pw2); // p1 - p0 __m128 px10 = _mm_sub_ps(px1, px0); __m128 py10 = _mm_sub_ps(py1, py0); __m128 pz10 = _mm_sub_ps(pz1, pz0); // p2 - p0 __m128 px20 = _mm_sub_ps(px2, px0); __m128 py20 = _mm_sub_ps(py2, py0); __m128 pz20 = _mm_sub_ps(pz2, pz0); // cross(p10, p20) __m128 normalx = _mm_sub_ps( _mm_mul_ps(py10, pz20), _mm_mul_ps(pz10, py20)); __m128 normaly = _mm_sub_ps( _mm_mul_ps(pz10, px20), _mm_mul_ps(px10, pz20)); __m128 normalz = _mm_sub_ps( _mm_mul_ps(px10, py20), _mm_mul_ps(py10, px20)); // normalize; note that areasq/area now contain 4 values, not just one __m128 areasq = _mm_add_ps( _mm_mul_ps(normalx, normalx), _mm_add_ps( _mm_mul_ps(normaly, normaly), _mm_mul_ps(normalz, normalz))); __m128 area = _mm_sqrt_ps(areasq); __m128 areanz = _mm_cmpneq_ps(area, _mm_setzero_ps()); normalx = _mm_and_ps(_mm_div_ps(normalx, area), areanz); normaly = _mm_and_ps(_mm_div_ps(normaly, area), areanz); normalz = _mm_and_ps(_mm_div_ps(normalz, area), areanz); __m128 distance = _mm_add_ps( _mm_mul_ps(normalx, px0), _mm_add_ps( _mm_mul_ps(normaly, py0), _mm_mul_ps(normalz, pz0))); __m128 negdistance = _mm_sub_ps(_mm_setzero_ps(), distance); // this computes the plane equations (a, b, c, d) for each of the 4 triangles __m128 plane0 = normalx; __m128 plane1 = normaly; __m128 plane2 = normalz; __m128 plane3 = negdistance; _MM_TRANSPOSE4_PS(plane0, plane1, plane2, plane3); 

أصبح الكود أطول قليلاً: نحن الآن نعالج أربعة مثلثات في كل تكرار ، ولم نعد بحاجة إلى تعليمات SSE4.1 لهذا الغرض. من الناحية النظرية ، ينبغي استخدام وحدات SIMD بشكل أكثر كفاءة. دعونا نرى كيف ساعد.



حسنا ، هذا جيد. تسارعت الشفرة بشكل طفيف للغاية ، على الرغم من أن وظيفة fillCellQuadrics تعمل الآن بمعدل أسرع مرتين تقريبًا من الوظيفة الأصلية بدون SIMD ، لكن من غير الواضح ما إذا كان هذا يبرر زيادة كبيرة في التعقيد. من الناحية النظرية ، يمكنك استخدام AVX2 ومعالجة 8 مثلثات لكل تكرار ، ولكن هنا ستحتاج إلى زيادة دوران الحلقة يدويًا (من الناحية المثالية ، يتم إنشاء كل هذا الرمز باستخدام ISPC ، ولكن محاولاتي الساذجة للحصول عليه لإنشاء رمز جيد لم تنجح: بدلاً من تسلسل التحميل / التخزين أصدر بإصرار مجموعة / مبعثر ، مما يؤدي إلى إبطاء التنفيذ بشكل كبير). لنجرب شيئًا آخر.

AVX2 = SSE2 + SSE2


AVX2 عبارة عن مجموعة من التعليمات. لديها سجلات عائمة بنطاق عريض 8 ، وتعليمات واحدة يمكن أن تؤدي 8 عمليات ؛ لكن في جوهرها ، لا يختلف مثل هذا التعليم عن إرشادات SSE2 التي يتم تنفيذها على نصفين من السجل (على حد علمي ، فك شفرة المعالجات الأولى التي تدعم AVX2 كل تعليمة في اثنين أو أكثر من عمليات التشغيل المصغرة ، وبالتالي فإن أداء الأداء كان محدودًا بمرحلة استخراج التعليمات). على سبيل المثال ، ينفذ _mm_dp_ps منتجًا عدديًا بين _mm256_dp_ps SSE2 ، وينتج عن _mm256_dp_ps بين نصفي _mm256_dp_ps ، وبالتالي فهو يقتصر على 4 عرض لكل نصف.

لهذا السبب ، غالبًا ما يختلف رمز AVX2 عن "SIMD عريض النطاق" الشامل ، ولكن هنا يعمل لصالحنا: بدلاً من محاولة تحسين التوازي عن طريق نقل المتجهات ذات الأربعة أحجام ، نعود إلى الإصدار الأول من SIMD ونضاعف الحلقة باستخدام تعليمات AVX2 بدلاً من SSE2 / SSE4. ما زلنا بحاجة إلى تحميل وتخزين المتجهات ذات 4 __m128 ، ولكن بشكل عام ، نقوم فقط بتغيير __m128 إلى __m256 و _mm_ إلى _mm256 في _mm256 مع العديد من الإعدادات:

 unsigned int i00 = indices[(i + 0) * 3 + 0]; unsigned int i01 = indices[(i + 0) * 3 + 1]; unsigned int i02 = indices[(i + 0) * 3 + 2]; unsigned int i10 = indices[(i + 1) * 3 + 0]; unsigned int i11 = indices[(i + 1) * 3 + 1]; unsigned int i12 = indices[(i + 1) * 3 + 2]; __m256 p0 = _mm256_loadu2_m128( &vertex_positions[i10].x, &vertex_positions[i00].x); __m256 p1 = _mm256_loadu2_m128( &vertex_positions[i11].x, &vertex_positions[i01].x); __m256 p2 = _mm256_loadu2_m128( &vertex_positions[i12].x, &vertex_positions[i02].x); __m256 p10 = _mm256_sub_ps(p1, p0); __m256 p20 = _mm256_sub_ps(p2, p0); __m256 normal = _mm256_sub_ps( _mm256_mul_ps( _mm256_shuffle_ps(p10, p10, yzx), _mm256_shuffle_ps(p20, p20, zxy)), _mm256_mul_ps( _mm256_shuffle_ps(p10, p10, zxy), _mm256_shuffle_ps(p20, p20, yzx))); __m256 areasq = _mm256_dp_ps(normal, normal, dp_xyz); __m256 area = _mm256_sqrt_ps(areasq); __m256 areanz = _mm256_cmp_ps(area, _mm256_setzero_ps(), _CMP_NEQ_OQ); normal = _mm256_and_ps(_mm256_div_ps(normal, area), areanz); __m256 distance = _mm256_dp_ps(normal, p0, dp_xyz); __m256 negdistance = _mm256_sub_ps(_mm256_setzero_ps(), distance); __m256 normalnegdist = _mm256_blend_ps(normal, negdistance, 0x88); __m256 Qx = _mm256_mul_ps(normal, normal); __m256 Qy = _mm256_mul_ps( _mm256_shuffle_ps(normal, normal, _MM_SHUFFLE(3, 2, 2, 1)), _mm256_shuffle_ps(normal, normal, _MM_SHUFFLE(3, 0, 1, 0))); __m256 Qz = _mm256_mul_ps(negdistance, normalnegdist); 

هنا يمكنك أن تأخذ كل نصف 128-bit من Qz Qx / Qy / Qz المستلمة وتشغيل نفس الكود الذي استخدمناه لإضافة الكواداليك. بدلاً من ذلك ، نفترض أنه إذا كان للمثلث ثلاثة رؤوس في خلية واحدة ( single_cell == true ) ، فمن المحتمل جدًا أن يكون لدى مثلث آخر ثلاثة رؤوس في خلية أخرى ، ونقوم بإجراء التجميع النهائي للرباعي باستخدام AVX2 أيضًا:

 unsigned int c00 = vertex_cells[i00]; unsigned int c01 = vertex_cells[i01]; unsigned int c02 = vertex_cells[i02]; unsigned int c10 = vertex_cells[i10]; unsigned int c11 = vertex_cells[i11]; unsigned int c12 = vertex_cells[i12]; bool single_cell = (c00 == c01) & (c00 == c02) & (c10 == c11) & (c10 == c12); if (single_cell) { area = _mm256_mul_ps(area, _mm256_set1_ps(3.f)); Qx = _mm256_mul_ps(Qx, area); Qy = _mm256_mul_ps(Qy, area); Qz = _mm256_mul_ps(Qz, area); Quadric& q00 = cell_quadrics[c00]; Quadric& q10 = cell_quadrics[c10]; __m256 q0x = _mm256_loadu2_m128(&q10.a00, &q00.a00); __m256 q0y = _mm256_loadu2_m128(&q10.a10, &q00.a10); __m256 q0z = _mm256_loadu2_m128(&q10.b0, &q00.b0); _mm256_storeu2_m128(&q10.a00, &q00.a00, _mm256_add_ps(q0x, Qx)); _mm256_storeu2_m128(&q10.a10, &q00.a10, _mm256_add_ps(q0y, Qy)); _mm256_storeu2_m128(&q10.b0, &q00.b0, _mm256_add_ps(q0z, Qz)); } else { // omitted for brevity } 

الشفرة الناتجة أبسط وموجزة وأسرع من طريقة SSE2 غير الناجحة:



بالطبع ، لم نحقق تسارع 8 مرات ، ولكن 2.45 مرة فقط. لا تزال عمليات التحميل والتخزين واسعة النطاق ، حيث إننا مضطرون للعمل مع تصميم ذاكرة غير مريح بسبب الفهرسة الديناميكية ، والحسابات ليست مثالية لـ SIMD. ولكن الآن لم تعد fillCellQuadrics عنق الزجاجة في خط أنابيب الملف الشخصي لدينا ، ويمكنك التركيز على وظائف أخرى.

جمع حولها ، والأطفال


لقد أنقذنا 4.8 ثانية من التشغيل التجريبي (48 مللي ثانية في كل تشغيل) ، والآن أصبح الدخيل الرئيسي هو countTriangles . يبدو أن هذه وظيفة بسيطة ، ولكنها لا تنفذ مرة واحدة ، بل خمس مرات:

 static size_t countTriangles(const unsigned int* vertex_ids, const unsigned int* indices, size_t index_count) { size_t result = 0; for (size_t i = 0; i < index_count; i += 3) { unsigned int id0 = vertex_ids[indices[i + 0]]; unsigned int id1 = vertex_ids[indices[i + 1]]; unsigned int id2 = vertex_ids[indices[i + 2]]; result += (id0 != id1) & (id0 != id2) & (id1 != id2); } return result; } 

تتكرر الوظيفة على كل المثلثات المصدر وتحسب عدد المثلثات غير المنحطة من خلال مقارنة معرفات الرؤوس. ليس من الواضح على الفور كيفية موازاة ذلك باستخدام SIMD ... إلا إذا كنت تستخدم إرشادات التجميع.

أضافت مجموعة تعليمات AVX2 عائلة إرشادات التجميع / الانتثار إلى x64 SIMD ؛ يأخذ كل منهم سجلًا متجهًا يحتوي على 4 أو 8 قيم - وفي نفس الوقت ينفذ عمليات تحميل أو حفظ 4 أو 8. إذا كنت تستخدم التجميع هنا ، فيمكنك تنزيل ثلاثة فهارس وتشغيل التجميع للجميع مرة واحدة (أو في مجموعات من 4 أو 8) ومقارنة النتائج. لقد كان تجميع معالجات Intel بطيئًا من الناحية التاريخية ، ولكن دعونا نجربها. للبساطة ، نقوم بتحميل البيانات الخاصة بـ 8 مثلثات ، ونحول المتجهات بشكل مشابه لمحاولتنا السابقة ، ونقارن العناصر المقابلة لكل متجه:

 for (size_t i = 0; i < (triangle_count & ~7); i += 8) { __m256 tri0 = _mm256_loadu2_m128( (const float*)&indices[(i + 4) * 3 + 0], (const float*)&indices[(i + 0) * 3 + 0]); __m256 tri1 = _mm256_loadu2_m128( (const float*)&indices[(i + 5) * 3 + 0], (const float*)&indices[(i + 1) * 3 + 0]); __m256 tri2 = _mm256_loadu2_m128( (const float*)&indices[(i + 6) * 3 + 0], (const float*)&indices[(i + 2) * 3 + 0]); __m256 tri3 = _mm256_loadu2_m128( (const float*)&indices[(i + 7) * 3 + 0], (const float*)&indices[(i + 3) * 3 + 0]); _MM_TRANSPOSE8_LANE4_PS(tri0, tri1, tri2, tri3); __m256i i0 = _mm256_castps_si256(tri0); __m256i i1 = _mm256_castps_si256(tri1); __m256i i2 = _mm256_castps_si256(tri2); __m256i id0 = _mm256_i32gather_epi32((int*)vertex_ids, i0, 4); __m256i id1 = _mm256_i32gather_epi32((int*)vertex_ids, i1, 4); __m256i id2 = _mm256_i32gather_epi32((int*)vertex_ids, i2, 4); __m256i deg = _mm256_or_si256( _mm256_cmpeq_epi32(id1, id2), _mm256_or_si256( _mm256_cmpeq_epi32(id0, id1), _mm256_cmpeq_epi32(id0, id2))); result += 8 - _mm_popcnt_u32(_mm256_movemask_epi8(deg)) / 4; } 

ماكرو _MM_TRANSPOSE8_LANE4_PS في _MM_TRANSPOSE8_LANE4_PS هو ما يعادل _MM_TRANSPOSE4_PS ، وهو غير موجود في الرأس القياسي ، ولكن يمكن عرضه بسهولة. يستغرق أربعة ناقلات AVX2 وينقل اثنين من المصفوفات 4 × 4 بشكل مستقل عن بعضها البعض:

 #define _MM_TRANSPOSE8_LANE4_PS(row0, row1, row2, row3) \ do { \ __m256 __t0, __t1, __t2, __t3; \ __t0 = _mm256_unpacklo_ps(row0, row1); \ __t1 = _mm256_unpackhi_ps(row0, row1); \ __t2 = _mm256_unpacklo_ps(row2, row3); \ __t3 = _mm256_unpackhi_ps(row2, row3); \ row0 = _mm256_shuffle_ps(__t0, __t2, _MM_SHUFFLE(1, 0, 1, 0)); \ row1 = _mm256_shuffle_ps(__t0, __t2, _MM_SHUFFLE(3, 2, 3, 2)); \ row2 = _mm256_shuffle_ps(__t1, __t3, _MM_SHUFFLE(1, 0, 1, 0)); \ row3 = _mm256_shuffle_ps(__t1, __t3, _MM_SHUFFLE(3, 2, 3, 2)); \ } while (0) 

نظرًا لبعض الميزات الموجودة في مجموعات التعليمات SSE2 / AVX2 ، يجب علينا استخدام عمليات تسجيل الفاصلة العائمة عند نقل المتجهات. نحن نقوم بتحميل البيانات قليلاً بلا مبالاة ؛ ولكن هذا لا يهم في الأساس ، لأن تجميع الأداء يحدنا:



countTriangles الآن حوالي 27٪ ويلاحظ مؤشر أسعار المستهلك الرهيب (دورات لكل تعليمة): نقوم بإرسال تعليمات أقل بمقدار أربعة أضعاف ، لكن التجميع يستغرق الكثير من الوقت. من الرائع تسريع العمل بشكل عام ، ولكن بالطبع ، فإن اكتساب الأداء محبط إلى حد ما. لقد نجحنا في تجاوز fillCellQuadrics في الملف الشخصي ، وهو ما يقودنا إلى الوظيفة الأخيرة في أعلى القائمة ، والتي لم نطلع عليها بعد.

الفصل 6 ، حيث يبدأ كل شيء في العمل كما ينبغي


آخر وظيفة هي computeVertexIds . في الخوارزمية ، يتم تنفيذها 6 مرات ، وبالتالي فهي تمثل أيضًا هدفًا ممتازًا للتحسين. لأول مرة ، نرى وظيفة يبدو أنها تم إنشاؤها لتحسين واضح في SIMD:

 static void computeVertexIds(unsigned int* vertex_ids, const Vector3* vertex_positions, size_t vertex_count, int grid_size) { assert(grid_size >= 1 && grid_size <= 1024); float cell_scale = float(grid_size - 1); for (size_t i = 0; i < vertex_count; ++i) { const Vector3& v = vertex_positions[i]; int xi = int(vx * cell_scale + 0.5f); int yi = int(vy * cell_scale + 0.5f); int zi = int(vz * cell_scale + 0.5f); vertex_ids[i] = (xi << 20) | (yi << 10) | zi; } } 

بعد التحسينات السابقة ، نعرف ما يجب القيام به: استرخ الدورة 4 أو 8 مرات ، حيث لا يوجد أي فائدة في محاولة تسريع تكرار واحد فقط ، ونقل مكونات المتجه وتشغيل الحسابات بالتوازي. لنقم بذلك باستخدام AVX2 ، ونقوم بمعالجة 8 رؤوس في كل مرة:

 __m256 scale = _mm256_set1_ps(cell_scale); __m256 half = _mm256_set1_ps(0.5f); for (size_t i = 0; i < (vertex_count & ~7); i += 8) { __m256 vx = _mm256_loadu2_m128( &vertex_positions[i + 4].x, &vertex_positions[i + 0].x); __m256 vy = _mm256_loadu2_m128( &vertex_positions[i + 5].x, &vertex_positions[i + 1].x); __m256 vz = _mm256_loadu2_m128( &vertex_positions[i + 6].x, &vertex_positions[i + 2].x); __m256 vw = _mm256_loadu2_m128( &vertex_positions[i + 7].x, &vertex_positions[i + 3].x); _MM_TRANSPOSE8_LANE4_PS(vx, vy, vz, vw); __m256i xi = _mm256_cvttps_epi32( _mm256_add_ps(_mm256_mul_ps(vx, scale), half)); __m256i yi = _mm256_cvttps_epi32( _mm256_add_ps(_mm256_mul_ps(vy, scale), half)); __m256i zi = _mm256_cvttps_epi32( _mm256_add_ps(_mm256_mul_ps(vz, scale), half)); __m256i id = _mm256_or_si256( zi, _mm256_or_si256( _mm256_slli_epi32(xi, 20), _mm256_slli_epi32(yi, 10))); _mm256_storeu_si256((__m256i*)&vertex_ids[i], id); } 

وإلقاء نظرة على النتائج:



تسارعنا computeVertexIdsمرتين. مع الأخذ في الاعتبار جميع التحسينات ، تم تخفيض إجمالي وقت تنفيذ البرنامج إلى حوالي 120 مللي ثانية ، وهو ما يتوافق مع حساب 50 مليون مثلث في الثانية.

قد يبدو أننا لم نحقق مرة أخرى نمو الإنتاجية المتوقع: computeVertexIdsألا ينبغي أن يتسارع أكثر من مرتين بعد الموازاة؟ للإجابة على هذا السؤال ، دعونا نحاول معرفة مقدار العمل الذي تؤديه هذه الوظيفة.

computeVertexIdsيتم تنفيذه ست مرات لبدء برنامج واحد: خمس مرات أثناء البحث الثنائي ومرة ​​واحدة في النهاية لحساب المعرفات النهائية التي يتم استخدامها للمعالجة الإضافية. في كل مرة تقوم هذه الوظيفة بمعالجة 3 ملايين رأس ، وقراءة 12 بايت لكل قمة وكتابة 4 بايت.

في المجموع ، أكثر من 100 يدير المبتكر ، هذه الوظيفة معالجة 1.8 مليار رأس ، وقراءة 21 غيغابايت وإعادة كتابة 7 غيغابايت. تتطلب معالجة 28 جيجابايت في 1.46 ثانية عرض نطاق ترددي للحافلة يبلغ 19 جيجابايت / ثانية. يمكننا التحقق من عرض النطاق الترددي الذاكرة عن طريق تشغيل memcmp(block1, block2, 512 MB). والنتيجة هي 45 مللي ثانية ، أي حوالي 22 جيجابايت / ثانية على لب واحد (على الرغم من أن معيار AIDA64 يُظهر سرعات قراءة على نظامي يصل إلى 31 جيجابايت / ثانية ، لكنه يستخدم عدة مراكز). في الواقع ، لقد اقتربنا من الحد الأقصى الممكن تحقيقه من الذاكرة ، وستتطلب زيادة أخرى في الأداء تعبئة أوثق لهذه القمم بحيث تشغل أقل من 12 بايت.

الخاتمة


لقد اتخذنا خوارزمية محسنة بشكل جيد للغاية تعمل على تبسيط الشبكات الكبيرة بسرعة 28 مليون مثلث في الثانية ، واستخدمنا مجموعات التعليمات SSE و AVX لتسريعها مرتين تقريبًا ، إلى 50 مليون مثلث في الثانية. خلال هذه الرحلة ، كان علينا أن نتعلم طرقًا مختلفة لاستخدام SIMD: سجلات لتخزين المتجهات 3 - واسعة ، وتعليمات SoA ، و AVX2 لتخزين متجهين عريضين - 3 ، وجمع الإرشادات لتسريع تحميل البيانات مقارنة بالتعليمات العددية ، وأخيراً قمنا بتطبيق AVX2 مباشرة لمعالجة الدفق.

غالبًا ما لا يكون SIMD هو أفضل نقطة بداية للتحسين: لقد مر عقلاني الشبكات بعدة تكرارات من التحسين الخوارزمي و microoptimization دون استخدام إرشادات خاصة بالمنصة. ولكن في مرحلة ما ، يتم استنفاد هذه الاحتمالات ، وإذا كان الأداء حاسمًا ، فإن SIMD أداة رائعة يمكن استخدامها إذا لزم الأمر.

لست متأكدًا من أي من هذه التحسينات ستندرج في الفرع الرئيسي meshoptimizer: في النهاية ، هذه مجرد تجربة لمعرفة مدى وفيركلوكيد الشفرة دون تغيير جذري في الخوارزميات. آمل أن يكون هذا المقال مفيدًا وأن يقدم لك أفكارًا لتحسين الشفرة. المصادر النهائية لهذا المقال هنا ؛ يستند هذا العمل إلى إصدار meshoptimizer 99ab49، ويتم نشر نموذج بوذا التايلاندي على Sketchfab .

Source: https://habr.com/ru/post/ar441536/


All Articles