🎯 ⚰️ ✋🏻 SIMD के साथ फ्लोट 4x4 मैट्रिक्स गुणन को गति दें 👨🏼‍⚕️ 🍣 🔤

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

उस समय, मैं विभिन्न प्रोसेसर पर समर्थन के साथ समस्याओं से बचना चाहता था। मैं अधिकतम उपलब्ध राशि पर अपने रेंडरर की जांच करने में सक्षम होना चाहता था। मेरे पास अभी भी पुराने AMD प्रोसेसर वाले दोस्त हैं, और उनकी छत SSE3 थी। इसलिए, उस समय मैंने खुद को अधिकतम SSE3 तक सीमित रखने का फैसला किया। तो एक वेक्टर गणितीय पुस्तकालय था, जो SSE3 के पहले एक दुर्लभ समावेश के साथ, SSE पर पूरी तरह से लागू होने से थोड़ा कम था। हालांकि, कुछ बिंदु पर मैंने सोचा कि मैं कितने महत्वपूर्ण प्रदर्शन कर सकता हूं जो कि कई महत्वपूर्ण वेक्टर गणित कार्यों के लिए प्रोसेसर से बाहर निचोड़ सकता है। ऐसा ही एक ऑपरेशन 4 मैट्रिसेस को 4 से गुणा करना है।

वास्तव में, मैंने इस व्यवसाय को मनोरंजन के लिए अधिक करने का फैसला किया। मैंने पहले से ही SSE पर अपने सॉफ्टवेयर प्रतिपादन के लिए मैट्रिक्स गुणा का उपयोग किया है और लिखा है और यह मेरे लिए पर्याप्त लगता है। लेकिन फिर मैंने यह देखने का फैसला किया कि मैं 2 फ्लोट 4x4 मैट्रिस को गुणा करने से सिद्धांत में कितने उपाय कर सकता हूं। मेरे वर्तमान SSE पर, ये 16 घड़ी चक्र हैं। सच है, IACA 3 के लिए हाल ही में संक्रमण 19 को दिखाना शुरू किया, क्योंकि मैंने 0 * के बजाय कुछ निर्देशों के लिए 1 * लिखना शुरू किया। स्पष्ट रूप से पहले यह केवल विश्लेषक में एक दोष था।

उपयोग की जाने वाली उपयोगिताओं के बारे में संक्षेप में

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

असेंबली के लिए मैं कंसोल से MSVS 2017 समुदाय का उपयोग करता हूं। मैं संस्करण में आंतरिक के साथ कोड लिखता हूं। आप एक बार लिखते हैं, और आमतौर पर यह तुरंत विभिन्न प्लेटफार्मों पर काम करता है। इसके अलावा, x64 VC ++ कंपाइलर इनलाइन असेंबलर का समर्थन नहीं करता है, लेकिन मैं चाहता हूं कि यह x64 के तहत भी काम करे।

चूंकि यह लेख पहले से ही SIMD प्रोग्रामिंग में शुरुआती स्तर से परे है, मैं सुंदर चित्रों के रजिस्टरों, निर्देशों, ड्रा (या दाढ़ी) का वर्णन नहीं करूंगा और SIMD निर्देशों का उपयोग करके प्रोग्रामिंग सीखने की कोशिश करूंगा। इंटेल वेबसाइट उत्कृष्ट, स्पष्ट और विस्तृत प्रलेखन से भरी है।

मैं सब कुछ आसान करना चाहता था। ... लेकिन यह हमेशा की तरह निकला

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

लेकिन विशेष रूप से मेट्रिसेस पर वापस। हमने कॉलम में भंडारण का निर्धारण किया। हालाँकि, यह और जटिल हो सकता है। SIMD रजिस्टरों में वैक्टर और मैट्रिक्स पंक्तियों के वरिष्ठ तत्वों को संग्रहीत करना मेरे लिए अधिक सुविधाजनक है ताकि x उच्चतम फ्लोट (अनुक्रमणिका 3) में हो, और माइनर (अनुक्रमणिका 0) में w । यहाँ, जाहिर है, हमें फिर से पीछे हटना होगा क्यों।

