Матрично-векторное и матрично-матричное умножение с использованием SSE
Мне нужно написать функции умножения матрицы на вектор и матрицы на матрицу, но я не могу обернуться вокруг команд SSE.
Размеры матриц и векторов всегда кратны 4.
Мне удалось написать функцию умножения вектора на вектор, которая выглядит так:
void vector_multiplication_SSE(float* m, float* n, float* result, unsigned const int size)
{
int i;
__declspec(align(16))__m128 *p_m = (__m128*)m;
__declspec(align(16))__m128 *p_n = (__m128*)n;
__declspec(align(16))__m128 *p_result = (__m128*)result;
for (i = 0; i < size / 4; ++i)
p_result[i] = _mm_mul_ps(p_m[i], p_n[i]);
// print the result
for (int i = 0; i < size; ++i)
{
if (i % 4 == 0) cout << endl;
cout << result[i] << '\t';
}
}
и теперь я пытаюсь реализовать умножение матрицы на вектор.
Вот что у меня так далеко:
void multiply_matrix_by_vector_SSE(float* m, float* v, float* result, unsigned const int vector_dims)
{
int i, j;
__declspec(align(16))__m128 *p_m = (__m128*)m;
__declspec(align(16))__m128 *p_v = (__m128*)v;
__declspec(align(16))__m128 *p_result = (__m128*)result;
for (i = 0; i < vector_dims; i += 4)
{
__m128 tmp = _mm_load_ps(&result[i]);
__m128 p_m_tmp = _mm_load_ps(&m[i]);
tmp = _mm_add_ps(tmp, _mm_mul_ps(tmp, p_m_tmp));
_mm_store_ps(&result[i], tmp);
// another for loop here?
}
// print the result
for (int i = 0; i < vector_dims; ++i)
{
if (i % 4 == 0) cout << endl;
cout << result[i] << '\t';
}
}
Эта функция выглядит совершенно неправильно. Я имею в виду не только то, что он работает неправильно, но и то, что я двигаюсь не в том направлении.
Может ли кто-нибудь помочь мне с реализацией векторно-матричного и матрично-матричного умножения? Я был бы очень признателен за пример кода и очень подробное объяснение.
Вот моя попытка № 2:
это терпит неудачу с Access reading violation исключение, но все еще чувствует себя ближе
void multiply_matrix_by_vector_SSE(float* m, float* v, float* result, unsigned const int vector_dims)
{
int i, j;
__declspec(align(16))__m128 *p_m = (__m128*)m;
__declspec(align(16))__m128 *p_v = (__m128*)v;
__declspec(align(16))__m128 *p_result = (__m128*)result;
for (i = 0; i < vector_dims; ++i)
{
p_result[i] = _mm_mul_ps(_mm_load_ps(&m[i]), _mm_load_ps1(&v[i]));
}
// print the result
for (int i = 0; i < vector_dims; ++i)
{
if (i % 4 == 0) cout << endl;
cout << result[i] << '\t';
}
}
void multiply_matrix_by_vector_SSE(float* m, float* v, float* result, unsigned const int vector_dims)
{
int i, j;
__declspec(align(16))__m128 *p_m = (__m128*)m;
__declspec(align(16))__m128 *p_v = (__m128*)v;
__declspec(align(16))__m128 *p_result = (__m128*)result;
for (i = 0; i < vector_dims; ++i)
{
for (j = 0; j < vector_dims * vector_dims / 4; ++j)
{
p_result[i] = _mm_mul_ps(p_v[i], p_m[j]);
}
}
for (int i = 0; i < vector_dims; ++i)
{
if (i % 4 == 0) cout << endl;
cout << result[i] << '\t';
}
cout << endl;
}
Решение
Без всяких ухищрений или умножения матрицы на вектор это просто набор точечных произведений между вектором и строкой матрицы. Ваш код на самом деле не имеет такой структуры. Написание это на самом деле как точечные продукты (не проверено):
for (int row = 0; row < nrows; ++row) {
__m128 acc = _mm_setzero_ps();
// I'm just going to assume the number of columns is a multiple of 4
for (int col = 0; col < ncols; col += 4) {
__m128 vec = _mm_load_ps(&v[col]);
// don't forget it's a matrix, do 2d addressing
__m128 mat = _mm_load_ps(&m[col + ncols * row]);
acc = _mm_add_ps(acc, _mm_mul_ps(mat, vec));
}
// now we have 4 floats in acc and they have to be summed
// can use two horizontal adds for this, they kind of suck but this
// isn't the inner loop anyway.
acc = _mm_hadd_ps(acc, acc);
acc = _mm_hadd_ps(acc, acc);
// store result, which is a single float
_mm_store_ss(&result[row], acc);
}
Есть несколько очевидных приемов, таких как обработка нескольких строк одновременно, повторное использование нагрузки из вектора и создание нескольких независимых цепочек зависимостей, чтобы вы могли лучше использовать пропускную способность. Также очень простой прием — использовать FMA для комбо mul / add, но поддержка пока еще не так широко распространена.
Из этого можно построить умножение матрицы на матрицу (если вы меняете место, куда идет результат), но это не оптимально.
Другие решения
Других решений пока нет …