Как мне более эффективно умножить транспонирование матриц?
Я продолжаю бить себя по голове этим. У меня есть SSE-алгоритм для умножения матрицы A по матрице B, Мне нужно также реализовать операции, для которых A, B или оба транспонированы. Я сделал наивную реализацию, матричный код 4×4, представленный ниже (я думаю, что это довольно стандартные операции SSE), но A*B^T Операция занимает примерно вдвое дольше, чем A*B, Реализация ATLAS возвращает аналогичные значения для A*Bи почти идентичные результаты для умножения на транспонирование, что говорит мне о том, что есть эффективный способ сделать это.
MM-Умножение:
m1 = (mat1.m_>>2)<<2;
n2 = (mat2.n_>>2)<<2;
n = (mat1.n_>>2)<<2;
for (k=0; k<n; k+=4) {
for (i=0; i<m1; i+=4) {
// fetch: get 4x4 matrix from mat1
// row-major storage, so get 4 rows
Float* a0 = mat1.el_[i]+k;
Float* a1 = mat1.el_[i+1]+k;
Float* a2 = mat1.el_[i+2]+k;
Float* a3 = mat1.el_[i+3]+k;
for (j=0; j<n2; j+=4) {
// fetch: get 4x4 matrix from mat2
// row-major storage, so get 4 rows
Float* b0 = mat2.el_[k]+j;
Float* b1 = mat2.el_[k+1]+j;
Float* b2 = mat2.el_[k+2]+j;
Float* b3 = mat2.el_[k+3]+j;
__m128 b0r = _mm_loadu_ps(b0);
__m128 b1r = _mm_loadu_ps(b1);
__m128 b2r = _mm_loadu_ps(b2);
__m128 b3r = _mm_loadu_ps(b3);
{ // first row of result += first row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a0+0), b0r), _mm_mul_ps(_mm_load_ps1(a0+1), b1r));
__m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a0+2), b2r), _mm_mul_ps(_mm_load_ps1(a0+3), b3r));
Float* c0 = this->el_[i]+j;
_mm_storeu_ps(c0, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c0)));
}
{ // second row of result += second row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a1+0), b0r), _mm_mul_ps(_mm_load_ps1(a1+1), b1r));
__m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a1+2), b2r), _mm_mul_ps(_mm_load_ps1(a1+3), b3r));
Float* c1 = this->el_[i+1]+j;
_mm_storeu_ps(c1, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c1)));
}
{ // third row of result += third row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a2+0), b0r), _mm_mul_ps(_mm_load_ps1(a2+1), b1r));
__m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a2+2), b2r), _mm_mul_ps(_mm_load_ps1(a2+3), b3r));
Float* c2 = this->el_[i+2]+j;
_mm_storeu_ps(c2, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c2)));
}
{ // fourth row of result += fourth row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a3+0), b0r), _mm_mul_ps(_mm_load_ps1(a3+1), b1r));
__m128 cX2 = _mm_add_ps(_mm_mul_ps(_mm_load_ps1(a3+2), b2r), _mm_mul_ps(_mm_load_ps1(a3+3), b3r));
Float* c3 = this->el_[i+3]+j;
_mm_storeu_ps(c3, _mm_add_ps(_mm_add_ps(cX1, cX2), _mm_loadu_ps(c3)));
}
}
// Code omitted to handle remaining rows and columns
}
Для умножения MT (матрицы, умноженной на транспонированную матрицу) я заменяю сохраненные значения от b0r до b3r с помощью следующих команд и соответствующим образом изменяю переменные цикла:
__m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]);
__m128 b1r = _mm_set_ps(b3[1], b2[1], b1[1], b0[1]);
__m128 b2r = _mm_set_ps(b3[2], b2[2], b1[2], b0[2]);
__m128 b3r = _mm_set_ps(b3[3], b2[3], b1[3], b0[3]);
Я подозреваю, что замедление происходит частично из-за разницы между вытягиванием строки за раз и необходимостью сохранять 4 значения каждый раз, чтобы получить столбец, но я чувствую, что есть другой способ сделать это, вытягивая строки B и затем умножение на столбец As просто сместит стоимость на хранение 4 столбцов результатов.
Я также попытался вытянуть ряды B как ряды, а затем с помощью _MM_TRANSPOSE4_PS(b0r, b1r, b2r, b3r); сделать транспозицию (я думал, что в этом макросе могут быть некоторые дополнительные оптимизации), но реального улучшения нет.
На первый взгляд, я чувствую, что это должно быть быстрее … вовлеченные точечные продукты будут по порядку строка за строкой, что кажется более эффективным, но попытка сделать точечные продукты прямо вверх просто приводит к необходимости делать то же самое сохранить результаты.
Что мне здесь не хватает?
Добавлено: Просто чтобы уточнить, я стараюсь не транспонировать матрицы. Я бы предпочел перебирать их. Проблема, насколько я могу судить, состоит в том, что команда _mm_set_ps намного медленнее, чем _mm_load_ps.
Я также попробовал вариант, в котором я сохранил 4 строки матрицы A, а затем заменил 4 сегмента с фигурными скобками, содержащих 1 загрузку, 4 умножения и 2 добавления, с 4 инструкциями умножения и 3 hadds, но безрезультатно. Время остается неизменным (и да, я попробовал его с помощью оператора отладки, чтобы убедиться, что код изменился в моей тестовой компиляции. Указанный оператор отладки был удален перед профилированием, конечно):
{ // first row of result += first row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a0r, b0r), _mm_mul_ps(a0r, b1r));
__m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a0r, b2r), _mm_mul_ps(a0r, b3r));
Float* c0 = this->el_[i]+j;
_mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
}
{ // second row of result += second row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a1r, b0r), _mm_mul_ps(a1r, b1r));
__m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a1r, b2r), _mm_mul_ps(a1r, b3r));
Float* c0 = this->el_[i+1]+j;
_mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
}
{ // third row of result += third row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a2r, b0r), _mm_mul_ps(a2r, b1r));
__m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a2r, b2r), _mm_mul_ps(a2r, b3r));
Float* c0 = this->el_[i+2]+j;
_mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
}
{ // fourth row of result += fourth row of mat1 * 4x4 of mat2
__m128 cX1 = _mm_hadd_ps(_mm_mul_ps(a3r, b0r), _mm_mul_ps(a3r, b1r));
__m128 cX2 = _mm_hadd_ps(_mm_mul_ps(a3r, b2r), _mm_mul_ps(a3r, b3r));
Float* c0 = this->el_[i+3]+j;
_mm_storeu_ps(c0, _mm_add_ps(_mm_hadd_ps(cX1, cX2), _mm_loadu_ps(c0)));
}
Обновить:
Право и перемещение загрузки строк a0r в a3r в фигурные скобки в попытке избежать сбоя регистра также не удалось.
Решение
Я думаю, что это один из немногих случаев, когда горизонтальное добавление полезно. Вы хотите C = AB ^ T, но B не сохраняется в памяти для транспонирования. Это проблема. Это магазин как AoS вместо SoA. В этом случае, я думаю, что транспонирование B и вертикальное добавление выполняется медленнее, чем горизонтальное. Это по крайней мере верно для Матрицывектор Эффективное матричное векторное умножение 4×4 с SSE: горизонтальное сложение и точечный продукт — какой смысл?. В коде ниже функции m4x4 не является продуктом матрицы SSE 4×4, m4x4_vec использует SSE,
m4x4T С = АB ^ T без SSE и m4x4T_vec С = АB ^ T usisg SSE. Думаю, последний — тот, который тебе нужен.
Примечание: для больших матриц я бы не использовал этот метод. В этом случае быстрее сначала взять транспонирование и использовать вертикальное добавление (с SSE / AVX вы делаете что-то более сложное, вы транспонируете полосы с шириной SSE / AVX). Это потому, что транспонирование идет как O (n ^ 2), а матричное произведение идет как O (n ^ 3), поэтому для больших матриц транспонирование незначительно. Тем не менее, для 4×4 транспонирование является значительным, поэтому горизонтальное добавление выигрывает.
Редактировать:
Я неправильно понял, что вы хотели. Вы хотите C = (AБ) ^ Т. Это должно быть так же быстро, как (AБ) и код почти такой же, вы просто меняете роли А и В.
Мы можем написать математику следующим образом:
C = A*B in Einstein notation is C_i,j = A_i,k * B_k,j.
Since (A*B)^T = B^T*A^T we can write
C = (A*B)^T in Einstein notation is C_i,j = B^T_i,k * A^T_k,j = A_j,k * B_k,i
Если сравнить два, единственное, что меняется, мы поменяемся ролями j и i. Я поставил код для этого в конце этого ответа.
#include "stdio.h"#include <nmmintrin.h>
void m4x4(const float *A, const float *B, float *C) {
for(int i=0; i<4; i++) {
for(int j=0; j<4; j++) {
float sum = 0.0f;
for(int k=0; k<4; k++) {
sum += A[i*4+k]*B[k*4+j];
}
C[i*4 + j] = sum;
}
}
}
void m4x4T(const float *A, const float *B, float *C) {
for(int i=0; i<4; i++) {
for(int j=0; j<4; j++) {
float sum = 0.0f;
for(int k=0; k<4; k++) {
sum += A[i*4+k]*B[j*4+k];
}
C[i*4 + j] = sum;
}
}
}
void m4x4_vec(const float *A, const float *B, float *C) {
__m128 Brow[4], Mrow[4];
for(int i=0; i<4; i++) {
Brow[i] = _mm_load_ps(&B[4*i]);
}
for(int i=0; i<4; i++) {
Mrow[i] = _mm_set1_ps(0.0f);
for(int j=0; j<4; j++) {
__m128 a = _mm_set1_ps(A[4*i +j]);
Mrow[i] = _mm_add_ps(Mrow[i], _mm_mul_ps(a, Brow[j]));
}
}
for(int i=0; i<4; i++) {
_mm_store_ps(&C[4*i], Mrow[i]);
}
}
void m4x4T_vec(const float *A, const float *B, float *C) {
__m128 Arow[4], Brow[4], Mrow[4];
for(int i=0; i<4; i++) {
Arow[i] = _mm_load_ps(&A[4*i]);
Brow[i] = _mm_load_ps(&B[4*i]);
}
for(int i=0; i<4; i++) {
__m128 prod[4];
for(int j=0; j<4; j++) {
prod[j] = _mm_mul_ps(Arow[i], Brow[j]);
}
Mrow[i] = _mm_hadd_ps(_mm_hadd_ps(prod[0], prod[1]), _mm_hadd_ps(prod[2], prod[3]));
}
for(int i=0; i<4; i++) {
_mm_store_ps(&C[4*i], Mrow[i]);
}
}
float compare_4x4(const float* A, const float*B) {
float diff = 0.0f;
for(int i=0; i<4; i++) {
for(int j=0; j<4; j++) {
diff += A[i*4 +j] - B[i*4+j];
printf("A %f, B %f\n", A[i*4 +j], B[i*4 +j]);
}
}
return diff;
}
int main() {
float *A = (float*)_mm_malloc(sizeof(float)*16,16);
float *B = (float*)_mm_malloc(sizeof(float)*16,16);
float *C1 = (float*)_mm_malloc(sizeof(float)*16,16);
float *C2 = (float*)_mm_malloc(sizeof(float)*16,16);
for(int i=0; i<4; i++) {
for(int j=0; j<4; j++) {
A[i*4 +j] = i*4+j;
B[i*4 +j] = i*4+j;
C1[i*4 +j] = 0.0f;
C2[i*4 +j] = 0.0f;
}
}
m4x4T(A, B, C1);
m4x4T_vec(A, B, C2);
printf("compare %f\n", compare_4x4(C1,C2));
}
Редактировать:
Вот скаляр и функция SSE, которые делают C = (AБ) ^ Т. Они должны быть такими же быстрыми, как их АБ версии.
void m4x4TT(const float *A, const float *B, float *C) {
for(int i=0; i<4; i++) {
for(int j=0; j<4; j++) {
float sum = 0.0f;
for(int k=0; k<4; k++) {
sum += A[j*4+k]*B[k*4+i];
}
C[i*4 + j] = sum;
}
}
}
void m4x4TT_vec(const float *A, const float *B, float *C) {
__m128 Arow[4], Crow[4];
for(int i=0; i<4; i++) {
Arow[i] = _mm_load_ps(&A[4*i]);
}
for(int i=0; i<4; i++) {
Crow[i] = _mm_set1_ps(0.0f);
for(int j=0; j<4; j++) {
__m128 a = _mm_set1_ps(B[4*i +j]);
Crow[i] = _mm_add_ps(Crow[i], _mm_mul_ps(a, Arow[j]));
}
}
for(int i=0; i<4; i++) {
_mm_store_ps(&C[4*i], Crow[i]);
}
}
Другие решения
Несколько рекомендаций, которые могут помочь:
- Не используйте невыровненную память (эти _mm_loadu * работают медленно).
- Вы не обращаетесь к памяти последовательно, что убивает кеш. Попробуйте транспонировать матрицу до фактического доступа к этой памяти, это позволит процессору извлекать и максимально использовать кэш. Таким образом, следующее
__m128 b0r = _mm_set_ps(b3[0], b2[0], b1[0], b0[0]); // and b1r, etc..не нужен Идея состоит в том, чтобы захватить все 4 компонента последовательно. Если вам нужно реорганизовать память до вызова кода SSE, сделайте это. - Вы загружаете во внутренний цикл это:
_mm_load_ps1(a0+0)(то же самое для a1, a2 и a3) но является постоянным для всех итераций во внутреннем цикле. Вы можете загрузить эти значения снаружи и сохранить несколько циклов. Следите за тем, что вы можете использовать из предыдущей итерации. - Профиль. использовать Intel VTune или что-то подобное, он скажет вам, где узкое место.