Mempercepat penggandaan matriks 4x4 dengan SIMD

Bertahun-tahun telah berlalu sejak saya berkenalan dengan instruksi MMX, SSE, dan kemudian AVX pada prosesor Intel. Pada suatu waktu, mereka tampak seperti semacam sihir dengan latar belakang assembler x86, yang telah lama menjadi sesuatu yang biasa. Mereka sangat mengaitkan saya sehingga beberapa tahun yang lalu saya memiliki ide untuk menulis perender perangkat lunak saya sendiri untuk satu permainan terkenal. Hal yang menjanjikan instruksi ini menjanjikan saya hal ini. Pada titik tertentu, saya bahkan berpikir untuk menulisnya. Tetapi menulis teks ternyata jauh lebih rumit daripada kode.

Saat itu, saya ingin menghindari masalah dengan dukungan pada prosesor yang berbeda. Saya ingin dapat memeriksa penyaji saya pada jumlah maksimum yang tersedia. Saya masih punya teman dengan prosesor AMD lama, dan langit-langit mereka adalah SSE3. Oleh karena itu, pada saat itu saya memutuskan untuk membatasi diri hingga maksimum SSE3. Jadi ada perpustakaan matematika vektor, sedikit kurang sepenuhnya diimplementasikan pada SSE, dengan inklusi langka sebelum SSE3. Namun, pada titik tertentu saya bertanya-tanya kinerja maksimum apa yang bisa saya peras dari prosesor untuk sejumlah operasi matematika vektor kritis. Salah satu operasi tersebut adalah dengan mengalikan float 4 matriks dengan 4.


Sebenarnya, saya memutuskan untuk melakukan bisnis ini lebih demi hiburan. Saya sudah menulis dan telah menggunakan perkalian matriks untuk perenderan perangkat lunak saya di SSE dan tampaknya sudah cukup bagi saya. Tetapi kemudian saya memutuskan untuk melihat berapa banyak ukuran yang bisa saya peras pada prinsipnya dengan mengalikan 2 matriks float4x4. Pada SSE saya saat ini, ini adalah 16 siklus clock. Benar, transisi terakhir ke IACA 3 mulai menunjukkan 19, karena saya mulai menulis 1 * untuk beberapa instruksi daripada 0 *. Rupanya sebelumnya itu hanya cacat pada alat analisa.

Secara singkat tentang utilitas yang digunakan


Untuk analisis kode, saya menggunakan utilitas Intel Architecture Code Analyzer yang terkenal. Untuk analisis saya menggunakan arsitektur Haswell (HSW), sebagai minimum dengan dukungan untuk AVX2. Untuk menulis kode juga sangat mudah digunakan: Intel Intrinsics Guide dan manual optimasi Intel .

Untuk perakitan saya menggunakan Komunitas MSVS 2017 dari konsol. Saya menulis kode dalam versi dengan intrinsik. Anda menulis sekali, dan biasanya itu langsung berfungsi pada platform yang berbeda. Selain itu, kompiler x64 VC ++ tidak mendukung inline assembler, tapi saya ingin itu bekerja di bawah x64 juga.

Karena artikel ini sudah agak melampaui level pemula dalam pemrograman SIMD, saya tidak akan menjelaskan register, instruksi, menggambar (atau mencukur) gambar yang indah dan mencoba untuk belajar pemrograman menggunakan instruksi SIMD. Situs web Intel penuh dengan dokumentasi yang luar biasa, jelas, dan terperinci.

Saya ingin membuat segalanya lebih mudah ... Tapi ternyata seperti biasa


Di sinilah saat dimulainya, yang menyulitkan banyak implementasi dan artikel. Karena itu, saya akan sedikit membahasnya. Tidak menarik bagi saya untuk menulis perkalian matriks dengan tata letak baris elemen standar. Siapa yang membutuhkannya, dan mereka belajar di universitas atau sendiri. Tujuan kami adalah produktivitas. Pertama, saya beralih ke tata letak kolom sejak lama. Perender perangkat lunak saya didasarkan pada API OpenGL dan oleh karena itu, untuk menghindari transposisi yang tidak perlu, saya mulai menyimpan elemen dalam kolom. Ini juga penting karena perkalian matriks tidak begitu kritis. Matikan 2-5-10 dikalikan dengan baik. Dan itu dia. Dan kemudian kita mengalikan matriks jadi dengan ribuan atau jutaan titik. Dan operasi ini jauh lebih kritis. Anda tentu saja dapat mengubah posisi setiap waktu. Tapi kenapa, kalau ini bisa dihindari.

