SSE编程练习
矩阵乘法的优化
问题描述
一个n行n列的矩阵a与一个n行n列的矩阵b的乘积，是一个n行n列的矩阵c。本实验分别使用了串行算法、Cache 优化、SSE 编程和分片策略四种算法实现了矩阵乘法，并就各算法的执行性能进行比较与分析。
算法设计实现与复杂性分析
串行算法
依据矩阵乘法规则，矩阵a的每一行与对应b的每一列分别相乘后求和，设计三层循环嵌套，按矩阵a的行主序处理。
具体代码如下：
void mul(float a[maxN][maxN], float b[maxN][maxN]){
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];
}
}
}
} 通过代码可以看出，串行算法的设计中有3个for循环，一个for循环是O(n)，所以串行算法的复杂度为O(n *n *n)=O(n^3)
cache优化
从上面串行算法的代码中可以看出，b [k] [j]在读取内存中的数据时是不连续的。在最底层的循环中，随着k不断加一，b[k] [j]不断地在内存中跳跃。这会引起cache命中率低，循环程序不断的把内存转移到缓存中，引起效率降低。
cache优化的思路便是通过将b矩阵转置，以行主序的方式来访问b矩阵，这样保证了b矩阵空间访问的局部性，效率得到提高。
具体代码如下：
void trans_mul(float a[maxN][maxN], float b[maxN][maxN]){
//转置
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];
}
}
}
//转置
for (int i = 0; i < n; ++i)
for (int j = 0; j < i; ++j) {
swap(b[i][j], b[j][i]);
}
} 通过代码可以看出，此算法的设计中有3个for循环，一个for循环是O(n)，所以cache优化算法的复杂度为O(n *n *n)=O(n^3)
SSE版本
矩阵乘法向量化，在cache优化的基础上，利用intrisic 中的函数，以4为步长，每次处理四个对应数据，求和后存入sum，两次相邻求和，最终sum为4维向量，均为同一个数。
具体代码如下：
void sse_mul(float a[maxN][maxN], float b[maxN][maxN]){
__m128 t1, t2, sum;
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();
for (int k = n - 4; k >= 0; k -= 4){
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);
//不能被4整除部分的处理
for (int k = (n % 4) - 1; k >= 0; --k){
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]);
}
} 此算法一次处理四个数据，算法时间复杂度为O((n^3)/4)
分片策略
将矩阵a和b分割成子矩阵，子矩阵分别计算，合并后结果与原结果相同。根据这一原理，将矩阵划分为T个矩阵块，n维矩阵元素进入Cache中的次数为n/T，即在SSE版本未分片的基础上减少循环次数.
具体代码如下：
void sse_tile(float a[maxN][maxN], float b[maxN][maxN]){
__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){
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]);
}
} 算法的时间复杂度仍然为O(n^3)，但因为充分利用了高速缓存，故计算性能也会有较显著的提升。
细节设计
计时方法
计时方面，我使用了课件中给出的Windows下精准计时方法，即QueryPerformance系列API。
以串行算法的计时为例，具体代码如下：
计时精度
本实验随机生成了30个矩阵，将计算过程重复了30次，足以匹配计算计时精度，再计算平均时间。
以串行算法为例，具体代码如下：
实验数据设计
本实验测试了不同规模的矩阵，由小至大，分别测试了规模为0、10、20、30……100的矩阵，随机生成矩阵元素值。
随机生成矩阵的过程为：
// 随机生成矩阵
void init(Matrix* dataset) {
srand(static_cast <unsigned> (time(0)));
for (int i = 0; i < 30; i++) {
for (int j = 0; j < n; j++) {
for (int k = 0; k < n; k++) {
dataset[i].a[j][k] = static_cast <float> (rand()) / (static_cast <float>(RAND_MAX / 1000));
dataset[i].b[j][k] = static_cast <float> (rand()) / (static_cast <float>(RAND_MAX / 1000));
}
}
}
}实验方法
为了降低误差，本实验重复了两次取平均值。
实验及结果分析
从图中可以看出，算法运行所需要的时间：串行算法>cache优化算法>sse优化>分片策略。
当矩阵规模较小时串行算法和cache优化算法的性能相差不大，但随着规模的增大，cache优化算法要优于串行算法。
可见，与串行算法相比，cache优化算法、sse版本以及分片策略提升了算法的效率，缩短了运行时间。
高斯消元法SSE并行化
问题描述
高斯消元法，是线性代数中的一个算法，可用来求解线性方程组。 高斯消元法的原理是：若用初等行变换将增广矩阵化为行阶梯矩阵 ，则AX = B与CX = D是同解方程组。所以我们可以用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解。
本实验分别使用了串行算法、SSE 编程分别实现了高斯消元法和回代求解的过程，并就各算法的执行性能进行比较与分析。
算法设计实现与复杂性分析
串行算法消元
按行处理，以主对角线元素为标准，主对角线右侧的值除以该行主对角线上的值，最后将对角线赋值为1。每一行处理完毕，下一行依次减去上行的值，此时值需要通过当前行的第一个不为0的元素进行还原。将对角线左侧的值均赋值为0。
具体代码如下：
void serialGauss(Matrix* dataset) {
for (int count = 0; count < 50; count++) {
for (int i = 0; i < n; i++) {
for (int j = i + 1; j < n; j++) {
float temp = dataset[count].a[j][i] / dataset[count].a[i][i];
for (int k = i; k < n; k++) {
dataset[count].a[j][k] = dataset[count].a[j][k] - temp * dataset[count].a[i][k];
}
}
}
}
} 通过代码可以看出，此算法的设计中有3个for循环，一个for循环是O(n)，所以串行算法消元的复杂度为O(n *n *n)=O(n^3)
sse消元
以4为处理单位，仍旧按行处理，每一行的数据可一次处理4个。
具体代码如下：
void sseGauss(Matrix* dataset) {
for (int count = 0; count < 50; count++) {
__m128 t1, t2, t3, t4;
for (int k = 0; k < n; k++) {
float tmp[4] = { dataset[count].a[k][k], dataset[count].a[k][k], dataset[count].a[k][k], dataset[count].a[k][k] };
t1 = _mm_loadu_ps(tmp);
for (int j = n - 4; j >= k; j -= 4) {
t2 = _mm_loadu_ps(dataset[count].a[k] + j);
t3 = _mm_div_ps(t2, t1);
_mm_storeu_ps(dataset[count].atemp[k] + j, t3);
}
if (k % 4 != (n % 4)) {
for (int j = k; j % 4 != (n % 4); j++) {
dataset[count].atemp[k][j] = dataset[count].a[k][j] / tmp[0];
}
}
for (int j = (n % 4) - 1; j >= 0; j--) {
dataset[count].atemp[k][j] = dataset[count].a[k][j] / tmp[0];
}
for (int i = k + 1; i < n; i++) {
float tmp[4] = { dataset[count].a[i][k], dataset[count].a[i][k], dataset[count].a[i][k], dataset[count].a[i][k] };
t1 = _mm_loadu_ps(tmp);
for (int j = n - 4; j > k; j -= 4) {
t2 = _mm_loadu_ps(dataset[count].a[i] + j);
t3 = _mm_loadu_ps(dataset[count].atemp[k] + j);
t4 = _mm_sub_ps(t2, _mm_mul_ps(t1, t3));
_mm_storeu_ps(dataset[count].a[i] + j, t4);
}
for (int j = k + 1; j % 4 != (n % 4); j++) {
dataset[count].a[i][j] = dataset[count].a[i][j] - dataset[count].a[i][k] * dataset[count].atemp[k][j];
}
dataset[count].a[i][k] = 0;
}
}
}
} 循环每次跳跃4个元素，算法时间复杂度为O((n^3)/4)
串行实现回代
将得到的上三角矩阵，从第 n 行开始，倒序回代到前面的每一行，即可求解 Xn 到 X1 的值。
具体代码如下：
for (int count = 0; count < 50; count++) {
for (int i = n - 1; i >= 0; i--) {
float sum = 0;
for (int j = n - 1; j > i; j--) {
sum += dataset[count].a[i][j] * dataset[count].x[j];
}
dataset[count].x[i] = (dataset[count].b[i] - sum) / dataset[count].a[i][i];
}
}
} 算法的时间复杂度为 O(n) = n^2
sse实现回代
和串行实现回代相比，向量化计算可以同时处理四个元素。
具体代码为：
void sseBack(Matrix* dataset) {
__m128 t1, t2, sumSSE;
float sum;
for (int count = 0; count < 50; count++) {
for (int i = n - 1; i >= 0; i--) {
if ((n - 1 - i) == 4)break;
float sum = 0;
for (int j = n - 1; j > i; j--) {
sum += dataset[count].a[i][j] * dataset[count].x[j];
}
dataset[count].x[i] = (dataset[count].b[i] - sum) / dataset[count].a[i][i];
}
for (int i = n - 5; i >= 0; i--) {
sumSSE = _mm_setzero_ps();
sum = 0;
int forSerial = (n - i - 1) % 4;
for (int j = i + forSerial + 1; j <= n - 4; j += 4) {
t1 = _mm_loadu_ps(dataset[count].a[i] + j);
t2 = _mm_loadu_ps(dataset[count].x + j);
t1 = _mm_mul_ps(t1, t2);
sumSSE = _mm_add_ps(sumSSE, t1);
}
sumSSE = _mm_hadd_ps(sumSSE, sumSSE);
sumSSE = _mm_hadd_ps(sumSSE, sumSSE);
_mm_store_ss(&sum, sumSSE);
for (int j = i + 1; j <= i + forSerial; j++) {
sum += dataset[count].a[i][j] * dataset[count].x[j];
}
dataset[count].x[i] = (dataset[count].b[i] - sum) / dataset[count].a[i][i];
}
}
} 算法时间复杂度 O(n) = ((n^2) / 4)
细节设计
计时方法
计时方面，我使用了课件中给出的Windows下精准计时方法，即QueryPerformance系列API。
以串行消元算法的计时为例，具体代码如下：
计时精度
本实验随机生成了50个矩阵，将计算过程重复了50次，足以匹配计算计时精度，再计算平均时间。
以串行算法为例，具体代码如下：
实验数据设计
本实验测试了不同规模的矩阵，由小至大，分别测试了规模为0、50、100、150……400的矩阵，随机生成矩阵元素值。
随机生成矩阵的过程为：
实验方法
为了降低误差，本实验重复了两次取平均值。
实验及结果分析
由图可知，SSE 算法的消元总体来说是运行效率是高于串行算法消元的;串行回代的效率也高于sse回代。
同时，通过 Excel 中趋势线的计算，可以得出消元算法的时间复杂度符合O(n^3)，验证了前面对时间复杂度的分析。同时也可以得出回代过程的时间复杂度为 O(n^2)
不同算法策略对性能的影响
相同算法对于不同问题规模的性能提升是否有影响，影响情况如何
从2.4实验结果的图表中可以看出，在问题规模较小时，串行消元和sse消元运行时间相差不大，sse消元法对于性能的提升不明显；但是当问题的规模逐渐变大sse算法明显要优于串行算法，对性能的提升更大。
回代过程可否向量化，有的话性能提升情况如何
回代过程可以向量化，根据下图可以看出回代过程采用向量化明显提升了运行效率，和串行回代相比缩短了运行时间。随着规模的增大，性能提升更显著。