तथ्य यह है कि एक वेक्टर में एक सॉफ्टवेयर रेंडरर में, आपको डब्ल्यू घटक को अधिक बार हेरफेर करना होगा ( 1 / z वहां संग्रहीत किया जाता है), और ऑपरेशन के _ss संस्करण ( एक्सएमएम रजिस्टर के निचले फ्लोट में घटक के साथ विशेष रूप से संचालन) के माध्यम से ऐसा करना बहुत आसान है, एक्स, वाई को छूने के बिना । z इसलिए, 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 ने जारी किया: x86 विधानसभा के लिए 3.65 घड़ी चक्र और x64 विधानसभा के लिए 2.97 घड़ी । यह मत पूछो कि भिन्नात्मक संख्याएँ क्यों हैं। मुझे नहीं पता 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 { //           union { struct { float _30, _20, _10, _00, _31, _21, _11, _01, _32, _22, _12, _02, _33, _23, _13, _03; }; __m128 r[4]; __m256 s[2]; vec4 v[4]; }; //    mtx4() {} mtx4( float i00, float i01, float i02, float i03, float i10, float i11, float i12, float i13, float i20, float i21, float i22, float i23, float i30, float i31, float i32, float i33) : _00(i00), _01(i01), _02(i02), _03(i03) , _10(i10), _11(i11), _12(i12), _13(i13) , _20(i20), _21(i21), _22(i22), _23(i23) , _30(i30), _31(i31), _32(i32), _33(i33) {} //      operator __m128 const* () const { return r; } operator __m128* () { return r; } //   bool operator == (mtx4 const& m) const { return v[0]==mv[0] && v[1]==mv[1] && v[2]==mv[2] && v[3]==mv[3]; } //  static mtx4 identity() { return mtx4( 1.f, 0.f, 0.f, 0.f, 0.f, 1.f, 0.f, 0.f, 0.f, 0.f, 1.f, 0.f, 0.f, 0.f, 0.f, 1.f); } static mtx4 zero() { return mtx4( 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f, 0.f); } };

परीक्षण के लिए संदर्भ समारोह

चूंकि मैट्रिक्स में तत्वों का स्वीकृत क्रम समझ को बहुत जटिल करता है, इसलिए हम संदर्भ स्पष्ट फ़ंक्शन से भी परेशान नहीं होंगे, जो भविष्य के कार्यान्वयन में दिखाएगा कि सब कुछ सही ढंग से काम करता है। हम इसके साथ बाद के परिणामों की तुलना करेंगे।

इसे बनाने के लिए, बस चक्र को लें और उसका विस्तार करें

 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 एल्गोरिथम

क्‍यों क्लासिक? क्योंकि यह लंबे समय से एमएसवीएस के हिस्से के रूप में एफवीईसी के कार्यान्वयन में है। शुरू करने के लिए, हम लिखेंगे कि हम एसएसई रजिस्टरों में मैट्रिक्स तत्व कैसे प्रस्तुत करते हैं। यह पहले से ही यहाँ सरल दिखता है। बस एक ट्रांसपोंड मैट्रिक्स।

 //     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 । हमारे पास यह कॉलम है, लेकिन हमें एक पंक्ति की आवश्यकता है। हां, और मी , एन गणना के लिए असुविधाजनक दिखते हैं। इसलिए, हम एल्गोरिथम की पंक्तियों को पुनर्व्यवस्थित करते हैं ताकि परिणाम आर पंक्ति-वार हो।

 //  ,     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}, m [1] = {01,11,21,31}, m [2] = {02,12,22,32}, m [3] = {} 03,13,23,33,

जो मैट्रिक्स n के एक ही तत्व द्वारा गुणा किए जाते हैं। उदाहरण के लिए, पहले समूह के लिए यह है: n._00, n._10, n._20, n._30 । और एल्गोरिथ्म की पंक्तियों के प्रत्येक समूह के लिए मैट्रिक्स के तत्व फिर से मैट्रिक्स की एक पंक्ति में झूठ बोलते हैं।

तब सब कुछ सरल है: हम केवल मैट्रिक्स मीटर की पंक्तियों को सूचकांक द्वारा लेते हैं, लेकिन तत्वों n के लिए , हम इसकी पंक्ति लेते हैं और इसे रजिस्टर के सभी 4 तत्वों के माध्यम से फेरबदल करते हैं, रजिस्टर में मैट्रिक्स मीटर की पंक्ति द्वारा गुणा करने के लिए। उदाहरण के लिए, तत्व 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;

बुनियादी एसएसई कार्यान्वयन

 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))))); }

