Матрица-вектор и умножение матрицы-матрицы с использованием SSE

Question

Мне нужно написать функции умножения матрицы-матрицы и матрицы-матрицы, но я не могу обернуть голову вокруг команд 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'; 
    } 
} 

Эта функция выглядит совершенно неправильно. Я имею в виду, что не только это работает неправильно, но также кажется, что я двигаюсь в неправильном направлении.

---

Может ли кто-нибудь помочь мне с внедрением матриц матричных матриц и матриц? Я бы очень признателен за кусок кода примера и очень подробное объяснение

Update

Вот моя попытка номер 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'; 
    } 
} 

Обновление 2

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; 
} 

Answer 1

Сено ut любые трюки или что-либо, умножение матрицы-вектора - это всего лишь куча точечных произведений между вектором и строкой матрицы. В вашем коде нет такой структуры. Дать это на самом деле, как продукты точечно (не проверено):

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, но поддержка пока не распространена.

Вы можете построить матричное матричное умножение из этого (если вы измените место, в котором результат), но это не оптимально.