Tetapi kembali ke matriks secara eksklusif. Kami menentukan penyimpanan dalam kolom. Namun, ini bisa menjadi lebih rumit. Lebih nyaman bagi saya untuk menyimpan elemen senior dari vektor dan baris matriks dalam register SIMD sehingga x berada di float tertinggi (indeks 3), dan w di bawah umur (indeks 0). Di sini, tampaknya, kita harus mundur lagi tentang alasannya.

Faktanya adalah bahwa dalam perender perangkat lunak dalam vektor, Anda harus memanipulasi komponen w lebih sering ( 1 / z disimpan di sana), dan sangat mudah untuk melakukan ini melalui versi _ss operasi (operasi secara eksklusif dengan komponen dalam float yang lebih rendah dari register xmm ), tanpa menyentuh x, y, z . Oleh karena itu, dalam register SSE, vektor disimpan dalam urutan yang dapat dimengerti x, y, z, w , dan dalam memori dalam w, z, y, x terbalik.

Lebih lanjut, semua opsi multiplikasi juga diterapkan oleh fungsi individual. Ini dilakukan karena saya menggunakannya untuk menggantikan opsi yang diinginkan tergantung pada jenis instruksi yang didukung. Dijelaskan dengan baik di sini.

Kami menerapkan fungsionalitas dasar


Perkalian dengan loop, baris dipesan


Opsi untuk tata letak elemen
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]; } } } 


Semuanya sederhana dan jelas di sini. Untuk setiap elemen kami melakukan 4 perkalian dan 3 tambahan. Secara total, ini adalah 64 perkalian dan 48 tambahan. Dan ini tanpa memperhitungkan pembacaan elemen rekaman.

Singkatnya, semuanya menyedihkan. Untuk opsi ini, untuk siklus internal, IACA mengeluarkan: 3.65 siklus clock untuk rakitan x86 dan 2.97 jam untuk rakitan x64 . Jangan tanya kenapa angka pecahan. Saya tidak tahu. IACA 2.1 tidak menderita dari ini. Bagaimanapun, angka-angka ini harus dikalikan dengan sekitar 4 * 4 * 4 = 64. Bahkan jika Anda mengambil x64, hasilnya adalah sekitar 192 tindakan. Jelas bahwa ini adalah perkiraan kasar. Saya tidak melihat titik mengevaluasi kinerja lebih akurat untuk opsi ini.

Implementasi looping, kolom dipesan


matriks transposited, mengatur ulang indeks baris dan kolom
 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]; } } } 


Penggandaan siklus, penyimpanan berorientasi SIMD


penyimpanan garis dalam urutan terbalik dalam memori ditambahkan
 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]; } } } 


Implementasi ini agak menyederhanakan pemahaman tentang apa yang terjadi di dalam, tetapi jelas tidak cukup.

Kelas pembantu


Untuk kenyamanan memahami dan menulis referensi dan kode debug, lebih mudah untuk mengimplementasikan beberapa kelas tambahan. Tidak lebih, semuanya hanya untuk pengertian. Saya perhatikan bahwa implementasi penuh vektor dan matriks kelas adalah pertanyaan sulit yang terpisah, dan tidak termasuk dalam topik artikel ini.

Kelas vektor dan matriks
 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); } }; 


Fungsi referensi untuk tes


Karena urutan elemen yang diterima dalam matriks mempersulit banyak pemahaman, kami juga tidak akan terganggu oleh fungsi referensi yang dapat dimengerti , yang akan menunjukkan dalam implementasi di masa depan bahwa semuanya bekerja dengan benar. Kami akan membandingkan hasil selanjutnya dengan itu.

Untuk membuatnya, cukup ambil dan perluas siklusnya
 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; } 


Algoritma klasik dilukis dengan jelas di sini, sulit untuk membuat kesalahan (tetapi Anda dapat :-)). Di atasnya, IACA mengeluarkan: ukuran x86 - 69,95, ukuran x64 - 64 . Berikut adalah sekitar 64 siklus dan kami akan melihat percepatan operasi ini di masa depan.