अब हम एल्गोरिथ्म के तत्वों को संबंधित फेरबदल में बदलते हैं , _mm_mul_ps द्वारा गुणा, _mm_add_ps द्वारा योग, और आप कर रहे हैं । यह काम करता है। कोड, हालांकि, एल्गोरिथ्म की तुलना में बहुत खराब लग रहा था। इस कोड के लिए, IACA ने जारी किया: x86 - 18.89, x64 - 16 चक्र । यह पिछले वाले की तुलना में 4 गुना तेज है। SSE रजिस्टर में 4th फ्लोट। लगभग रैखिक संबंध।

एसएसई कार्यान्वयन को सजाने

फिर भी, कोड में यह भयानक लग रहा है। हम थोड़ा सा संश्लिष्ट चीनी लिखकर इसे सुधारने का प्रयास करेंगे।

संचालक और सुधारक

 //        ( -    namespace) __m128 operator + (__m128 const a, __m128 const b) { return _mm_add_ps(a, b); } __m128 operator - (__m128 const a, __m128 const b) { return _mm_sub_ps(a, b); } __m128 operator * (__m128 const a, __m128 const b) { return _mm_mul_ps(a, b); } __m128 operator / (__m128 const a, __m128 const b) { return _mm_div_ps(a, b); } //_mm_shuffle_ps(u, v, _MM_SHUFFLE(3,2,1,0))   shuf<3,2,1,0>(u, v) template <int a, int b, int c, int d> __m128 shuf(__m128 const u, __m128 const v) { return _mm_shuffle_ps(u, v, _MM_SHUFFLE(a, b, c, d)); } template <int a, int b, int c, int d> __m128 shuf(__m128 const v) { return _mm_shuffle_ps(v, v, _MM_SHUFFLE(a, b, c, d)); } //    template <int i> __m128 shuf(__m128 const u, __m128 const v) { return _mm_shuffle_ps(u, v, _MM_SHUFFLE(i, i, i, i)); } template <int i> __m128 shuf(__m128 const v) { return _mm_shuffle_ps(v, v, _MM_SHUFFLE(i, i, i, i)); } //  float       , //    ,    template <int a, int b, int c, int d> __m128 shufd(__m128 const v) { return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(v), _MM_SHUFFLE(a, b, c, d))); } template <int i> __m128 shufd(__m128 const v) { return _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(v), _MM_SHUFFLE(i, i, i, i))); }

कंपाइलर इन कार्यों को पूरी तरह से इनलाइन कर सकता है (हालांकि कभी-कभी __फोर्सिनलाइन के बिना किसी भी तरह से)।

इसलिए कोड बदल रहा है ...

 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 (ए, बी, सी) = ए * बी + सी ।

कई-कई कार्यान्वयन

 __m128 mad(__m128 const a, __m128 const b, __m128 const c) { return _mm_add_ps(_mm_mul_ps(a, b), c); }

यह क्यों आवश्यक है? भविष्य के अनुकूलन के लिए सबसे पहले। उदाहरण के लिए, आप बस उसी मैक्रोज़ के माध्यम से तैयार कोड में 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 विकल्प

अगला, हम अनुकूलन के अगले चरण पर जाते हैं। 4 वें फ्लोट को एसएसई रजिस्टर में शामिल किया गया है, और एवीएक्स में यह पहले से ही 8. है। अर्थात्, प्रदर्शन किए गए संचालन की संख्या को कम करने के लिए एक सैद्धांतिक मौका है, और उत्पादकता में वृद्धि अगर आधे से नहीं, तो कम से कम 1.5 गुना। लेकिन कुछ मुझे बताता है कि एवीएक्स के संक्रमण के साथ सब कुछ इतना सरल नहीं होगा। क्या हम डबल रजिस्टर से आवश्यक डेटा प्राप्त कर सकते हैं?

आइए इसे जानने की कोशिश करते हैं। फिर, हम अपने गुणन एल्गोरिथ्म का उपयोग ऊपर लिखते हैं। आप ऐसा नहीं कर सकते, लेकिन जब सब कुछ पास हो तो कोड से निपटना अधिक सुविधाजनक होगा और आपको आधे पृष्ठ तक स्क्रॉल नहीं करना पड़ेगा।

 //    : 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 संस्करण में हमारे एल्गोरिथ्म को स्तंभों के लिए सामान्यीकृत किया गया था, तो अब हमें इसे पंक्तियों के लिए सामान्यीकृत करने की आवश्यकता है। तो एसएसई विकल्प के मामले में काम करने के लिए काम नहीं करेगा।

