SSE умножение матрицы на матрицу

Question

У меня возникают проблемы, делая умножение матрицы на матрицу с SSE в С.

Вот что я получил до сих пор:

#define N 1000 

void matmulSSE(int mat1[N][N], int mat2[N][N], int result[N][N]) { 
    int i, j, k; 
    __m128i vA, vB, vR; 

    for(i = 0; i < N; ++i) { 
    for(j = 0; j < N; ++j) { 
     vR = _mm_setzero_si128(); 
     for(k = 0; k < N; k += 4) { 
      //result[i][j] += mat1[i][k] * mat2[k][j]; 
      vA = _mm_loadu_si128((__m128i*)&mat1[i][k]); 
      vB = _mm_loadu_si128((__m128i*)&mat2[k][j]); //how well does the k += 4 work here? Should it be unrolled? 
      vR = _mm_add_epi32(vR, _mm_mul_epi32(vA, vB)); 
     } 
     vR = _mm_hadd_epi32(vR, vR); 
     vR = _mm_hadd_epi32(vR, vR); 
     result[i][j] += _mm_extract_epi32(vR, 0); 
    } 
    } 
} 

Я не могу это сделать дать правильные результаты. Я что-то упускаю? и поиск доцент, кажется, очень помогает - каждый результат либо только делает матрицы 4x4, мат-VEC или какая-то особая магия то не очень читаемый и трудно понять ...

Update: Woho! Я наконец-то понял. Помимо ошибок в моей логике (спасибо за помощь Питера Кордеса), также возникла проблема с _mm_mul_epi32(), которая не работает, как я и думал, - я должен был использовать _mm_mullo_epi32() вместо этого!

Я знаю, что это не самый эффективный код, но он был создан, чтобы заставить его работать должным образом в первую очередь - теперь я могу перейти к его оптимизации.

void matmulSSE(int mat1[N][N], int mat2[N][N], int result[N][N]) { 
    int i, j, k; 
    __m128i vA, vB, vR, vSum; 

    for(i = 0; i < N; ++i) { 
     for(j = 0; j < N; ++j) { 
      vR = _mm_setzero_si128(); 
      for(k = 0; k < N; k += 4) { 
       //result[i][j] += mat1[i][k] * mat2[k][j]; 
       vA = _mm_loadu_si128((__m128i*)&mat1[i][k]); 
       vB = _mm_insert_epi32(vB, mat2[k][j], 0); 
       vB = _mm_insert_epi32(vB, mat2[k + 1][j], 1); 
       vB = _mm_insert_epi32(vB, mat2[k + 2][j], 2); 
       vB = _mm_insert_epi32(vB, mat2[k + 3][j], 3); 
       vR = _mm_mullo_epi32(vA, vB); 
       vR = _mm_hadd_epi32(vR, vR); 
       vR = _mm_hadd_epi32(vR, vR); 
       result[i][j] += _mm_extract_epi32(vR, 0); 

       //DEBUG 
       //printf("vA: %d, %d, %d, %d\n", vA.m128i_i32[0], vA.m128i_i32[1], vA.m128i_i32[2], vA.m128i_i32[3]); 
       //printf("vB: %d, %d, %d, %d\n", vB.m128i_i32[0], vB.m128i_i32[1], vB.m128i_i32[2], vB.m128i_i32[3]); 
       //printf("vR: %d, %d, %d, %d\n", vR.m128i_i32[0], vR.m128i_i32[1], vR.m128i_i32[2], vR.m128i_i32[3]); 
       //printf("\n"); 
      } 
     } 
    } 
} 

Обновление 2: конвертированы пример Peters к I-к-J версии порядка цикла. Требуется дополнительная нагрузка для vR и перемещение в хранилище во внутренний цикл, но установка vA может быть перемещена по петле. Выключилось быстрее.

