并行程序设计 S5课上实验报告

并行程序设计

发布日期: 2022-10-22

文章字数: 2.8k

OpenMP 编程练习

实验内容

分别实现课件中的梯形积分法的 Pthread、OpenMP 版本，熟悉并掌握 OpenMP 编程方法，探讨两种编程方式的异同。
对于课件中“多个数组排序”的任务不均衡案例进行 OpenMP 编程实现（规模可自己调整），并探索不同循环调度方案的优劣。提示：可从任务分块的大小、线程数的多少、静态动态多线程结合等方面进行尝试，探索规律。
附加题：实现高斯消去法解线性方程组的 OpenMP 编程，与 SSE/AVX 编程结合，并探索优化任务分配方法。

梯形积分法

问题描述

对于下图中给出的梯形积分法，实现多线程编程。每个线程实现 a 到 b 区间中的某一段的梯形积分的计算。

算法设计与实现

使用Pthread实现

对于每个线程，先求局部和最后全局求和，避免过多通信。同时注意给全局和加锁保证结果的正确性。代码实现如下：

void calculate_differential(int i) {
    double x = a + i * h;
    double currArea = func(x) * h;
    pthread_mutex_lock(&barrier_mutex);
    totalPthread += currArea;
    pthread_mutex_unlock(&barrier_mutex);
}

void *Cal_pthread(void *parm) {
    int task = 0;
    while (true) {
        pthread_mutex_lock(&barrier_mutex);
        next_task += seg;
        task = next_task;
        pthread_mutex_unlock(&barrier_mutex);
        if (task >= n + seg)
            break;
        else {
            for (int i = task - seg; i < (task < n ? task : n); i++) {
                calculate_differential(i);
            }
        }
    }
    pthread_exit(nullptr);
}

使用OpenMP实现

利用 OpenMP 的相关方法实现多线程并行编程，思路与 Pthread 编程类似。代码实现如下：

void Trap(double a, double b, double h, int n, double *global_result_p) {
    double x, my_result;
    double local_a, local_b;
    int i, local_n;
    int my_rank = omp_get_thread_num();
    int thread_count = omp_get_num_threads();
    local_n = n / thread_count;
    local_a = a + my_rank * local_n * h;
    local_b = local_a + local_n * h;
    my_result = (func(local_a) + func(local_b)) / 2.0;
    for (i = 1; i <= local_n; i++) {
        x = local_a + i * h;
        my_result += func(x);
    }
    my_result = my_result * h;
# pragma omp critical
    *global_result_p += my_result;
}

实验结果分析

将待积分的函数设置为2*x^2-x，积分区间从2到10，将整个图形划分为2000个小梯形。根据实验结果，可以看出多线程和OpenMP实验结果相同，梯形面积都为613.333。

Pthread 和 OpenMP 两种编程方式异同比较

Pthread 在程序启动时创建线程，再将工作分配到线程上。然而，这种方法需要相当多的线程指定代码，而且不能保证能够随着可用处理器的数量而合理地进行扩充。OpenMP 不需要指定数量，在有循环的地方加上代码，修改设置文件即可。OpenMP 非常方便，因为它不会将软件锁定在事先设定的线程数量中，但是相对的查错更难也更麻烦。
OpenMP 和 Pthread 之间的区别主要在编译的方式上，OpenMP 的编译需要添加编译器预处理指令#pragma，创建线程等后续工作要编译器来完成。而 pthread 就是一个库，所有的并行线程创建都需要我们自己完成，较 OpenMP 麻烦一点。但如果开发人员需要精细纹理的控制，Pthread 能够提供更大范围的原函数，属于更优的选择。
OpenMP 的编译指示还有另一项重要优势：通过禁用 OpenMP 支持，代码可用作为单一线程应用进行编译。当调试程序时，以这样的方式编译代码拥有巨大优势。如果没有这种选择，开发人员会经常发现很难说明复杂的代码是否能够正确工作，因为线程问题或因为与线程无关的设计错误。

数组排序的任务分配不均衡案例复现

问题描述

对 ARR_NUM 个长度为 ARR_LEN 的一维数组进行排序，使用 OpenMP 多线程编程，每个线程处理一部分数组。若数组分布不均衡，如前二分之一的数组全部升序，后二分之一的数组全部逆序，则每个线程分得的任务负载可能也会不均衡。

算法设计与实现

均匀随机初始化数组

完全随机生成的各个数组内数据顺序完全随机，基本无差异。将数组列数设置为 2000，行数设置为 2000，线程数设置为 4，生成的数据大小范围在 0~10000 之间。

void init_uniform(void) {
    srand(unsigned(time(nullptr)));
    for (int i = 0; i < ARR_NUM; i++) {
        for (int j = 0; j < ARR_LEN; j++)
            arr[i][j] = rand() % 10000;
    }
}

不均匀随机初始化数组

不均匀随机初始化待排序数组，使得第一段数组完全逆序，第二段数组中 1/2 逆序，1/2 升序，第三段数组中 1/4 逆序，3/4 升序，第四段数组中完全升序。

// 不均匀初始化数组
void init_uneven(void){
    int ratio;
    srand(unsigned(time(nullptr)));
    for (int i = 0; i < ARR_NUM; i++){
        if(i < seg) ratio = 0;
        else if(i < seg * 2) ratio = 32;
        else if(i < seg * 3) ratio = 64;
        else ratio = 128;
        if((rand() & 127) < ratio){
            for(int j = 0; j < ARR_LEN; j++)
                arr[i][j] = ARR_LEN - j;
        } else{
            for(int j = 0; j < ARR_LEN; j++)
                arr[i][j] = j;
        }
    }
}