यदि हम रजिस्टरों ymm में मैट्रिक्स m पर विचार करते हैं, तो हम क्रमशः ymm = {m0: m1} और ymm = {m2: m3} प्राप्त करते हैं। पहले, हमारे पास रजिस्टर में केवल मैट्रिक्स कॉलम थे, और अब कॉलम और पंक्तियाँ हैं।

यदि आप पहले की तरह कार्य करने का प्रयास करते हैं, तो आपको ymm = {m0: m1} को रजिस्टर ymm = {n00, n00, n00, n00} से गुणा करना होगा : {n10, n10, n10, n10} । चूंकि n00 और n01 मैट्रिक्स n की एक ही पंक्ति में हैं, एवीएक्स निर्देशों के उपलब्ध सेट को देखते हुए, उन्हें यम द्वारा बिखेरना महंगा होगा। शफल और परमिट दोनों यम रजिस्टरों के अंदर एक फ्लोट (उच्च और निम्न एक्सएमएम ) के दो चार में से प्रत्येक के लिए अलग से काम करते हैं।

यदि हम मैट्रिक्स n से ymm लेते हैं, तो हम दोनों तत्वों n00 और n10 को ymm रजिस्टर के अंदर अधिकतम 2 xmm में प्राप्त करते हैं। {n00, n10, n20, n30}: {n01, n11, n21, n31} । आमतौर पर, मौजूदा निर्देशों के लिए सूचकांक 0 से 3 तक होता है। और यह पता चलता है कि केवल एक xmm रजिस्टर के अंदर दो ymm रजिस्टर में से तैरता है । सस्ते में पुराने xmm से n10 को सस्ते में स्थानांतरित करना संभव नहीं है। और फिर इस फोकस को कई बार दोहराया जाना चाहिए। हम इस तरह के नुकसान से नहीं जूझ सकते। कुछ और के साथ आना जरूरी है।

हम स्तंभों को सामान्य करते थे, लेकिन अब पंक्तियों को। इसलिए, हम थोड़ा अलग तरीके से जाने की कोशिश करेंगे। हमें {r0: r1} में परिणाम प्राप्त करने की आवश्यकता है। इसका मतलब यह है कि एल्गोरिथ्म को एल्गोरिथ्म की अलग-अलग लाइनों में नहीं, बल्कि एक साथ दो में सुधार किया जाना चाहिए। और यहाँ, क्या फेरबदल और परमिट के काम में एक शून्य था जो हमारे लिए एक प्लस होगा। हम देखते हैं कि जब हम मैट्रिक्स एन पर विचार करते हैं तो हमारे पास यम रजिस्टरों में क्या होगा।

 n0n1 = {00, 10, 20, 30} : {01, 11, 21, 31} n2n3 = {02, 12, 22, 32} : {03, 13, 23, 33}

हाँ, हम देखते हैं कि ymm register के विभिन्न xmm भागों में हमारे पास 00 और 01 तत्व हैं। उन्हें {_00, _00, _00, _00} में परमिट कमांड के माध्यम से रजिस्टर करके गुणा किया जा सकता है : {_ 01, _01, _01, _01} , दोनों xmm भागों के लिए केवल एक इंडेक्स 3 को दर्शाता है । यह वही है जो हमें चाहिए। दरअसल, गुणांक का उपयोग विभिन्न लाइनों में भी किया जाता है। गुणन के लिए संबंधित यम रजिस्टर में अब केवल {m0: m0} , यानी मैट्रिक्स m की डुप्लीकेटेड पहली पंक्ति को रखना आवश्यक होगा।

इसलिए, हम एल्गोरिथ्म को अधिक विस्तार से चित्रित करते हैं। हम ymm रजिस्टरों में मैट्रिक्स m की दोहरी पंक्तियों को पढ़ते हैं:

 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);//3,3,3,3 __m256 y2 = _mm256_permute_ps(n0n1, 0xAA);//2,2,2,2 __m256 y3 = _mm256_permute_ps(n0n1, 0x55);//1,1,1,1 __m256 y4 = _mm256_permute_ps(n0n1, 0x00);//0,0,0,0 y1 = _mm256_mul_ps(y1, mm0); y2 = _mm256_mul_ps(y2, mm1); y3 = _mm256_mul_ps(y3, mm2); y4 = _mm256_mul_ps(y4, mm3); y1 = _mm256_add_ps(y1, y2); y3 = _mm256_add_ps(y3, y4); y1 = _mm256_add_ps(y1, y3); __m256 n2n3 = _mm256_load_ps(&n[2].m128_f32[0]); __m256 y5 = _mm256_permute_ps(n2n3, 0xFF); __m256 y6 = _mm256_permute_ps(n2n3, 0xAA); __m256 y7 = _mm256_permute_ps(n2n3, 0x55); __m256 y8 = _mm256_permute_ps(n2n3, 0x00); y5 = _mm256_mul_ps(y5, mm0); y6 = _mm256_mul_ps(y6, mm1); y7 = _mm256_mul_ps(y7, mm2); y8 = _mm256_mul_ps(y8, mm3); y5 = _mm256_add_ps(y5, y6); y7 = _mm256_add_ps(y7, y8); y5 = _mm256_add_ps(y5, y7); _mm256_stream_ps(&r[0].m128_f32[0], y1); _mm256_stream_ps(&r[2].m128_f32[0], y5); }

