Мне нужно написать функции умножения матрицы на вектор и матрицы на матрицу, но я не могу обернуться вокруг команд 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, но поддержка пока еще не так широко распространена.
Из этого можно построить умножение матрицы на матрицу (если вы меняете место, куда идет результат), но это не оптимально.
Других решений пока нет …