使用OpenMP对某块数组进行排序

指定task = next_arr += seg;在同一时间只能被一条线程执行，每次对seg块的数组进行排序。完成后线程再次领取新的任务，直至所有数组均已排序完成。

// OpenMP排序
void OpenMP_sort() {
    int task = 0;
    while (true) {
#pragma omp critical
        task = next_arr += seg;
        if (task >= ARR_NUM + seg)
            break;
        for (int i = task - seg; i < (task < ARR_NUM ? task : ARR_NUM); i++)
            sort_bubbling(arr[i]);
    }
}

控制粒度粗细

通过改变seg的大小，控制排序时每个线程分得的行数大小，从而改变粒度粗细。

    // 粗细粒度变化
    THREAD_NUM = 4;
    init_uniform();
    rollover(arr, arrtemp);
    for (seg = 10; seg <= 100; seg += 10) {
        rollover(arrtemp, arr);
        next_arr = 0;
        gettimeofday(&startTime, NULL);
#pragma omp parallel num_threads(THREAD_NUM)
        OpenMP_sort();
        gettimeofday(&stopTime, NULL);
    }

实验结果分析

均匀随机分布的数组排序

粒度变化

程序输出样例如下：

经过分析处理后得到的图表结果展示如下：

seg	10	20	30	40	50	60	70	80	90	100
time	2234.3	2212.9	2208.4	2252.2	2202.4	2341.8	2290.4	2442.8	2324.2	2193.8

从上述实验结果可以看出，在数组完全随机，数组内数据分布均匀的情况下，粒度从40增加到90时，耗时波动起伏较大；在粒度大小为80时，耗时达到峰值。同时可以看出，在粒度大小等于50时，耗时处于低谷。

证明在数组均匀随机分布的情况下，粒度大小的变化会对耗时产生较大影响，选择合适粒度能够显著提高排序速度。

线程变化

程序输出样例如下：

经过分析处理后得到的图表结果展示如下：

thread	2	4	6	8	10	12	14	16	18	20
time	4206.8	2253.7	1551.8	1216.4	1045.5	909.5	824.9	784.9	791.1	788.9

从实验结果可以看出，在数组均匀随机分布的情况下，随着线程数增加，实验耗时逐渐减小。当线程数大于12时，耗时趋于平缓。证明在线程数等于12时，排序速度已经达到较高水平，无需再增加线程数。

不均匀随机分布的数组排序

粒度变化

程序输出样例如下：

经过分析处理后得到的图表结果展示如下：

seg	10	20	30	40	50	60	70	80	90	100
time	2221.5	2235.3	2211	2277.9	2229.2	2340.2	2324.2	2382.6	2358.1	2219.2

从上述实验结果可以看出，在数组完全随机，数组内数据分布不均匀的情况下，粒度从40增加到100时，耗时波动起伏较大；在粒度大小为80时，耗时达到峰值。同时可以看出，在粒度大小等于30时，耗时处于低谷。

证明在数组不均匀随机分布的情况下，粒度大小的变化会对耗时产生较大影响，选择合适粒度能够显著提高排序速度。

线程变化

程序输出样例如下：

经过分析处理后得到的图表结果展示如下：

thread	2	4	6	8	10	12	14	16	18	20
time	1096.3	638.9	448.5	411.3	359.3	345.8	335.9	337.1	279.9	278.8

从实验结果可以看出，在数组不均匀随机分布的情况下，随着线程数增加，实验耗时逐渐减小。当线程数大于6时，耗时趋于平缓。证明在线程数等于6时，排序速度已经达到较高水平，无需再增加线程数。

高斯消去法解线性方程组的 OpenMP 多线程编程实现

问题描述

对于给定的方程组 Ax=b，使用 OpenMP多线程编程，同时结合 SSE 向量计算，完成对于线性方程组的高斯消元和回代求解。本实验采取静态划分、动态划分、粗粒度动态划分多个角度，探索任务分块大小、线程数多少、静态动态多线程结合等方面因素对高斯消去法解线性方程组效率的影响。并探索多线程分配任务的最优方案。

算法设计与实现

使用OpenMP

使用SSE，以4个元素为一个单位进行消元，剩余最后三个元素串行处理。

// s使用SSE对j行消元
void SSE_elimination(int i, int j) {
    float tep;
    __m128 div, t1, t2, sub;
    tep = a[j][i] / a[i][i];
    div = _mm_set1_ps(tep);
    // 每次处理4个元素
    int k;
    for (k = n - 3; k >= i + 1; k -= 4) {
        t1 = _mm_loadu_ps(a[i] + k);
        t2 = _mm_loadu_ps(a[j] + k);
        sub = _mm_sub_ps(t2, _mm_mul_ps(t1, div));
        _mm_storeu_ps(a[j] + k, sub);
    }
    for (k += 3; k >= i + 1; k--) {
        a[j][k] -= a[i][k] * tep;
    }
    a[j][i] = 0.00;
}

使用 OpenMP 的提供的相关方法，实现高斯消元。并结合 sse 指令向量化求解，提高效率。具体实现代码如下：

void OpenMP_func(int next_task) {
    int task = 0;
    while (true) {
#pragma omp critical
        {
            task = next_task;
            next_task += seg;
        }
        if (task >= n)
            break;
        else {
            for (int i = task; i < (task + seg < n ? task + seg : n); i++)
                SSE_elimination(line, i);
        }
    }
}