यहाँ IACA से रोचक संख्याएँ हैं: x86 - 12.53, x64 - 12 । हालांकि, निश्चित रूप से, मैं बेहतर चाहता था। कुछ याद आया।

एवीएक्स ऑप्टिमाइज़ेशन प्लस सिंथैटिक शुगर

ऐसा लगता है कि उपरोक्त कोड में, AVX को इसकी पूरी क्षमता के लिए इस्तेमाल नहीं किया गया था। हम पाते हैं कि यम रजिस्टर में दो समान रेखाएँ सेट करने के बजाय, हम प्रसारण का उपयोग कर सकते हैं, जो दो समान मिमी मानों के साथ यम रजिस्टर भर सकते हैं। इसके अलावा, रास्ते में, 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 , जो बहुत बेहतर दिखता है, लेकिन फिर भी 2 बार नहीं।

AVX + FMA विकल्प (अंतिम)

चलिए एक और प्रयास करते हैं। अब एएमएक्स के बाद इसे प्रोसेसर में जोड़ा गया था, फिर से एफएमए निर्देश को फिर से याद करना तर्कसंगत होगा। बस अलग-अलग mul + को एक ऑपरेशन के लिए जोड़ें । यद्यपि हम अभी भी कंपाइलर को अनुकूलन के लिए अधिक अवसर देने के लिए गुणन के निर्देश का उपयोग करते हैं, और प्रोसेसर गुणन के समानांतर निष्पादन के लिए। आमतौर पर मैं असेंबलर में उत्पन्न कोड को यह सुनिश्चित करने के लिए देखता हूं कि कौन सा विकल्प बेहतर है।

इस मामले में, हमें एक * बी + सी * डी + ई * एफ + जी * एच की गणना करने की आवश्यकता है। आप यह माथा कर सकते हैं: fma (a, b, fma (c, d, fma (e, f, g * h))) । लेकिन, जैसा कि हम देखते हैं, पिछले एक को पूरा किए बिना यहां एक ऑपरेशन करना असंभव है। और इसका मतलब है कि हम युग्मित गुणा करने की क्षमता का उपयोग नहीं कर पाएंगे, क्योंकि SIMD पाइपलाइन हमें करने की अनुमति देती है। यदि हम गणनाओं को थोड़ा बदल देते हैं तो fma (a, b, c * d) + fma (e, f, g * h) , हम देखेंगे कि हम गणना को समानांतर कर सकते हैं। पहले दो स्वतंत्र गुणन करते हैं, और फिर दो स्वतंत्र 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 । अब यह बहुत अच्छा है कोई शायद कहेगा कि क्या बेहतर किया जा सकता है, लेकिन मैं नहीं जानता कि कैसे।

मानक

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

सामग्री की तालिका

समारोह: समारोह का नाम। एस पर समाप्त - एक स्ट्रीमिंग के साथ कार्य करता है, अलग-अलग सामान्य Mov (स्ट्रीमिंग के बिना)। स्पष्टता के लिए जोड़ा गया, क्योंकि यह पर्याप्त महत्वपूर्ण है।
IACA चक्र: IACA द्वारा गणना की गई प्रति फ़ंक्शन टिक की संख्या
मापा चक्र: मापी गई संख्या की माप (कम अधिक है)
IACA स्पीडअप: एक शून्य रेखा में उपायों की संख्या / एक पंक्ति में उपायों की संख्या
मापा गति: शून्य लाइन में उपायों की संख्या / लाइन में उपायों की संख्या (अधिक बेहतर)