Implementasi SSE


Algoritma SSE Klasik


Mengapa klasik? Karena sudah lama dalam implementasi FVec sebagai bagian dari MSVS. Untuk mulai dengan, kami akan menulis bagaimana kami menyajikan elemen matriks dalam register SSE. Ini sudah terlihat lebih sederhana di sini. Hanya matriks transpos.

 //     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] 

Kami mengambil kode membuka gulungan varian di atas. Entah bagaimana dia tidak ramah untuk SSE. Kelompok baris pertama terdiri dari hasil untuk kolom matriks yang dihasilkan: r._00, r._01, r._02, r._03 . Kami memiliki kolom ini, tetapi kami membutuhkan satu baris. Ya, dan m , n terlihat tidak nyaman untuk perhitungan. Oleh karena itu, kami mengatur ulang garis-garis algoritma sehingga hasilnya r adalah baris-bijaksana.

 //  ,     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; 

Tapi ini sudah jauh lebih baik. Sebenarnya, apa yang kita lihat? Menurut kolom algoritma di setiap kelompok, kami memiliki baris matriks m yang terlibat:
 m [0] = {00,10,20,30}, m [1] = {01,11,21,31}, m [2] = {02,12,22,32}, m [3] = {03,13,23,33},
yang dikalikan dengan elemen yang sama dari matriks n . Misalnya, untuk grup pertama adalah: n._00, n._10, n._20, n._30 . Dan elemen-elemen dari matriks n untuk setiap kelompok baris dari algoritma lagi terletak pada satu baris dari matriks.

Kemudian semuanya sederhana: kita cukup mengambil baris matriks m berdasarkan indeks, tetapi untuk elemen n , kita mengambil barisnya dan mengocoknya melalui keempat elemen register melalui instruksi shuffle , untuk dikalikan dengan baris matriks m dalam register. Misalnya, untuk elemen n._00 (ingat bahwa perubahannya dalam register memiliki indeks 3), itu akan menjadi:
  _mm_shuffle_ps (n [0], n [0], _MM_SHUFFLE (3,3,3,3)) 

Dalam bentuk yang disederhanakan, algoritmenya terlihat seperti ini:

 //   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; 

Implementasi SSE dasar
 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))))); } 


Sekarang kita mengubah elemen n dalam algoritma ke shuffle yang sesuai, perkalian dengan _mm_mul_ps , jumlah dengan _mm_add_ps , dan Anda selesai. Itu bekerja. Akan tetapi, kode tersebut terlihat jauh lebih buruk daripada algoritma yang terlihat. Untuk kode ini, IACA mengeluarkan: x86 - 18.89, x64 - 16 siklus . Ini 4 kali lebih cepat dari yang sebelumnya. Dalam daftar SSE float ke-4. Hubungannya hampir linier.

Hiasi implementasi SSE


Namun, dalam kode itu terlihat mengerikan. Kami akan mencoba memperbaikinya dengan menulis sedikit gula sintaksis.

Operator dan Peningkatan
 //        ( -    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))); } 


Kompiler dapat dengan sempurna menyejajarkan fungsi-fungsi ini (meskipun kadang-kadang tanpa __forceinline dengan cara apa pun).

Jadi kodenya berubah ...
 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]); } 


Dan itu sudah jauh lebih baik dan lebih mudah dibaca. Untuk ini, IACA menghasilkan sekitar hasil yang diharapkan: x86 - 19 (dan mengapa tidak fraksional?), X64 - 16 . Sebenarnya, kinerjanya tidak berubah, tetapi kodenya jauh lebih indah dan dapat dimengerti.

Kontribusi kecil untuk optimasi di masa depan


Mari kita perkenalkan satu peningkatan lagi di tingkat fungsi yang baru-baru ini muncul di versi besi. Operasi multiple-add (fma) . fma (a, b, c) = a * b + c .

Implementasi multi-tambah
 __m128 mad(__m128 const a, __m128 const b, __m128 const c) { return _mm_add_ps(_mm_mul_ps(a, b), c); } 