void matmulSSE_2(int mat1[N][N], int mat2[N][N], int result[N][N]) { 
    int i, j, k; 
    __m128i vA, vB, vR; 

    for(i = 0; i < N; ++i) { 
     for(k = 0; k < N; ++k) { 
      vA = _mm_set1_epi32(mat1[i][k]); 
      for(j = 0; j < N; j += 4) { 
       //result[i][j] += mat1[i][k] * mat2[k][j]; 
       vB = _mm_loadu_si128((__m128i*)&mat2[k][j]); 
       vR = _mm_loadu_si128((__m128i*)&result[i][j]); 
       vR = _mm_add_epi32(vR, _mm_mullo_epi32(vA, vB)); 
       _mm_storeu_si128((__m128i*)&result[i][j], vR); 

       //DEBUG 
       //printf("vA: %d, %d, %d, %d\n", vA.m128i_i32[0], vA.m128i_i32[1], vA.m128i_i32[2], vA.m128i_i32[3]); 
       //printf("vB: %d, %d, %d, %d\n", vB.m128i_i32[0], vB.m128i_i32[1], vB.m128i_i32[2], vB.m128i_i32[3]); 
       //printf("vR: %d, %d, %d, %d\n", vR.m128i_i32[0], vR.m128i_i32[1], vR.m128i_i32[2], vR.m128i_i32[3]); 

       //printf("\n"); 
      } 
     } 
    } 
} 

Answer 1

Вы правы, ваша vB проблема. Вы загружаете 4 последовательных целых числа, но mat2[k+0..3][j] не смежны. Вы на самом деле получаете mat2[k][j+0..3].

---

Я забыл, что хорошо работает для matmul. Иногда он работает хорошо, чтобы производить 4 результата параллельно, вместо того, чтобы делать горизонтальную сумму для каждого результата.

Транспонирование одной из ваших входных матриц хорошо работает и стоит O (N^2). Это стоит того, потому что это означает, что макет O (N^3) может использовать последовательный доступ, а ваша текущая структура цикла становится дружественной к SIMD.

Есть еще лучшие способы, переносящие небольшие блоки прямо перед использованием, поэтому они все еще горят в кеше L1, когда вы читаете их снова. Блокировка кэша, например, черепица, является одним из ключевых показателей хорошей производительности матрицы.

Много написано об оптимизации умножения матрицы, с SIMD и с блокировкой кеша. Я предлагаю вам это сделать. Большинство, если это, вероятно, говорит о FP, но все это относится и к целому числу.

(за исключением, что SSE/AVX имеет только FMA для FP, а не для 32-битных чисел, и 8 и 16-битные инструкции входные PMADD сделать горизонтальный добавляет пар.)

---

На самом деле я думаю, что вы можете произвести 4 результатов параллельно здесь:

void matmulSSE(int mat1[N][N], int mat2[N][N], int result[N][N]) { 

    for(int i = 0; i < N; ++i) { 
    for(int j = 0; j < N; j+=4) { // vectorize over this loop 
     __m128i vR = _mm_setzero_si128(); 
     for(int k = 0; k < N; k++) { // not this loop 
      //result[i][j] += mat1[i][k] * mat2[k][j]; 
      __m128i vA = _mm_set1_epi32(mat1[i][k]); // load+broadcast is much cheaper than MOVD + 3 inserts (or especially 4x insert, which your new code is doing) 
      __m128i vB = _mm_loadu_si128((__m128i*)&mat2[k][j]); // mat2[k][j+0..3] 
      vR = _mm_add_epi32(vR, _mm_mullo_epi32(vA, vB)); 
     } 
     _mm_storeu_si128((__m128i*)&result[i][j], vR)); 
    } 
    } 
} 

широковещательной передачи нагрузки (или отдельный груз + радиопередачу без AVX) по-прежнему намного дешевле, чем собрать.

Ваш текущий код собирает с 4 вставками вместо того, чтобы разрушать цепочку зависимостей по предыдущему значению итерации, используя MOVD для первого элемента, так что это еще хуже.Но даже лучшая сборка из 4 рассеянных элементов довольно плохая по сравнению с нагрузкой + PSHUFD. Не говоря уже о том, что для вашего кода нужен SSE4.1. Хотя это так или иначе, для _mm_mullo_epi32 вместо расширения PMULDQ (_mm_mul_epi32).