लूप_म के लिए, लेख के टिक्स को 64 से गुणा किया गया था। अर्थात्, यह एक बहुत ही अनुमानित मूल्य है। वास्तव में, यह इस तरह से निकला।

i3-3770:

समारोह	IACA चक्र	मापा चक्र	IACA स्पीडअप	मापा गया स्पीडअप
unroll_m	70.00	50.75	1.00	1.00
loop_m	233.60	119.21	0.30	0.43
sse_v1m	18.89	27.51	3.70	1.84
sse_v2m	19.00	27.61	3.68	1.84
sse_v3m	18.89	27.22	3.70	1.86
sse_v4s	18.89	27.18	3.70	1.87
avx_v1m	13.00	19.21	5.38	2.64
avx_v1s	13.00	20.03	5.38	2.53
avx_v2m	10.00	12.91	6.99	3.93
avx_v2s	10.00	17.34	6.99	2.93

Function	IACA cycles	Measured cycles	IACA speedup	Measured speedup
unroll_m	70	68.60	1.00	1.00
loop_m	233.60	119.37	0.30	0.57
sse_v1m	18.89	21.98	3.70	3.12
sse_v2m	19.00	21.09	3.68	3.25
sse_v3m	18.89	22.19	3.70	3.09
sse_v4s	18.89	22.39	3.70	3.06
avx_v1m	13.00	9.61	5.38	7.13
avx_v1s	13.00	16.90	5.38	4.06
avx_v2m	10.00	9.20	6.99	7.45
avx_v2s	10.00	14.64	6.99	4.68

i7-8700K:

Function	IACA cycles	Measured cycles	IACA speedup	Measured speedup
unroll_m	69.95	40.25	1.00	1.00
loop_m	233.60	79.49	0.30	0.51
sse_v1m	18.89	19.31	3.70	2.09
sse_v2m	19.00	19.98	3.68	2.01
sse_v3m	18.89	19.69	3.70	2.04
sse_v4s	18.89	19.67	3.70	2.05
avx_v1m	13.00	14.22	5.38	2.83
avx_v1s	13.00	14.13	5.38	2.85
avx_v2m	10.00	11.73	6.99	3.43
avx_v2s	10.00	11.81	6.99	3.41
AVX+FMAm	9.21	10.38	7.60	3.88
AVX+FMAs	9.21	10.32	7.60	3.90

Function	IACA cycles	Measured cycles	IACA speedup	Measured speedup
unroll_m	69.95	57.11	1.00	1.00
loop_m	233.60	75.73	0.30	0.75
sse_v1m	18.89	15.83	3.70	3.61
sse_v2m	19.00	17.22	3.68	3.32
sse_v3m	18.89	15.92	3.70	3.59
sse_v4s	18.89	16.18	3.70	3.53
avx_v1m	13.00	7.03	5.38	8.12
avx_v1s	13.00	12.98	5.38	4.40
avx_v2m	10.00	5.40	6.99	10.57
avx_v2s	10.00	11.39	6.99	5.01
AVX+FMAm	9.21	9.73	7.60	5.87
AVX+FMAs	9.21	9.81	7.60	5.82

स्रोत में टेस्ट कोड। यदि उनमें सुधार करने के बारे में उचित सुझाव हैं, तो टिप्पणियों में लिखें।

कल्पना के दायरे से अलग

वास्तव में, यह कल्पना के दायरे से है क्योंकि अगर मैंने AVX512 के लिए प्रोसेसर देखा, तो शायद तस्वीरों में। हालांकि, मैंने एल्गोरिथ्म को लागू करने की कोशिश की। यहां मैं कुछ भी नहीं समझाऊंगा, एवीएक्स + एफएमए के साथ एक पूर्ण सादृश्य। एल्गोरिथ्म एक ही है, केवल कम संचालन।

जैसा कि वे कहते हैं, मैं इसे यहाँ छोड़ दूँगा

 __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 (SKX वास्तुकला के साथ IACA)। इस तथ्य के बावजूद कि एल्गोरिथ्म में 64 गुणा और 48 जोड़ हैं।

लेख से PS कोड

https://github.com/truthfinder/mul_m4_m4

यह एक लेख लिखने का मेरा पहला अनुभव है। आपकी टिप्पणियों के लिए आप सभी का धन्यवाद। वे कोड और लेख को बेहतर बनाने में मदद करते हैं।

SIMD के साथ फ्लोट 4x4 मैट्रिक्स गुणन को गति दें