neon指令集加速arm架构下的浮点型矩阵乘法

前言

最近遇到一个问题, fpga进行卷积后的结果需要乘上一个pca矩阵, 本来打算用fpga做这个矩阵乘法的,考虑到量化后会有精度损失,打算先用arm实现,如果效果还不错的话,那就不用fpga做矩阵乘法。下面是c++写矩阵乘法的优化过程。

约定:
矩阵M(1x2048) : fpga_output_float_value;
矩阵N(2048x128) : pca_array_row / pca_array_col , 按行/列存储, 矩阵的内容一致
v1版本: 直接按照矩阵乘法的定义编写代码,所写即所想
std::vector<float> matrix_multi_result;
for (int col = 0; col < 128; col++)
{
float temp = 0.0;
for (int row = 0; row < 2048; row++)
{
temp += fpga_output_float_value[row] * pca_array_row[row * 128 + col];
}

matrix_multi_result.push_back(temp);
}
v2版本: 改写for循环使得sum[0], sunm[1], sum[2], sum[3]各自运算过程相互独立,让编译器自动优化, 记得编译选项中加上优化参数 -O3
float sum[4] = {0.0, 0.0, 0.0, 0.0};
std::vector<float> matrix_multi_result;
for (int col = 0; col < 128; ++col)
{
float temp = 0.0;
for (int row = 0; row < 2048; row += 4)
{
sum[0] = fpga_output_float_value[row + 0] * pca_array_row[(row + 0) * 128 + i];
sum[1] = fpga_output_float_value[row + 1] * pca_array_row[(row + 1) * 128 + i];
sum[2] = fpga_output_float_value[row + 2] * pca_array_row[(row + 2) * 128 + i];
sum[3] = fpga_output_float_value[row + 3] * pca_array_row[(row + 3) * 128 + i];

temp += sum[0] + sum[1] + sum[2] + sum[3];
}

matrix_multi_result.push_back(temp);
}
v3版本: 使用neon指令集优化运算过程
// neon 指令集实现的"乘加"过程
float dot_product_intrinsic(float * __restrict vec1, float * __restrict vec2, int length)
{
float32x4_t vec1_q, vec2_q;
float32x4_t sum_q = {0.0, 0.0, 0.0, 0.0};
float32x2_t tmp[2];
float result = 0.0;
for (int i = 0; i < (length & ~3); i += 4)
{
vec1_q = vld1q_f32(&vec1[i]);
vec2_q = vld1q_f32(&vec2[i]);
sum_q = vmlaq_f32(sum_q, vec1_q, vec2_q);
}
tmp[0] = vget_high_f32(sum_q);
tmp[1] = vget_low_f32(sum_q);
tmp[0] = vpadd_f32(tmp[0], tmp[1]);
tmp[0] = vpadd_f32(tmp[0], tmp[0]);
result = vget_lane_f32(tmp[0], 0);

return result;
}

// 主循环直接调用函数 dot_product_intrinsic
std::vector<float> matrix_multi_result;
const int length = 2048;
for (int col = 0; col < 128; ++col)
{
float temp = 0.0;
temp = dot_product_intrinsic(fpga_output_float_value.data(), const_cast<float *> (pca_array_col[col]), length);
matrix_multi_result.push_back(temp);
}

下面是三个版本的耗时统计, 可以看到使用neon指令集优化后已经不需要fpga做这个矩阵乘法了。