Вычисление матричного произведения намного медленнее с SSE, чем с прямым алгоритмом
Я хочу умножить две матрицы, один раз, используя прямой алгоритм:
template <typename T>
void multiplicate_straight(T ** A, T ** B, T ** C, int sizeX)
{
T ** D = AllocateDynamicArray2D<T>(sizeX, sizeX);
transpose_matrix(B, D,sizeX);
for(int i = 0; i < sizeX; i++)
{
for(int j = 0; j < sizeX; j++)
{
for(int g = 0; g < sizeX; g++)
{
C[i][j] += A[i][g]*D[j][g];
}
}
}
FreeDynamicArray2D<T>(D);
}
и один раз через использование функций SSE. Для этого я создал две функции:
template <typename T>
void SSE_vectormult(T * A, T * B, int size)
{
__m128d a;
__m128d b;
__m128d c;
#ifdef linux
double A2[2], B2[2], C[2] __attribute__ ((aligned(16)));
#endif
#ifdef _WIN32
__declspec(align(16)) double A2[2], B2[2], C[2];
#endif
for(int i = 0; i < size; i+=2)
{
//std::cout << "In SSE_vectormult: i is: " << i << '\n';
A2[0] = A[i];
B2[0] = B[i];
A2[1] = A[i+1];
B2[1] = B[i+1];
//std::cout << "Values from A and B written to A2 and B2\n";
a = _mm_load_pd(A2);
b = _mm_load_pd(B2);
//std::cout << "Values converted to a and b\n";
c = _mm_mul_pd(a,b);
_mm_store_pd(C, c);
A[i] = C[0];
A[i+1] = C[1];
};
}
а также
template <typename T>
void multiplicate_SSE(T ** A, T ** B, T ** C, int sizeX)
{
// std::cout << "Entered SSE-Function\n";
T ** D = AllocateDynamicArray2D<T>(sizeX, sizeX);
T * tmp = AllocateDynamicArray1D<T>(sizeX);
T * tmp2 = AllocateDynamicArray1D<T>(sizeX);
//std::cout << "Matrices allocated\n";
transpose_matrix<T>(B, D,sizeX);
//std::cout << "Matrix B transposed\n";
for(int i = 0; i < sizeX; i++)
{
for(int j = 0; j < sizeX; j++)
{
extract_row<T>(A,tmp, i, sizeX);
// std::cout << "Row from A extracted\n";
//print_vector(tmp, sizeX);
extract_row<T>(D, tmp2, j, sizeX);
// std::cout << "Row from D extracted\n";
//print_vector(tmp2, sizeX);
SSE_vectormult<T>(tmp, tmp2, sizeX);
// std::cout << "Vectors multiplicated\n";
//print_vector(tmp, sizeX);
C[i][j] = add_vector(tmp, sizeX);
// std::cout << "Written value to C\n";
// std::cout << "j is " << j << " and i is " << i << '\n';
}
}
// std::cout << "Loop finished\n";
FreeDynamicArray2D<T>(D);
//std::cout << "Freed D\n";
//FreeDynamicArray1D<T>(tmp);????
// std::cout << "Freed tmp\n";
FreeDynamicArray1D<T>(tmp2);
// std::cout << "Everything freed, returning\n";
}
Но тогда у меня возникает несколько проблем: С одной стороны, когда я хочу освободить массив tmp в multiplicate_SSE(), отмеченном несколькими вопросительными знаками, я получаю ошибку "_BLOCK_TYPE_IS_VALID". Я думал о возможности освобождения одного и того же пространства дважды, поэтому я раскомментировал это (но, как я полагаю, я получаю утечку памяти из-за этого?). Теперь, когда я сравниваю производительность обеих функций с одинаковыми матрицами, функция SSE требует примерно в четыре раза больше для двух матриц 1024x1024, чем прямой метод.
Как я могу переписать мою SSE-функцию, чтобы получить лучшую производительность (я никогда раньше не работал с SSE), и как я могу исправить утечку памяти?
Спасибо!
2 ответа
У вас есть правильная идея, взяв транспонирование в скалярный код, но вы не хотите точно транспонировать при использовании SSE.
Давайте придерживаться, чтобы плавать (SGEMM). То, что вы хотите сделать с SSE, - это делать четыре продукта одновременно. Ты хочешь C = A*B, Давайте посмотрим на матрицу 8x8. Давайте предположим B является:
(0 1 2 3) ( 4 5 6 7)
(8 9 10 11) (12 13 14 15)
(16 17 18 19) (20 21 22 23)
(24 25 26 27) (28 29 30 31)
(32 33 34 35) (36 37 38 39)
(40 41 42 43) (44 45 46 47)
(48 49 50 51) (52 53 54 55)
(56 57 58 59) (60 61 62 63)
Итак, с SSE вы делаете:
C[0][0] C[0][1] C[0][2] C[0][3] =
A[0][0]*(0 1 2 3) + A[0][1]*(8 9 10 11) + A[0][2]*(16 17 18 19)...+ A[0][7]*(56 57 58 59)
Это дает вам четыре продукта сразу. Проблема в том, что вы должны двигаться вниз по столбцу в B и значения не находятся в одной строке кэша. Было бы лучше, если бы каждый столбец шириной четыре был непрерывным в памяти. Таким образом, вместо транспонирования каждого элемента вы перемещаете полосы шириной 4 примерно так:
(0 1 2 3)( 8 9 10 11)(16 17 18 19)(24 25 26 27)(32 33 34 35)(40 41 42 43)(48 49 50 51)(56 57 58 59)
(4 5 6 7)(12 13 14 15)(20 21 22 23)(28 29 30 31)(36 37 38 39)(44 45 46 47)(52 53 54 55)(60 61 62 63)
Если вы считаете каждое из четырех значений в скобках одной единицей, это эквивалентно переносу матрицы 8x2 в матрицу 2x8. Обратите внимание, что столбцы шириной четыре B смежны в памяти. Это гораздо более дружественно к кешу. Для матрицы 8x8 это не проблема, но, например, для матрицы 1024x1024. Смотрите код ниже, чтобы узнать, как это сделать. Для AVX вы переставляете полосы шириной 8 (что означает, что вам нечего делать с матрицей 8x8). Для удвоения ширина составляет два с SSE и четыре с AVX.
Это должно быть в четыре раза быстрее, чем скалярный код, при условии, что матрицы помещаются в кэш. Однако для больших матриц этот метод все еще будет ограничен памятью, и поэтому ваш код SSE может быть не намного быстрее, чем скалярный код (но он не должен быть хуже).
Однако, если вы используете мозаичное построение циклов и перестановку матрицы в мозаичных элементах (которые помещаются в кэш L2), а не для всей матрицы, умножение матрицы ограничивается вычислением, а не памятью даже для очень больших матриц, которые не помещаются в кэш L3., Это другая тема.
Изменить: некоторый (непроверенный) код для сравнения с вашим скалярным кодом. Я развернул петлю на 2.
void SGEMM_SSE(const float *A, const float *B, float *C, const int sizeX) {
const int simd_width = 4;
const int unroll = 2;
const int strip_width = simd_width*unroll
float *D = (float*)_mm_malloc(sizeof(float)*sizeX*sizeX, 16);
transpose_matrix_strip(B, D,sizeX, strip_width); //tranpose B in strips of width eight
for(int i = 0; i < sizeX; i++) {
for(int j = 0; j < sizeX; j+=strip_width) {
float4 out_v1 = 0; //broadcast (0,0,0,0)
float4 out_V2 = 0;
//now calculate eight dot products
for(int g = 0; g < sizeX; g++) {
//load eight values rrom D into two SSE registers
float4 vec4_1.load(&D[j*sizeX + strip_width*g]);
float4 vec4_2.load(&D[j*sizeX + strip_width*g + simd_width]);
out_v1 += A[i][g]*vec4_v1;
out_v2 += A[i][g]*vec4_v2;
}
//store eight dot prodcuts into C
out_v1.store(&C[i*sizeX + j]);
out_v2.store(&C[i*sizeX + j + simd_width]);
}
}
_mm_free(D);
}
void transpose_matrix_strip(const float* A, float* B, const int N, const int strip_width) {
//#pragma omp parallel for
for(int n=0; n<N*N; n++) {
int k = strip_width*(n/N/strip_width);
int i = (n/strip_width)%N;
int j = n%strip_width;
B[n] = A[N*i + k + j];
}
}
Заметить, что j увеличивается на восемь сейчас. Более раскрутка может помочь. Если вы хотите использовать встроенные функции, вы можете использовать _mm_load_ps, _mm_store_ps, _mm_set1_ps (например, для трансляций _mm_set1_ps(A[i][g])), _mm_add_ps, а также _mm_mul_ps, Вот и все.
Я считаю, что это должно делать то же самое, что и первый цикл с SSE, при условии, что sizeX кратен двум, а память выровнена по 16 байтов.
Вы можете повысить производительность, развернув цикл и используя несколько временных переменных, которые вы добавляете в конце. Вы также можете попробовать AVX и новую инструкцию Fused Multiply Add.
template <typename T>
void multiplicate_SSE2(T ** A, T ** B, T ** C, int sizeX)
{
T ** D = AllocateDynamicArray2D<T>(sizeX, sizeX);
transpose_matrix(B, D,sizeX);
for(int i = 0; i < sizeX; i++)
{
for(int j = 0; j < sizeX; j++)
{
__m128d temp = _mm_setzero_pd();
for(int g = 0; g < sizeX; g += 2)
{
__m128d a = _mm_load_pd(&A[i][g]);
__m128d b = _mm_load_pd(&D[j][g]);
temp = _mm_add_pd(temp, _mm_mul_pd(a,b));
}
// Add top and bottom half of temp together
temp = _mm_add_pd(temp, _mm_shuffle_pd(temp, temp, 1));
_mm_store_sd(temp, &C[i][j]); // Store one value
}
}
FreeDynamicArray2D<T>(D);
}