Mengapa ini perlu? Pertama-tama, untuk optimasi di masa depan. Misalnya, Anda cukup mengganti mad dengan fma dalam kode selesai melalui makro yang sama seperti yang Anda inginkan. Tapi kami akan meletakkan dasar untuk optimasi sekarang:

Varian dengan multiple-add
 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 . Lagi pecahan. Meski begitu, IACA terkadang menghasilkan hasil yang aneh. Kode tidak banyak berubah. Mungkin bahkan sedikit lebih buruk. Tetapi optimasi terkadang membutuhkan pengorbanan seperti itu.

Kami beralih ke menabung melalui _mm_stream


Berbagai panduan pengoptimalan merekomendasikan sekali lagi untuk tidak menarik cache untuk operasi penyimpanan massal. Ini biasanya dibenarkan ketika Anda memproses simpul yang ribuan atau lebih. Tetapi untuk matriks, ini mungkin tidak begitu penting. Namun, saya tetap akan menambahkannya.

Opsi penghematan aliran
 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]))); } 


Tidak ada yang berubah waktu di sini, dari kata sama sekali. Tetapi, menurut rekomendasi, kita sekarang tidak menyentuh cache sekali lagi.

Implementasi AVX


Opsi AVX dasar



Selanjutnya, kita beralih ke tahap optimasi berikutnya. Pelampung ke-4 termasuk dalam register SSE, dan dalam AVX sudah 8. Artinya, ada peluang teoritis untuk mengurangi jumlah operasi yang dilakukan, dan meningkatkan produktivitas jika tidak setengahnya, maka setidaknya 1,5 kali. Tetapi ada sesuatu yang memberitahu saya bahwa tidak semuanya akan menjadi sangat sederhana dengan transisi ke AVX. Bisakah kita mendapatkan data yang diperlukan dari register ganda?

Mari kita coba mencari tahu. Sekali lagi, kami menuliskan algoritma multiplikasi yang digunakan di atas. Anda tidak bisa melakukan ini, tetapi lebih mudah untuk menangani kode ketika semuanya ada di dekat Anda dan Anda tidak perlu menggulir ke atas setengah halaman.

 //    : 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 

Pada output, kami berharap mendapatkan hasilnya dalam ymm = {r0: r1} dan ymm = {r2: r3} . Jika dalam versi SSE, algoritma kami digeneralisasikan ke kolom, sekarang kami perlu menggeneralisasinya menjadi baris. Jadi untuk bertindak seperti dalam kasus opsi SSE tidak akan berfungsi.

Jika kita mempertimbangkan matriks m dalam register ymm , kita mendapatkan ymm = {m0: m1} dan ymm = {m2: m3}, masing-masing. Sebelumnya, kami hanya memiliki kolom matriks dalam register, dan sekarang kolom dan baris.

Jika Anda mencoba bertindak seperti sebelumnya, maka Anda harus mengalikan ymm = {m0: m1} dengan register ymm = {n00, n00, n00, n00}: {n10, n10, n10, n10, n10} . Karena n00 dan n01 berada di baris yang sama dari matriks n , menilai dengan set instruksi AVX yang tersedia, menyebarkannya dengan ymm akan mahal. Acak dan permutasi bekerja secara terpisah untuk masing-masing dari dua merangkak (tinggi dan rendah xmm ) di dalam register ymm .

Jika kita mengambil ymm dari matriks n , maka kita mendapatkan kedua elemen n00 dan n10 di tertinggi 2 xmm di dalam register ymm . {n00, n10, n20, n30}: {n01, n11, n21, n31} . Biasanya, indeks untuk instruksi yang ada adalah dari 0 hingga 3. Dan itu hanya alamat float dalam satu register xmm dari dua di dalam register ymm . Tidak mungkin mentransfer n10 dari xmm yang lebih lama ke yang lebih muda dengan murah . Dan kemudian fokus ini harus diulang beberapa kali. Kami tidak tahan dengan kehilangan langkah-langkah tersebut. Perlu untuk datang dengan sesuatu yang lain.

