Ответ 1
У вас есть правильная идея, взяв транспонирование в скалярном коде, но вы не хотите точно транспонировать при использовании SSE.
Пусть придерживается float (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), а не для всего матричного умножения матрицы, это вычисление связано, а не связано с памятью даже для очень больших матриц, t в кеш-памяти 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 теперь увеличивается на восемь. Больше разворачивания может помочь. Если вы хотите использовать intrinsics, вы можете использовать _mm_load_ps, _mm_store_ps, _mm_set1_ps (для трансляций, например _mm_set1_ps(A[i][g])), _mm_add_ps и _mm_mul_ps. Что это.