利用原子指令加速矩阵运算
C++里面有个原子指令库,不需要通过内嵌汇编就可以调用cpu内部SIMD的指令,头文件<xmmintrin.h> 包含SSE库,<emmintrin.h> 包含SSE2库,<pmmintrin.h> 包含SSE3库。
原子指令可以利用cpu内部的128位寄存器,同时进行多个数据运算,比如可以同时计算4个float类型数,或2个double类型的数。
运行环境:
win10 下ubuntu18.04
g++ 编译的时候要加关键选项 -msse3 否则识别不了_mm_hadd_pd。非常关键,百度找不到,最终还是借助bing英文搜索。
g++ -msse3 test.cpp -o test
# include <iostream>
# include <cstdlib>
# include <ctime>
# include <pmmintrin.h>
using namespace std;
//show result
template<typename T>
void show_matrix(int size,T ** M){
for(int i=0;i<size;i++){
for(int j = 0;j<size;j++){
cout<<M[i][j]<<" ";
}
cout<<"\n";
}
}
// roll over matrix b
template<typename T>
void Transpose(int size,T** m)
{
for(int i=0;i<size;++i){
for(int j=i+1;j<size;++j){
std::swap(m[i][j],m[j][i]);
}
}
}
//method 1
template<typename T>
void SeqMatrixMult1(int size, T** m1, T** m2, T** result)
{
Transpose(size, m2);
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
// temp parameter can reduce memory access, which is very important
T c = 0;
for (int k = 0; k < size; k++) {
c += m1[i][k] * m2[j][k];
}
result[i][j] = c;
}
}
Transpose(size, m2);
}
//method 2
//template<typename T>
void SeqMatrixMult2(int size, double** m1, double** m2, double** result)
{
Transpose(size, m2);
for (int i = 0; i < size; i++) {
for (int j = 0; j < size; j++) {
__m128d c = _mm_setzero_pd();
for (int k = 0; k < size; k += 2) {
c = _mm_add_pd(c, _mm_mul_pd(_mm_load_pd(&m1[i][k]), _mm_load_pd(&m2[j][k])));
}
// horizontal add of the single register
c = _mm_hadd_pd(c, c);
_mm_store_sd(&result[i][j], c);
}
}
Transpose(size, m2);
}
int main(){
int n = 500;
double **a = new double*[n];
double **b = new double*[n];
double **result1 = new double*[n];
double **result2 = new double*[n];
for(int i=0;i<n;i++){
a [i]=new double[n];
b [i]=new double[n];
result1[i]=new double[n];
result2[i]=new double[n];
}
for(int i=0;i<n;++i)
for(int j=0;j<n;++j){
a[i][j]=rand()%100/10.0;
b[i][j]=rand()%100/10.0;
result1[i][j]=0;
result2[i][j]=0;
}
// show_matrix(n,a);
// show_matrix(n,b);
clock_t start,end;
// method 1
start=clock();
SeqMatrixMult1(n,a,b,result1);
end = clock();
cout<<end-start<<"\n";
// method 2
start=clock();
SeqMatrixMult2(n,a,b,result2);
end = clock();
cout<<end-start<<"\n";
// show_matrix(n,result1);
// show_matrix(n,result2);
//delete;
for(int i=0;i<n;++i){
delete[]a[i];
delete[]b[i];
delete[]result1[i];
delete[]result2[i];
}
delete[]a;
delete[]b;
delete[]result1;
delete[]result2;
return 0;
}
参考:
Optimize Your Code: Matrix Multiplication
'mm_hadd_ps' was not declared in this scope
网友评论