Kami dulu menggeneralisasi kolom, tetapi sekarang baris. Karena itu, kami akan mencoba untuk pergi dengan cara yang sedikit berbeda. Kita perlu mendapatkan hasilnya di {r0: r1} . Ini berarti bahwa algoritme harus ditingkatkan bukan dalam garis algoritma yang terpisah, tetapi dalam dua sekaligus. Dan di sini, apa yang minus dalam karya shuffle dan permutasi , akan menjadi nilai tambah bagi kita. Kita melihat apa yang akan kita miliki di register ymm ketika kita mempertimbangkan matriks n .

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

Ya, kami perhatikan bahwa di bagian xmm berbeda dari register ymm kami memiliki elemen 00 dan 01 . Mereka dapat dikalikan case dengan mendaftar melalui perintah permute di {_00, _00, _00, _00}: {_ 01, _01, _01, _01} , menunjukkan hanya satu indeks 3 untuk kedua bagian xmm . Inilah yang kita butuhkan. Memang, koefisien juga digunakan pada garis yang berbeda. Hanya sekarang di register ymm yang sesuai untuk perkalian, perlu untuk menjaga {m0: m0} , yaitu baris pertama yang diduplikasi dari matriks m .

Jadi, kami mengecat algoritme lebih detail. Kita membaca baris ganda dari matriks m dalam register ymm :

 mm[0] = {m0:m0} mm[1] = {m1:m1} mm[2] = {m2:m2} mm[3] = {m3:m3} 

Dan kemudian kita akan menghitung perkaliannya sebagai:

 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) 

Kami menulis ulang dengan lebih jelas:

 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> 

Atau dalam bentuk yang disederhanakan:

 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> 

Segalanya tampak jelas.

Tinggal menulis implementasi saja
 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); } 


Berikut adalah angka menarik dari IACA: x86 - 12.53, x64 - 12 . Meskipun, tentu saja, saya ingin yang lebih baik. Melewatkan sesuatu.

Optimasi AVX plus gula sintaksis


Tampaknya dalam kode di atas, AVX tidak digunakan secara maksimal. Kami menemukan bahwa alih-alih menetapkan dua garis identik dalam register ymm , kita dapat menggunakan siaran , yang dapat mengisi register ymm dengan dua nilai xmm identik. Juga, di sepanjang jalan, tambahkan beberapa "gula sintaksis" untuk fungsi AVX.

Peningkatan Implementasi 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])); } 


Dan di sini hasilnya sudah lebih menarik. IACA menghasilkan angka: x86 - 10, x64 - 8.58 , yang terlihat jauh lebih baik, tetapi masih belum 2 kali.

Opsi AVX + FMA (final)


Mari kita coba satu lagi. Sekarang akan logis untuk mengingat set instruksi FMA lagi, karena itu ditambahkan ke prosesor setelah AVX. Cukup ubah mul + individual tambahkan untuk satu operasi. Meskipun kami masih menggunakan instruksi perkalian untuk memberi kompiler lebih banyak peluang untuk optimisasi, dan prosesor untuk eksekusi paralel perkalian. Biasanya saya melihat kode yang dihasilkan di assembler untuk memastikan opsi mana yang lebih baik.

Dalam hal ini, kita perlu menghitung a * b + c * d + e * f + g * h . Anda dapat melakukan ini dahi: fma (a, b, fma (c, d, fma (e, f, g * h)))) . Tapi, seperti yang kita lihat, tidak mungkin untuk melakukan satu operasi di sini tanpa menyelesaikan yang sebelumnya. Dan ini berarti bahwa kita tidak akan dapat menggunakan kemampuan untuk melakukan penggandaan berpasangan, seperti yang dimungkinkan oleh pipa SIMD. Jika kita sedikit mengubah perhitungan fma (a, b, c * d) + fma (e, f, g * h) , kita akan melihat bahwa kita dapat memparalelkan perhitungan. Pertama lakukan dua perkalian independen, dan kemudian dua operasi fma independen.

Implementasi 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 . Sekarang sangat bagus. Seseorang mungkin akan mengatakan apa yang bisa dilakukan dengan lebih baik, tetapi saya tidak tahu caranya.

Tingkatan yang dicapai


Saya segera mencatat bahwa angka-angka ini tidak boleh dianggap sebagai kebenaran tertinggi. Bahkan dengan tes tetap, mereka berenang dalam batas-batas tertentu.Dan terlebih lagi mereka berperilaku berbeda pada platform yang berbeda. Dengan optimasi apa pun, lakukan pengukuran khusus untuk kasus Anda.

