मेषोप्टाइज़र के विकास के दौरान, अक्सर यह सवाल उठता है: "क्या यह एल्गोरिथ्म SIMD का उपयोग कर सकता है?"

पुस्तकालय प्रदर्शन-उन्मुख है, लेकिन SIMD हमेशा महत्वपूर्ण गति लाभ प्रदान नहीं करता है। दुर्भाग्य से, SIMD कोड को कम पोर्टेबल और कम रखरखाव योग्य बना सकता है। इसलिए, प्रत्येक मामले में, समझौता करना आवश्यक है। जब प्रदर्शन सर्वोपरि होता है, तो आपको SSE और NEON निर्देश सेट के लिए अलग-अलग SIMD कार्यान्वयन विकसित करने और बनाए रखना होगा। अन्य मामलों में, आपको यह समझने की आवश्यकता है कि SIMD का उपयोग करने का प्रभाव क्या है। आज हम SSEn / AVXn निर्देश सेट का उपयोग करते हुए हाल ही में पुस्तकालय में जोड़े गए एक नए एल्गोरिथ्म मैला ढोने वाले सरलीकरण को गति देने का प्रयास करेंगे।



हमारे बेंचमार्क के लिए, हम थाई बुद्ध मॉडल को 6 मिलियन त्रिकोण से इस संख्या के 0.1% तक सरल करते हैं। हम लक्ष्य x64 आर्किटेक्चर के लिए Microsoft Visual Studio 2019 कंपाइलर का उपयोग करते हैं। स्केलर एल्गोरिथ्म एक इंटेल कोर i7-8700K थ्रेड (~ 4.4 हर्ट्ज पर) में लगभग 210 एमएस में इस तरह के युक्तिकरण का प्रदर्शन कर सकता है। यह प्रति सेकंड लगभग 28.5 मिलियन त्रिकोण से मेल खाती है। शायद यह अभ्यास में पर्याप्त है, लेकिन मैं उपकरणों की अधिकतम क्षमताओं का पता लगाने के लिए उत्सुक था।

कुछ मामलों में, ग्रिड को टुकड़ों में विभाजित करके प्रक्रिया को समानांतर किया जा सकता है, लेकिन इसके लिए सीमाओं को बनाए रखने के लिए कनेक्टिविटी का एक अतिरिक्त विश्लेषण करना आवश्यक है, इसलिए अब हम खुद को विशुद्ध रूप से SIMD अनुकूलन के लिए प्रतिबंधित करेंगे।

सात माप

अनुकूलन की संभावनाओं को समझने के लिए, हम Intel VTune का उपयोग करके प्रोफाइलिंग करेंगे। यह सुनिश्चित करने के लिए प्रक्रिया को 100 बार चलाएं कि पर्याप्त डेटा है।



यहां मैंने प्रत्येक फ़ंक्शन के निष्पादन समय को ठीक करने के लिए, साथ ही बाधाओं को खोजने के लिए माइक्रोआर्किटेक्चर मोड चालू किया। हम देखते हैं कि कार्यों के एक सेट का उपयोग करके युक्तिकरण किया जाता है, जिनमें से प्रत्येक को निश्चित संख्या में चक्र की आवश्यकता होती है। कार्यों की सूची समय के अनुसार क्रमबद्ध होती है। एल्गोरिथ्म को समझने में आसान बनाने के लिए वे निष्पादन के क्रम में हैं:

