并行程序设计实验

并行程序设计实验——矩阵相乘

一、问题描述

实现课件中四种矩阵乘法程序，并对比运行效率
1）串行算法
2）Cache优化
3）SSE版本
4）分片策略

二、算法设计与实现

I. 矩阵乘法算法设计

1. 用矩阵相乘的法则实现普通串行算法的设计，时间复杂度$O(n^3)$

/**
 * @param n 表示矩阵的长度
 */
void mul(int n, float a[][maxN], float b[][maxN], float c[][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];
            }
        }
    }
}

2. 根据Cache的存储和取出原理，设计Cache优化版本，时间复杂度$O(n^3)$

void trans_mul(int n, float a[][maxN], float b[][maxN], float c[][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. 利用并行的思想，在Cache优化的基础上，每4个同时相乘并且相加设计算法，时间复杂度$O(n^3)$

void sse_mul(int n, float a[][maxN], float b[][maxN], float c[][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) { // 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%4 elements
                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]);
}

4. 利用分块的思想，将矩阵每T个分为一个小矩阵块，每一个小矩阵块并行执行矩阵的乘法，时间复杂度$O(n^3)$

void sse_tile(int n, float a[][maxN], float b[][maxN], float c[][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) { //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]);
}

II. 主函数算法设计

1. 计时函数

• 计时函数采用sys/time.h库中的gettimeofday函数，该函数的精确度可以打到微妙级别，所以对函数的运行时间可以很好的估量
• 我还设置了多组循环实验，每一次循环将产生不同的矩阵，每一组矩阵乘法算法重复运行Counts次，然后得到平均估值，这样可以很好地减少一定量的误差

2. 规模自增

• 为了更好地对比四种算法的运行效率，我设计了一个GapSize值，用于矩阵规模每一次都自增GapSize，直到达到最大值。
• 另外，这里的区间个数自由度非常高，我们可以任意设置区间个数SizeCounts，以便帮助我们更好地比较算法之间的差异

const int maxN = 1024; // 矩阵的最大值
const int T = 64; // tile size
int GapSize = 128;//设置间隙大小，每一个矩阵规模自增GapSize
int SizeCounts = 10;//设置区间个数，可以自由调整
int Counts = 50;//设置每次循环的次数
float a[maxN][maxN];
float b[maxN][maxN];
float c[maxN][maxN];

//用于矩阵改变数值,为防止数据溢出,随机数的区间为1~9
void change(int n, float a[][maxN], float b[][maxN]) {
    srand(time(0));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            a[i][j] = rand() % 10;
            b[i][j] = rand() % 10;
        }
    }
}

int main(int arg, char *argv[]) {
	//设置初始时间和结束时间
    struct timeval startTime, stopTime;
    //各个矩阵规模的时间数组,为了方便画echarts表格
    double mul_times[SizeCounts], trans_mul_times[SizeCounts], sse_mul_times[SizeCounts], sse_tile_times[SizeCounts];
    for (int temp = maxN - GapSize * (SizeCounts - 1); temp <= maxN; temp += GapSize) {
        //设置每一个矩阵规模的总时间,每一个循环都加入到改变量中
        double mul_time = 0, trans_mul_time = 0, sse_mul_time = 0, sse_tile_time = 0;	
        //循环Counts次
        for (int i = 0; i < Counts; ++i) {
        	// 每一个均改变矩阵数值,控制唯一变量
            change(temp, a, b);
            gettimeofday(&startTime, NULL);
            mul(temp, a, b, c);
            gettimeofday(&stopTime, NULL);
            // 计算一个矩阵相乘的时间,并且加入的总时间中
            mul_time +=
                    (stopTime.tv_sec - startTime.tv_sec) * 1000 +
                    (double) (stopTime.tv_usec - startTime.tv_usec) * 0.001;

             gettimeofday(&startTime, NULL);
             trans_mul(temp, a, b, c);
             gettimeofday(&stopTime, NULL);
             trans_mul_time +=
                     (stopTime.tv_sec - startTime.tv_sec) * 1000 + (stopTime.tv_usec - startTime.tv_usec) * 0.001;
           
           
             gettimeofday(&startTime, NULL);
             sse_mul(temp, a, b, c);
             gettimeofday(&stopTime, NULL);
             sse_mul_time +=
                     (stopTime.tv_sec - startTime.tv_sec) * 1000 + (stopTime.tv_usec - startTime.tv_usec) * 0.001;
           
           
             gettimeofday(&startTime, NULL);
             sse_tile(temp, a, b, c);
             gettimeofday(&stopTime, NULL);
             sse_tile_time +=
                     (stopTime.tv_sec - startTime.tv_sec) * 1000 + (stopTime.tv_usec - startTime.tv_usec) * 0.001;
            
        }
        // 将每一个总和都加入到相对应的数组中
        mul_times[GapSize - 1 - (maxN - temp) / GapSize] = mul_time / Counts;
        trans_mul_times[GapSize - 1 - (maxN - temp) / GapSize] = trans_mul_time / Counts;
        sse_mul_times[GapSize - 1 - (maxN - temp) / GapSize] = sse_mul_time / Counts;
        sse_tile_times[GapSize - 1 - (maxN - temp) / GapSize] = sse_tile_time / Counts;
    }
}

三、实现结果和分析

以上是我的echarts表格,代码可见index.html.

• 我将所有实验都重复进行50次,运行时间取平均值
• 每一次矩阵的规模均匀地增长

分析

• 由图可以看出,串行的运行时间最长,且图像增长类似 $x^3$ 函数,可以看出时间复杂度为$O(n^3)$
• 由上图可以看出其余三个算法的增长都类似于 $x^3$ 函数,由此可以看出时间复杂度为$O(n^3)$
  
• 首先比较串行算法和Cache优化，可以看出，在矩阵规模达到1024的情况下，cache优化算法的运行时间基本上是串行算法的 $\frac{1}{6}$，由此可见cache优化已经大大提高了算法的运行效率，且这种优势随着矩阵规模的扩大愈发明显
  
• 比较三种优化算法，可见在实验的矩阵规模下，SSE版本的运行时间最短。可以SSE版本的效率要更胜一筹。除此之外，SSE算法是在Cache优化的基础上又进行了优化，可见此优化方法将运行效率有提高了一倍

四、源代码

源代码