Daftar isi


  • Fungsi: nama fungsi. Berakhir pada s - fungsi dengan streaming, bergerak normal berbeda (tanpa streaming). Ditambahkan untuk kejelasan, karena ini cukup penting.
  • Siklus IACA: jumlah kutu per fungsi yang dihitung oleh IACA
  • Siklus terukur: jumlah langkah yang diukur (lebih sedikit lebih banyak)
  • IACA speedup: jumlah tindakan dalam garis nol / jumlah tindakan dalam satu garis
  • Speedup terukur: jumlah pengukuran di garis nol / jumlah pengukuran di garis (semakin banyak semakin baik)


Untuk loop_m, kutu dari artikel dikalikan dengan 64. Artinya, ini adalah nilai yang sangat perkiraan. Bahkan, ternyata seperti itu.

i3-3770:


x86
FungsiSiklus IACASiklus terukurIACA speedupSpeedup terukur
buka gulungan_m70.0050.751,001,00
loop_m233.60119.210,300,43
sse_v1m18.8927.513.701,84
sse_v2m19.0027.613.681,84
sse_v3m18.8927.223.701.86
sse_v4s18.8927.183.701.87
avx_v1m13.0019.215.382.64
avx_v1s13.0020.035.382.53
avx_v2m10.0012.916.993.93
avx_v2s10.0017.346.992.93

x64
FunctionIACA cyclesMeasured cyclesIACA speedupMeasured speedup
unroll_m7068.601.001.00
loop_m233.60119.370.300.57
sse_v1m18.8921.983.703.12
sse_v2m19.0021.093.683.25
sse_v3m18.8922.193.703.09
sse_v4s18.8922.393.703.06
avx_v1m13.009.615.387.13
avx_v1s13.0016.905.384.06
avx_v2m10.009.206.997.45
avx_v2s10.0014.646.994.68

i7-8700K:


x86
FunctionIACA cyclesMeasured cyclesIACA speedupMeasured speedup
unroll_m69.9540.251.001.00
loop_m233.6079.490.300.51
sse_v1m18.8919.313.702.09
sse_v2m19.0019.983.682.01
sse_v3m18.8919.693.702.04
sse_v4s18.8919.673.702.05
avx_v1m13.0014.225.382.83
avx_v1s13.0014.135.382.85
avx_v2m10.0011.736.993.43
avx_v2s10.0011.816.993.41
AVX+FMAm9.2110.387.603.88
AVX+FMAs9.2110.327.603.90

x64
FunctionIACA cyclesMeasured cyclesIACA speedupMeasured speedup
unroll_m69.9557.111.001.00
loop_m233.6075.730.300.75
sse_v1m18.8915.833.703.61
sse_v2m19.0017.223.683.32
sse_v3m18.8915.923.703.59
sse_v4s18.8916.183.703.53
avx_v1m13.007.035.388.12
avx_v1s13.0012.985.384.40
avx_v2m10.005.406.9910.57
avx_v2s10.0011.396.995.01
AVX+FMAm9.219.737.605.87
AVX+FMAs9.219.817.605.82

Kode uji dalam sumber. Jika ada saran yang masuk akal tentang cara memperbaikinya, tulis di komentar.

BONUS dari ranah fiksi


Sebenarnya, ini dari ranah fiksi karena jika saya melihat prosesor dengan dukungan untuk AVX512, maka mungkin dalam gambar. Namun, saya mencoba mengimplementasikan algoritma. Di sini saya tidak akan menjelaskan apa pun, analogi lengkap dengan AVX + FMA. Algoritmanya sama, hanya operasi yang lebih sedikit.

Seperti yang mereka katakan, saya hanya akan meninggalkannya di sini
 __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])); } 


Jumlahnya fantastis: x86 - 4,79, x64 - 5,42 (IACA dengan arsitektur SKX). Ini terlepas dari kenyataan bahwa algoritma ini memiliki 64 perkalian dan 48 tambahan.

Kode PS dari artikel




Ini adalah pengalaman pertama saya menulis artikel. Terima kasih atas komentar Anda. Mereka membantu membuat kode dan artikel menjadi lebih baik.

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


All Articles