• rescalePositions एक ही क्यूब में सभी कोने की स्थिति को सामान्य rescalePositions है, जो computeVertexIds का उपयोग करके परिमाणीकरण के लिए तैयार computeVertexIds
• computeVertexIds किसी दिए गए आकार के एक समान ग्रिड पर प्रत्येक शीर्ष के लिए एक 30-बिट परिमाणित पहचानकर्ता की गणना करता है, जहां प्रत्येक अक्ष को ग्रिड पर निर्धारित किया जाता है (ग्रिड आकार 10 बिट, इसलिए पहचानकर्ता 30 है)
• countTriangles त्रिकोण की अनुमानित संख्या की गणना करता है जो एक ग्रिड ग्रिड में सभी कोने के मिलन को देखते हुए इनोवेटर किसी दिए गए ग्रिड के आकार का निर्माण करेगा।
• fillVertexCells एक तालिका भरता है जो संबंधित कोशिकाओं के सभी कोने को fillVertexCells करता है; एक ही आईडी के साथ सभी कोने एक सेल के अनुरूप हैं
• fillCellQuadrics सेल ज्यामिति के बारे में समग्र जानकारी को प्रतिबिंबित करने के लिए प्रत्येक कोशिका के लिए Quadric संरचना ( Quadric सममित मैट्रिक्स) भरता है
• fillCellRemap प्रत्येक सेल के लिए वर्टेक्स इंडेक्स की गणना करता है, इस सेल में एक कोने का चयन करता है, और fillCellRemap कम करता है
• filterTriangles पहले बनाए गए वर्टेक्स-सेल-वर्टेक्स टेबल के अनुसार त्रिकोण का अंतिम सेट प्रदर्शित करता है; अनुभवहीन अनुवाद औसतन ~ 5% डुप्लिकेट त्रिकोण का उत्पादन कर सकता है, इसलिए फ़ंक्शन डुप्लिकेट को फ़िल्टर करता है।

computeVertexIds और countTriangles को कई बार निष्पादित किया जाता है: एल्गोरिथ्म countTriangles लिए मेष आकार को निर्धारित करता है, त्रिकोण की लक्ष्य संख्या (हमारे मामले में 6000) को प्राप्त करने के लिए एक त्वरित बाइनरी खोज का प्रदर्शन करता है और त्रिकोण की संख्या की गणना करता है जो प्रत्येक मेष आकार प्रत्येक पुनरावृत्ति पर उत्पन्न करेगा। अन्य कार्यों को एक बार लॉन्च किया जाता है। हमारी फ़ाइल में, पांच खोज पास 40 ³ का लक्ष्य जाल आकार देते हैं।

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 का उपयोग करके समानांतर किया जा सकता है। सबसे पहले, हम प्रत्येक वेक्टर को 4-व्यापी SIMD वेक्टर के रूप में Quadric हैं, और Quadric संरचना को 10 के बजाय 12 फ़्लोटिंग पॉइंट नंबरों में बदलते हैं, ताकि यह तीन SIMD रजिस्टरों में सटीक रूप से फिट हो जाए (आकार में वृद्धि प्रदर्शन को प्रभावित नहीं करती) और फ़ील्ड्स के क्रम को बदल दें ताकि गणना करें quadricFromPlane अधिक समान हो गया:

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

यहां, कुछ गणना, विशेष रूप से स्केलर उत्पाद, एसएसई के पुराने संस्करणों के अनुरूप नहीं हैं। सौभाग्य से, एक स्केलर उत्पाद के लिए एक निर्देश 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 } } } 

इस कोड में विशेष रूप से दिलचस्प कुछ भी नहीं है; हम बहुतायत से अनलग्ड लोड / स्टोर निर्देशों का उपयोग कर रहे हैं। हालाँकि, वेक्टर 3 के इनपुट को संरेखित किया जा सकता है, लेकिन बिना पढ़े लिखे के लिए कोई उल्लेखनीय दंड नहीं लगता है। कृपया ध्यान दें कि फ़ंक्शन के पहले छमाही में वैक्टर का उपयोग नहीं किया जाता है, जो अच्छा है - हमारे वैक्टर में तीन घटक होते हैं, और कुछ मामलों में केवल एक (क्षेत्रों / क्षेत्र / दूरी की गणना देखें), जबकि प्रोसेसर समानांतर में 4 ऑपरेशन करता है। किसी भी मामले में, आइए देखें कि समानांतरकरण ने कैसे मदद की है।



fillCellQuadrics की एक सौ शुरुआत अब 9.8 के बजाय 5.3 एस में चलती है, जो प्रत्येक ऑपरेशन पर लगभग 45 एमएस बचाता है - बुरा नहीं है, लेकिन बहुत प्रभावशाली नहीं है। कई निर्देशों में, हम चार घटकों के बजाय तीन का उपयोग करते हैं, साथ ही सटीक गुणा भी करते हैं, जो काफी महत्वपूर्ण देरी देता है। यदि आपने पहले SIMD के लिए निर्देश लिखे थे, तो आप जानते हैं कि स्केलर उत्पाद को सही तरीके से कैसे किया जाए।

