我需要编写矩阵向量和矩阵矩阵乘法函数，但是我无法将头围在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组合，但是支持还不那么广泛（虽然不是在2015年，但现在在2020年相当普遍）。
您可以据此建立矩阵-矩阵乘法（如果您改变结果所在的位置），但这并不是最佳选择（请参阅下文）。>
一次处理四行（未测试）：
for (int row = 0; row < nrows; row += 4) { __m128 acc0 = _mm_setzero_ps(); __m128 acc1 = _mm_setzero_ps(); __m128 acc2 = _mm_setzero_ps(); __m128 acc3 = _mm_setzero_ps(); for (int col = 0; col < ncols; col += 4) { __m128 vec = _mm_load_ps(&v[col]); __m128 mat0 = _mm_load_ps(&m[col + ncols * row]); __m128 mat1 = _mm_load_ps(&m[col + ncols * (row + 1)]); __m128 mat2 = _mm_load_ps(&m[col + ncols * (row + 2)]); __m128 mat3 = _mm_load_ps(&m[col + ncols * (row + 3)]); acc0 = _mm_add_ps(acc0, _mm_mul_ps(mat0, vec)); acc1 = _mm_add_ps(acc1, _mm_mul_ps(mat1, vec)); acc2 = _mm_add_ps(acc2, _mm_mul_ps(mat2, vec)); acc3 = _mm_add_ps(acc3, _mm_mul_ps(mat3, vec)); } acc0 = _mm_hadd_ps(acc0, acc1); acc2 = _mm_hadd_ps(acc2, acc3); acc0 = _mm_hadd_ps(acc0, acc2); _mm_store_ps(&result[row], acc0); }现在，每4个FMA只有5个负载，而在未行展开的版本中，每1个FMA有2个负载。另外，还有4个独立的FMA，或没有FMA收缩的添加/多对，这两种方式都增加了流水线/同时执行的可能性。实际上，您可能想展开得更多，例如Skylake可以每个周期启动2个独立FMA，并且它们需要4个周期才能完成，因此要完全占用两个FMA单元，您需要8个独立FMA。作为奖励，最后这3个水平加法相对较好地用于水平求和。
最初，不同的数据布局似乎是一个缺点，不再可能简单地从矩阵和向量中进行向量加载并将它们相乘（这会使第一个矩阵的微小行向量乘以微小的[[行
第二个矩阵的向量，这是错误的）。但是，完全矩阵-矩阵乘法可以利用以下事实：它实际上是将矩阵与许多独立的向量相乘，这充满了需要完成的独立工作。水平和也可以很容易地避免。因此，实际上它比矩阵向量乘法更方便。关键是从矩阵A提取一点列向量，从矩阵B提取一点行向量，然后将它们相乘成一个小的矩阵。与您惯常的习惯相比，这听起来可能是相反的，但是通过SIMD这样做会更好，因为计算始终保持独立且无水平操作。例如（未经测试，假设矩阵的尺寸可被展开因子整除，需要x64，否则将用尽寄存器）
for (size_t i = 0; i < mat1rows; i += 4) {
for (size_t j = 0; j < mat2cols; j += 8) {
float* mat1ptr = &mat1[i * mat1cols];
__m256 sumA_1, sumB_1, sumA_2, sumB_2, sumA_3, sumB_3, sumA_4, sumB_4;
sumA_1 = _mm_setzero_ps();
sumB_1 = _mm_setzero_ps();
sumA_2 = _mm_setzero_ps();
sumB_2 = _mm_setzero_ps();
sumA_3 = _mm_setzero_ps();
sumB_3 = _mm_setzero_ps();
sumA_4 = _mm_setzero_ps();
sumB_4 = _mm_setzero_ps();
for (size_t k = 0; k < mat2rows; ++k) {
auto bc_mat1_1 = _mm_set1_ps(mat1ptr[0]);
auto vecA_mat2 = _mm_load_ps(mat2 + m2idx);
auto vecB_mat2 = _mm_load_ps(mat2 + m2idx + 4);
sumA_1 = _mm_add_ps(_mm_mul_ps(bc_mat1_1, vecA_mat2), sumA_1);
sumB_1 = _mm_add_ps(_mm_mul_ps(bc_mat1_1, vecB_mat2), sumB_1);
auto bc_mat1_2 = _mm_set1_ps(mat1ptr[N]);
sumA_2 = _mm_add_ps(_mm_mul_ps(bc_mat1_2, vecA_mat2), sumA_2);
sumB_2 = _mm_add_ps(_mm_mul_ps(bc_mat1_2, vecB_mat2), sumB_2);
auto bc_mat1_3 = _mm_set1_ps(mat1ptr[N * 2]);
sumA_3 = _mm_add_ps(_mm_mul_ps(bc_mat1_3, vecA_mat2), sumA_3);
sumB_3 = _mm_add_ps(_mm_mul_ps(bc_mat1_3, vecB_mat2), sumB_3);
auto bc_mat1_4 = _mm_set1_ps(mat1ptr[N * 3]);
sumA_4 = _mm_add_ps(_mm_mul_ps(bc_mat1_4, vecA_mat2), sumA_4);
sumB_4 = _mm_add_ps(_mm_mul_ps(bc_mat1_4, vecB_mat2), sumB_4);
m2idx += 8;
mat1ptr++;
}
_mm_store_ps(&result[i * mat2cols + j], sumA_1);
_mm_store_ps(&result[i * mat2cols + j + 4], sumB_1);
_mm_store_ps(&result[(i + 1) * mat2cols + j], sumA_2);
_mm_store_ps(&result[(i + 1) * mat2cols + j + 4], sumB_2);
_mm_store_ps(&result[(i + 2) * mat2cols + j], sumA_3);
_mm_store_ps(&result[(i + 2) * mat2cols + j + 4], sumB_3);
_mm_store_ps(&result[(i + 3) * mat2cols + j], sumA_4);
_mm_store_ps(&result[(i + 3) * mat2cols + j + 4], sumB_4);
}
}
该代码的重点是易于将计算安排为非常易于SIMD使用，并具有许多独立的算法来使浮点单元饱和，同时使用较少的负载（否则可能会变得一个瓶颈，甚至抛开它们可能会错过L1缓存，只是因为它们太多了。]
您甚至可以使用此代码，但是它与Intel MKL相比没有竞争力。特别是对于中型或大型矩阵，其中平铺极为重要。将其升级到AVX很容易。它根本不适合微小的矩阵，例如将两个4x4矩阵相乘，请参见Efficient 4x4 matrix multiplication。