1 问题描述
本次实验分别使用串行算法、Cache优化算法、SSE编程和分片策略算法实现了矩阵乘法运算,实验采用同一个样本,即矩阵大小为512个元素,元素值为由时间生成的随机数,每个算法对此样本运行十次,并记录每次运行时间和十次运算的平均运行时间。实验环境:计算机apple macbook pro2015、系统macOS High Sierra10.13.5、编辑器vscode&C/C++ Extension、gcc编译。
2 算法设计与实现
首先在矩阵样本设计方面,考虑大越大规模的数据越能体现不同算法在运行上的差距,因而将矩阵大小设置为512个元素,并且为了减少元素数据对实验的影响,采取了随机数赋值的方法。
srand((unsigned)time(NULL));
float a[maxN][maxN];
float b[maxN][maxN];
for(int i = 0; i < n; i++)
{
for(int j=0; j < n; j++){
a[i][j] = rand();
b[i][j] = rand();
}
}
其次考虑每个算法的具体实现:
2.1 串行算法
由于矩阵乘法规则为乘积AB的第i个分量等于矩阵A的第i个行向量与列向量B的对应分量乘积之和,因而对矩阵A的每一行的元素,矩阵B的每一列的元素都要与之相乘并相加,所以需要先遍历矩阵A一行中的元素的同时遍历矩阵B一列中的元素,再遍历矩阵B中的每一列,再遍历矩阵A中的每一行,三层遍历嵌套,时间复杂度O(n) = n^3,用clock_t记录时间。
具体代码如下:
double mul(int n,float a[][maxN],float b[][maxN],float c[][maxN]){ //串行算法 O(n) = n^3
double time;
clock_t start,end;
start = clock();
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
{
c[i][j]= 0.0;
for(int k = 0; k < n; k++){
c[i][j] += a[i][k] * b[k][j];
}
}
}
end = clock();
time = (double)(end-start);
cout<<"串行算法耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
2.2 Cache优化
由于寄存器从内存中取值时,需要从内存中寻址,上述串行算法b[k][j]
需要先对矩阵B的行进行遍历,又因为在c++中,二维数组相邻的元素地址相近,因而先对行进行遍历造成了地址跳跃,需要更长的寻址时间。将矩阵B转置,就能够实现对连续的地址进行读写操作,节省大量寻址时间。算法时间复杂度O(n) = n^3。
具体代码如下:
double trans_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]){ // Cache优化 O(n) = n^3
double time;
clock_t start,end;
start = clock();
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
c[i][j] = 0.0;
for (int k = 0; k < n; ++k) {
c[i][j] += a[i][k] * b[j][k];
}
}
} // transposition
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
end = clock();
time = (double)(end-start);
cout<<"Cache优化耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
2.3 SSE版本
基于SIMD编程,将多个元素同时进行寄存器加载和存储器存值等操作能够提高并行效率。使用SSE指令能够实现该功能。时间复杂度为O(n) = (n^3)/4,同样对矩阵B进行转置。
具体代码如下:
double sse_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]){ //SSE版本 O(n) = (n^3)/4
double time;
clock_t start,end;
start = clock();
__m128 t1, t2, sum; //__m128 == float
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
for (int i = 0; i < n; ++i){
for (int j = 0; j < n; ++j){
c[i][j] = 0.0;
sum = _mm_setzero_ps(); //Initialize
for (int k = n - 4; k >= 0; k -= 4){ // sum every 4 elements
t1 = _mm_loadu_ps(a[i] + k);
t2 = _mm_loadu_ps(b[j] + k);
t1 = _mm_mul_ps(t1, t2);
sum = _mm_add_ps(sum, t1);
}
sum = _mm_hadd_ps(sum, sum);
sum = _mm_hadd_ps(sum, sum);
_mm_store_ss(c[i] + j, sum);
for (int k = (n % 4) - 1; k >= 0; --k){ // handle the last n%4elements
c[i][j] += a[i][k] * b[j][k];
}
}
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
end = clock();
time = (double)(end-start);
cout<<"SSE版本耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
其中_mm_loadu_ps为同时加载四个单精度浮点数值的指令,因而循环每次跳跃四个元素,再将四个元素的乘积相加。对于矩阵n不能被4整除的情况,对剩余元素做类似Cache优化算法的操作。
2.4 分片策略
先针对cpu线程数设置分片数,将矩阵分为多个小矩阵,对每个小矩阵进行SSE指令计算。时间复杂度为O(n) = n^3,同样对矩阵B进行转置。
具体代码如下:
double sse_tile(int n, float a[][maxN], float b[][maxN], float c[][maxN]){ //分片策略 O(n) = n^3
double time;
clock_t start,end;
start = clock();
__m128 t1, t2, sum;
float t;
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
for (int r = 0; r < n / T; ++r)
for (int q = 0; q < n / T; ++q) {
for (int i = 0; i < T; ++i)
for (int j = 0; j < T; ++j)
c[r * T + i][q * T +j] = 0.0;
for (int p = 0; p < n / T; ++p) {
for (int i = 0; i < T; ++i)
for (int j = 0; j < T; ++j) {
sum = _mm_setzero_ps();
for (int k = 0; k < T; k += 4){ //sum every 4th elements
t1 = _mm_loadu_ps(a[r * T + i] + p * T + k);
t2 = _mm_loadu_ps(b[q * T + j] + p * T + k);
t1 = _mm_mul_ps(t1, t2);
sum = _mm_add_ps(sum, t1);
}
sum = _mm_hadd_ps(sum, sum);
sum = _mm_hadd_ps(sum, sum);
_mm_store_ss(&t, sum);
c[r * T + i][q * T + j] += t;
}
}
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
end = clock();
time = (double)(end-start);
cout<<"分片策略耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
3 实验结果及结果分析
每个算法的十次运行时间结果如下:
第1次运行结果:
串行算法耗时: 0.753588s
Cache优化耗时: 0.453238s
SSE版本耗时: 0.230628s
分片策略耗时: 0.279179s
第2次运行结果:
串行算法耗时: 0.77117s
Cache优化耗时: 0.465396s
SSE版本耗时: 0.239504s
分片策略耗时: 0.293299s
第3次运行结果:
串行算法耗时: 0.746869s
Cache优化耗时: 0.451797s
SSE版本耗时: 0.228981s
分片策略耗时: 0.282906s
第4次运行结果:
串行算法耗时: 0.696163s
Cache优化耗时: 0.453274s
SSE版本耗时: 0.233736s
分片策略耗时: 0.277309s
第5次运行结果:
串行算法耗时: 0.878417s
Cache优化耗时: 0.490044s
SSE版本耗时: 0.271495s
分片策略耗时: 0.378412s
第6次运行结果:
串行算法耗时: 1.08664s
Cache优化耗时: 0.491328s
SSE版本耗时: 0.243387s
分片策略耗时: 0.291609s
第7次运行结果:
串行算法耗时: 0.777393s
Cache优化耗时: 0.470352s
SSE版本耗时: 0.233467s
分片策略耗时: 0.285942s
第8次运行结果:
串行算法耗时: 0.689719s
Cache优化耗时: 0.450596s
SSE版本耗时: 0.224129s
分片策略耗时: 0.284757s
第9次运行结果:
串行算法耗时: 0.692528s
Cache优化耗时: 0.450607s
SSE版本耗时: 0.229762s
分片策略耗时: 0.282219s
第10次运行结果:
串行算法耗时: 0.708097s
Cache优化耗时: 0.45446s
SSE版本耗时: 0.223255s
分片策略耗时: 0.28649s
平均耗时:
串行算法耗时:0.780058s
Cache优化耗时:0.463109s
SSE版本耗时:0.235834s
分片策略耗时:0.294212s
图表展示如下:
算法 | 1 | 2 | 3 | 4 | 5 |
---|---|---|---|---|---|
串行算法(s) | 0.75 | 0.77 | 0.75 | 0.69 | 0.88 |
Cache优化(s) | 0.45 | 0.47 | 0.45 | 0.45 | 0.49 |
SSE版本(s) | 0.23 | 0.24 | 0.23 | 0.23 | 0.27 |
分片策略(s) | 0.28 | 0.29 | 0.28 | 0.28 | 0.38 |
算法 | 6 | 7 | 8 | 9 | 10 |
---|---|---|---|---|---|
串行算法(s) | 1.09 | 0.77 | 0.69 | 0.69 | 0.71 |
Cache优化(s) | 0.49 | 0.47 | 0.45 | 0.45 | 0.45 |
SSE版本(s) | 0.24 | 0.23 | 0.22 | 0.23 | 0.22 |
分片策略(s) | 0.29 | 0.29 | 0.28 | 0.28 | 0.29 |
从图中可以明显看出,耗时最多的算法时串行算法,也是最原始的算法,Cache优化对串行算法有很明显的优化效果,时间减短大约40% ,但用时仍比SSE编程要长,SSE版本算法运行时间最短,采用分片策略后时间反而有点增长。
由于串行算法没做任何优化,因而运行时间最长。Cache优化算法对串行算法的取址时间进行了压缩。而想要更短的运行时间就需要用到SIMD编程。分片策略比SSE版本用时较长的原因大概是分片数量不理想,并行成本消耗过高。
4 源码
本次实验代码如下:
#include<iostream>
#include<stdlib.h>
#include <pmmintrin.h>
#include <stdio.h>
#include <algorithm>
#include<cmath>
using namespace std;
const int maxN = 512;
const int T = 64; // tile size
double mul(int n,float a[][maxN],float b[][maxN],float c[][maxN]);
double trans_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]);
double sse_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]);
double sse_tile(int n, float a[][maxN], float b[][maxN], float c[][maxN]);
int main(){
int n = 512;
srand((unsigned)time(NULL));
float a[maxN][maxN];
float b[maxN][maxN];
float c[maxN][maxN];
double mul_sum = 0;
double trans_mul_sum = 0;
double sse_mul_sum = 0;
double sse_tile_sum = 0;
for(int i = 0; i < n; i++)
{
for(int j=0; j < n; j++){
a[i][j] = rand();
b[i][j] = rand();
}
}
for(int i=0;i<10;i++){
cout<<"第"<<i+1<<"次运行结果:"<<endl;
mul_sum += mul(n,a,b,c);
trans_mul_sum += trans_mul(n,a,b,c);
sse_mul_sum += sse_mul(n,a,b,c);
sse_tile_sum += sse_tile(n,a,b,c);
cout<<endl;
}
cout<<"平均耗时:"<<endl;
cout<<"串行算法耗时:"<<mul_sum/10<<"s"<<endl;
cout<<"Cache优化耗时:"<<trans_mul_sum/10<<"s"<<endl;
cout<<"SSE版本耗时:"<<sse_mul_sum/10<<"s"<<endl;
cout<<"分片策略耗时:"<<sse_tile_sum/10<<"s"<<endl;
}
double mul(int n,float a[][maxN],float b[][maxN],float c[][maxN]){ //串行算法 O(n) = n^3
double time;
clock_t start,end;
start = clock();
for(int i = 0; i < n; i++)
{
for(int j = 0; j < n; j++)
{
c[i][j]= 0.0;
for(int k = 0; k < n; k++){
c[i][j] += a[i][k] * b[k][j];
}
}
}
end = clock();
time = (double)(end-start);
cout<<"串行算法耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
double trans_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]){ // Cache优化 O(n) = n^3
double time;
clock_t start,end;
start = clock();
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
c[i][j] = 0.0;
for (int k = 0; k < n; ++k) {
c[i][j] += a[i][k] * b[j][k];
}
}
} // transposition
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
end = clock();
time = (double)(end-start);
cout<<"Cache优化耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
double sse_mul(int n, float a[][maxN], float b[][maxN], float c[][maxN]){ //SSE版本 O(n) = (n^3)/4
double time;
clock_t start,end;
start = clock();
__m128 t1, t2, sum; //__m128 == float
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
for (int i = 0; i < n; ++i){
for (int j = 0; j < n; ++j){
c[i][j] = 0.0;
sum = _mm_setzero_ps(); //Initialize
for (int k = n - 4; k >= 0; k -= 4){ // sum every 4 elements
t1 = _mm_loadu_ps(a[i] + k);
t2 = _mm_loadu_ps(b[j] + k);
t1 = _mm_mul_ps(t1, t2);
sum = _mm_add_ps(sum, t1);
}
sum = _mm_hadd_ps(sum, sum);
sum = _mm_hadd_ps(sum, sum);
_mm_store_ss(c[i] + j, sum);
for (int k = (n % 4) - 1; k >= 0; --k){ // handle the last n%4elements
c[i][j] += a[i][k] * b[j][k];
}
}
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
end = clock();
time = (double)(end-start);
cout<<"SSE版本耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}
double sse_tile(int n, float a[][maxN], float b[][maxN], float c[][maxN]){ //分片策略 O(n) = n^3
double time;
clock_t start,end;
start = clock();
__m128 t1, t2, sum;
float t;
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
for (int r = 0; r < n / T; ++r)
for (int q = 0; q < n / T; ++q) {
for (int i = 0; i < T; ++i)
for (int j = 0; j < T; ++j)
c[r * T + i][q * T +j] = 0.0;
for (int p = 0; p < n / T; ++p) {
for (int i = 0; i < T; ++i)
for (int j = 0; j < T; ++j) {
sum = _mm_setzero_ps();
for (int k = 0; k < T; k += 4){ //sum every 4th elements
t1 = _mm_loadu_ps(a[r * T + i] + p * T + k);
t2 = _mm_loadu_ps(b[q * T + j] + p * T + k);
t1 = _mm_mul_ps(t1, t2);
sum = _mm_add_ps(sum, t1);
}
sum = _mm_hadd_ps(sum, sum);
sum = _mm_hadd_ps(sum, sum);
_mm_store_ss(&t, sum);
c[r * T + i][q * T + j] += t;
}
}
}
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j)
swap(b[i][j], b[j][i]);
end = clock();
time = (double)(end-start);
cout<<"分片策略耗时: "<<(time/CLOCKS_PER_SEC)<<"s"<<endl;
return time/CLOCKS_PER_SEC;
}