ऐसा करने के लिए, आपको एक बार में चार वैक्टर करने चाहिए। एक SIMD रजिस्टर में एक पूर्ण वेक्टर को संग्रहीत करने के बजाय, हम तीन रजिस्टरों का उपयोग करते हैं - एक में हम x चार घटकों को स्टोर करते हैं, दूसरे में हम और तीसरे z में स्टोर करेंगे। यहां एक बार में काम के लिए चार वैक्टर की जरूरत होती है: इसका मतलब है कि हम एक साथ चार त्रिकोणों की प्रक्रिया करेंगे।

डायनेमिक इंडेक्सिंग के साथ हमारे पास कई एरेज़ हैं। आमतौर पर, यह x / y / z घटकों के तैयार सरणियों में डेटा स्थानांतरित करने में मदद करता है (या यों कहें कि आमतौर पर छोटे SIMD रजिस्टरों का उपयोग किया जाता है, उदाहरण के लिए, float x[8], y[8], z[8] , इनपुट में 8 में से प्रत्येक के लिए। डेटा: इसे AoSoA (सरणी संरचनाओं का सरणियाँ) कहा जाता है और कैश लोकेलिटी और SIMD रजिस्टरों में लोड करने में आसानी के बीच एक अच्छा संतुलन देता है), लेकिन यहां डायनामिक इंडेक्सिंग बहुत अच्छी तरह से काम नहीं करेगी, इसलिए हमेशा की तरह चार त्रिकोणीय डेटा लोड करें, और एक सुविधाजनक उपयोग करके वैक्टर को स्थानांतरित करें। मैक्रो _MM_TRANSPOSE

सैद्धांतिक रूप से, आपको अपने स्वयं के SIMD रजिस्टर में चार परिमित चतुष्कोणों के प्रत्येक घटक की गणना करने की आवश्यकता है (उदाहरण के लिए, हमारे पास __m128 Q_a00 साथ चार परिमित __m128 Q_a00 के चार घटक a00 )। इस मामले में, क्वाड्रिक्स पर ऑपरेशन 4-वाइड 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 निर्देशों का एक थोड़ा सा idiosyncratic सेट है। इसमें 8-वाइड फ्लोटिंग पॉइंट रजिस्टर हैं, और एक निर्देश 8 ऑपरेशन कर सकता है; लेकिन संक्षेप में, इस तरह के एक निर्देश रजिस्टर के दो हिस्सों पर निष्पादित दो SSE2 निर्देशों से भिन्न नहीं होता है (जहाँ तक मैं समझता हूं, AVX2 समर्थन वाले पहले प्रोसेसर ने प्रत्येक निर्देश को दो या अधिक माइक्रोफ़ाररेशन में डिकोड किया है, इसलिए प्रदर्शन लाभ निर्देश निकालने के चरण तक सीमित था)। उदाहरण के लिए, _mm_dp_ps दो SSE2 रजिस्टरों के बीच एक स्केलर उत्पाद करता है, और _mm256_dp_ps दो _mm256_dp_ps रजिस्टरों के दो हिस्सों के बीच दो स्केलर उत्पाद बनाता है, इसलिए यह प्रत्येक आधे के लिए 4-चौड़ा तक सीमित है।

इस वजह से, AVX2 कोड अक्सर सार्वभौमिक "8-वाइड SIMD" से भिन्न होता है, लेकिन यहां यह हमारे पक्ष में काम करता है: 4-वाइड वैक्टर का प्रत्यारोपण करके वैश्वीकरण में सुधार करने की कोशिश करने के बजाय, हम SIMD के पहले संस्करण में लौटते हैं और SSE2 के बजाय AVX2 निर्देशों का उपयोग करते हुए लूप को दोगुना करते हैं। / एसएसई ४। हमें अभी भी 4-वाइड वैक्टर को लोड करने और संग्रहीत करने की आवश्यकता है, लेकिन सामान्य तौर पर, हम _mm256 में __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); 

यहां आप प्राप्त Qx / Qz / Qz प्रत्येक 128-बिट आधे हिस्से को ले सकते हैं और उसी कोड को चला सकते हैं जिसे हम क्वाड्रिक्स में जोड़ते थे। इसके बजाय, हम मानते हैं कि अगर एक त्रिकोण में एक सेल में तीन कोने होते हैं ( 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 बार। लोड और स्टोर संचालन अभी भी 4-चौड़ा है, क्योंकि हमें डायनामिक इंडेक्सिंग के कारण असुविधाजनक मेमोरी लेआउट के साथ काम करने के लिए मजबूर किया जाता है, और गणना 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 के समूहों में) के लिए इकट्ठा रन कर सकते हैं और परिणामों की तुलना कर सकते हैं। इंटेल प्रोसेसर पर इकट्ठा ऐतिहासिक रूप से बहुत धीमा रहा है, लेकिन चलो इसे आज़माएं। सादगी के लिए, हम 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_TRANSPOSE8_LANE4_PS के बराबर है, जो मानक शीर्ष लेख में नहीं पाया जाता है, लेकिन आसानी से प्रदर्शित होता है। यह चार एवीएक्स 2 वैक्टर लेता है और एक दूसरे से स्वतंत्र रूप से दो 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 से आगे 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 जीबी वापस लिखना। 1.46 सेकंड में 28 जीबी की प्रोसेसिंग के लिए 19 जीबी / एस के बस बैंडविड्थ की आवश्यकता होती है। हम मेमोरी बैंडविड्थ को चलाकर चेक कर सकते हैं memcmp(block1, block2, 512 MB)। परिणाम 45 एमएस है, यानी, एक कोर पर लगभग 22 जीबी / एस (हालांकि एआईडीए 64 बेंचमार्क मेरे सिस्टम पर 31 जीबी / एस तक रीड गति दिखाता है, लेकिन यह कई कोर का उपयोग करता है)। वास्तव में, हम अधिकतम प्राप्त करने योग्य मेमोरी सीमा के करीब आ गए हैं, और प्रदर्शन में और वृद्धि के लिए इन कोने की करीब पैकिंग की आवश्यकता होगी ताकि वे 12 बाइट्स से कम पर कब्जा कर लें।

निष्कर्ष

हमने एक बहुत अच्छी तरह से अनुकूलित एल्गोरिथ्म लिया, जो 28 मिलियन त्रिकोण प्रति सेकंड की गति से बहुत बड़ी ग्रिड को सरल करता है, और SSE और AVX निर्देश सेट का उपयोग करके इसे लगभग दो बार, 50 मिलियन त्रिकोण प्रति सेकंड तक गति देता है। इस यात्रा के दौरान, हमें SIMD का उपयोग करने के लिए अलग-अलग तरीके सीखने पड़े: 3-वाइड वैक्टर को स्टोर करने के लिए रजिस्टर, SoA ट्रांसपोज़, AVX2 दो 3-वाइड वैक्टर स्टोर करने के लिए निर्देश, स्केलर निर्देशों की तुलना में डेटा लोडिंग को गति देने के लिए निर्देश इकट्ठा करना और आखिरकार। हमने सीधे स्ट्रीमिंग प्रोसेसिंग के लिए AVX2 लगाया।

SIMD अक्सर अनुकूलन के लिए सबसे अच्छा शुरुआती बिंदु नहीं होता है: प्लेटफ़ॉर्म-विशिष्ट निर्देशों का उपयोग किए बिना अल्गोरिथमिक ऑप्टिमाइज़ेशन और माइक्रोपोप्टिमाइज़ेशन के कई पुनरावृत्तियों के माध्यम से मेष तर्कसंगतता गई है। लेकिन कुछ बिंदु पर, ये संभावनाएं समाप्त हो जाती हैं, और यदि प्रदर्शन महत्वपूर्ण है, तो SIMD एक शानदार उपकरण है जिसका उपयोग यदि आवश्यक हो तो किया जा सकता है।

मुझे यकीन नहीं है कि इनमें से कौन सी अनुकूलन मुख्य शाखा में आएगी meshoptimizer: अंत में, यह केवल यह देखने के लिए एक प्रयोग है कि कोड एल्गोरिदम में मूलभूत परिवर्तन के बिना कितना ओवरक्लॉकिंग है। मुझे उम्मीद है कि लेख जानकारीपूर्ण निकला और आपको कोड के अनुकूलन के लिए विचार देगा। इस लेख के अंतिम स्रोत यहाँ हैं ; यह कार्य मेज़पोटीमाइज़र 99ab49 के संस्करण पर आधारित है, और थाई बुद्ध मॉडल स्केचफैब पर प्रकाशित हुआ है ।