OpenMP并行编程计算π值及PSRS排序

关注

原创最新推荐文章于 2025-07-23 14:54:58 发布 1.3w 阅读 AI 写同款· GEO 优化 ›

《并行计算》并行计算性能评测并行模型与并行算法呜~ 就隔了一段时间没看并行计算，发现作业贼难顶，不得不写篇博客来记录一下复习（预习）的内容。并行计算性能评测并行机的一些基本性能指标对并行计算机的性能关注点还是落在了CPU和存储器上，毕竟CPU和存储器决定了计算机处理问题速度的上限。在https://www.top500.org/收录了世界上前500台计算能力最强的计算机可以看到，超算top 5领域美国就居了3个，我国神威·太湖之光排名4，看看神威的详细配置：这么多核核内存orz… 不过好像处理器的频数并没有那么高，其实也就1.45G. 阅读详情

一. OpenMP简介

OpenMP是一个共享存储并行系统上的应用程序接口。它规范了一系列的编译制导、运行库例程和环境变量。

它提供了C/C++和FORTRAN等的应用编程接口，已经应用到UNIX、Windows NT等多种平台上。

OpenMP使用FORK-JOIN并行执行模型。所有的OpenMP程序开始于一个单独的主线程（Master Thread）。主线程会一直串行地执行，直到遇到第一个并行域（Parallel Region）才开始并行执行。

①FORK：主线程创建一队并行的线程，然后，并行域中的代码在不同的线程队中并行执行；②JOIN：当诸线程在并行域中执行完之后，它们或被同步或被中断，最后只有主线程在执行。

作用域：
编译制导语句的作用域(Scopng)可分为静态范围、孤幼制导和动态范围。

OpenMP的环境变量：
OMP_SCHEDULE：控制for循环任务分配结构的调度
OMP_NUM_THREADS：设置默认线程的个数

OpenMP的库函数
int omp_get_num_threads(void)：返回当前使用的线程个数，如果在并行区域外则返回1；
int omp_set_num_threads(int i)：设置要使用的线程个数，它可以覆盖OMP_NUM_THREADS；
int omp_get_thread_num(void)：返回当前线程号，0代表主线程；
int omp_get_num_procs(void)：返回可用的处理核（处理器）个数，对于支持超线程技术的处理器被算作两个处理核。

OpenMP的编译
windows平台 intel C++编译器：命令 icl /Qopenmp
linux平台 intel C++编译器：命令 icl -openmp
gcc：命令 gcc -fopenmp

Visual Studio 对OpenMP的支持:
属性-C/C++-语言-OpenMP支持

Linux环境对OpenMP的支持：
在Linux上编译和运行OpenMP程序
编译OpenMP程序： gcc a.c –fopenmp –o a
运行OpenMP程序： ./a

二.利用积分法计算π

在这里简单说明一下求π的积分方法，使用公式arctan(1)=π/4以及(arctan(x))’=1/(1+x^2).
在求解arctan(1)时使用矩形法求解：
求解arctan(1)是取a=0, b=1.

1.串行代码：

#include <stdio.h>
static long num_steps = 100000;//越大值越精确
double step;
void main(){
    int i;
    double x, pi, sum = 0.0;
    step = 1.0/(double)num_steps;
    for(i=1;i<= num_steps;i++){
        x = (i-0.5)*step;
        sum=sum+4.0/(1.0+x*x);
    }
    pi=step*sum;
    printf("%lf\n",pi);
}

2.使用并行域并行化的程序:

#include <stdio.h>
#include <omp.h>
static long num_steps = 100000;
double step;
#define NUM_THREADS 2
void main ()
{ 
    int i;
    double x, pi, sum[NUM_THREADS];
    step = 1.0/(double) num_steps;
    omp_set_num_threads(NUM_THREADS);  //设置2线程
 #pragma omp parallel private(i)  //并行域开始，每个线程(0和1)都会执行该代码
{
    double x;
    int id;
    id = omp_get_thread_num();
    for (i=id, sum[id]=0.0;i< num_steps; i=i+NUM_THREADS){
        x = (i+0.5)*step;
        sum[id] += 4.0/(1.0+x*x);
    }
}
    for(i=0, pi=0.0;i<NUM_THREADS;i++)  pi += sum[i] * step;
    printf("%lf\n",pi);
 }  
//共2个线程参加计算，其中线程0进行迭代步0,2,4,...线程1进行迭代步1,3,5,....

3.使用共享任务结构并行化的程序:

#include <stdio.h>
#include <omp.h>
static long num_steps = 100000;
double step;
#define NUM_THREADS 2
void main ()
{ 
    int i;
    double x, pi, sum[NUM_THREADS];
    step = 1.0/(double) num_steps;
    omp_set_num_threads(NUM_THREADS);  //设置2线程
 #pragma omp parallel  //并行域开始，每个线程(0和1)都会执行该代码
{
    double x;
    int id;
    id = omp_get_thread_num();
    sum[id]=0;
#pragma omp for  //未指定chunk，迭代平均分配给各线程（0和1），连续划分
    for (i=0;i< num_steps; i++){
        x = (i+0.5)*step;
        sum[id] += 4.0/(1.0+x*x);
    }
}
    for(i=0, pi=0.0;i<NUM_THREADS;i++)  pi += sum[i] * step;
    printf("%lf\n",pi);
}//共2个线程参加计算，其中线程0进行迭代步0~49999，线程1进行迭代步50000~99999.

4.使用private子句和critical部分并行化的程序:

#include <stdio.h>
#include <omp.h>
static long num_steps = 100000;
double step;
#define NUM_THREADS 2
void main ()
{ 
    int i;
    double pi=0.0;
    double sum=0.0;
    double x=0.0;
    step = 1.0/(double) num_steps;
    omp_set_num_threads(NUM_THREADS);  //设置2线程
#pragma omp parallel private(x,sum,id) //该子句表示x,sum变量对于每个线程是私有的(!!!注意id要设为私有变量)
{
    int id;
    id = omp_get_thread_num();
    for (i=id, sum=0.0;i< num_steps; i=i+NUM_THREADS){
        x = (i+0.5)*step;
        sum += 4.0/(1.0+x*x);
    }
#pragma omp critical  //指定代码段在同一时刻只能由一个线程进行执行
    pi += sum*step;
}
    printf("%lf\n",pi);
}   //共2个线程参加计算，其中线程0进行迭代步0,2,4,...线程1进行迭代步1,3,5,....当被指定为critical的代码段  正在被0线程执行时，1线程的执行也到达该代码段，则它将被阻塞知道0线程退出临界区。

5.使用并行规约的并行程序：

#include <stdio.h>
#include <omp.h>
static long num_steps = 100000;
double step;
#define NUM_THREADS 2
void main ()
{ 
    int i;
    double pi=0.0;
    double sum=0.0;
    double x=0.0;
    step = 1.0/(double) num_steps;
    omp_set_num_threads(NUM_THREADS);  //设置2线程
#pragma omp parallel for reduction(+:sum) private(x) //每个线程保留一份私有拷贝sum，x为线程私有，最后对线程中所以sum进行+规约，并更新sum的全局值
    for(i=1;i<= num_steps; i++){
        x = (i-0.5)*step;
        sum += 4.0/(1.0+x*x);
    }
    pi = sum * step;
    printf("%lf\n",pi);
}   //共2个线程参加计算，其中线程0进行迭代步0~49999，线程1进行迭代步50000~99999.

三.PSRS算法的设计

1.PSRS算法简介
PSRS排序可分为8个部分：
均匀划分：将n个元素A[1..n]均匀划分成p段，每个pi处理A[(i-1)n/p+1..in/p]；
局部排序：pi调用串行排序算法对A[(i-1)n/p+1..in/p]排序；
正则采样：pi从其有序子序列A[(i-1)n/p+1..in/p]中选取p个样本元素；
采样排序：用一台处理器对p2个样本元素进行串行排序；
选择主元：用一台处理器从排好序的样本序列中选取p-1个主元，并
播送给其他pi；
主元划分：pi按主元将有序段A[(i-1)n/p+1..in/p]划分成p段；
全局交换：各处理器将其有序段按段号交换到对应的处理器中；
归并排序：各处理器对接收到的元素进行归并排序。

2.PSRS并行算法设计思想

首先确定算法中调用的串行排序为归并排序，总共需要调用3次。
接下来对算法的各个部分进行并行设计，伪代码如下：

定义处理器处理段个数数组变量、划分段的大小、正则采样数p、主元数组、处理器各部分排序的三维数组空间。
设置并行线程后，添加第一个并行域：
\#pragma omp parallel shared(base,array,n,i,pivot,count) private(id)
{
　　每个处理器对所在的段进行局部串行归并排序；
　　并行域嵌套
　　#pragma omp critical
　　{
　　　　每个处理器选出p个样本；
   　　 }
　　设置路障#pragma omp barrier
　　主线程并行域#pragma omp master
　　{
   　　　　 选出num_threads-1个主元
　　}
　　设置路障#pragma omp barrier
　　{
    　　　　根据主元对每一个cpu数据进行划分
　　}
}
第二个并行域：
＃pragma omp parallel shared(pivot_array,count)
{
　　向各个线程发送数据，各个线程自己排序；
　　打印输出数据；
}
主函数测试数组排序。

3.PSRS排序并行算法实现：

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<omp.h>

#define num_threads 3

int *L,*R;

//Merge函数合并两个子数组形成单一的已排好序的字数组
//并代替当前的子数组A[p..r] 
void Merge(int *a, int p, int q, int r)
{
    int i,j,k;
    int n1 = q - p + 1;
    int n2 = r - q;
    L = (int*)malloc((n1+1)*sizeof(int));
    R = (int*)malloc((n2+1)*sizeof(int));
    for(i=0; i<n1; i++)
        L[i] = a[p+i];
    L[i] = 65536;
    for(j=0; j<n2; j++)
        R[j] = a[q+j+1];
    R[j] = 65536;
    i=0,j=0;
    for(k=p; k<=r; k++){
        if(L[i]<=R[j]){
            a[k] = L[i];
            i++;
        }
        else{
            a[k] = R[j];
            j++;
        }
    }
} 
//归并排序
void MergeSort(int *a, int p, int r)
{
    if(p<r){
        int q = (p+r)/2;
        MergeSort(a,p,q);
        MergeSort(a,q+1,r);
        Merge(a,p,q,r);
    }
} 

void PSRS(int *array, int n)
{
    int id;
    int i=0;
    int count[num_threads][num_threads] = { 0 };    //每个处理器每段的个数
    int base = n / num_threads;     //划分的每段段数
    int p[num_threads*num_threads]; //正则采样数为p
    int pivot[num_threads-1];       //主元
    int pivot_array[num_threads][num_threads][50]={0};  //处理器数组空间

    omp_set_num_threads(num_threads);
#pragma omp parallel shared(base,array,n,i,p,pivot,count) private(id)
    {
        id = omp_get_thread_num();
        //每个处理器对所在的段进行局部串行归并排序
        MergeSort(array,id*base,(id+1)*base-1);

#pragma omp critical
        //每个处理器选出P个样本，进行正则采样
        for(int k=0; k<num_threads; k++)
            p[i++] = array[(id-1)*base+(k+1)*base/(num_threads+1)];
//设置路障，同步队列中的所有线程
#pragma omp barrier
//主线程对采样的p个样本进行排序
#pragma omp master
        {
            MergeSort(p,0,i-1);
            //选出num_threads-1个主元
            for(int m=0; m<num_threads-1; m++)
                pivot[m] = p[(m+1)*num_threads];
        }
#pragma omp barrier
        //根据主元对每一个cpu数据段进行划分
        for(int k=0; k<base; k++)
        {
            for(int m=0; m<num_threads; m++)
            {
                if(array[id*base+k] < pivot[m])
                {
                    pivot_array[id][m][count[id][m]++] = array[id*base+k];
                    break;
                }
                else if(m == num_threads-1) //最后一段
                {
                    pivot_array[id][m][count[id][m]++] = array[id*base+k];
                }
            }
        }
    }
//向各个线程发送数据，各个线程自己排序
#pragma omp parallel shared(pivot_array,count)
    {
        int id=omp_get_thread_num();
        for(int k=0; k<num_threads; k++)
        {
            if(k!=id)
            {
                memcpy(pivot_array[id][id]+count[id][id],pivot_array[k][id],sizeof(int)*count[k][id]);
                count[id][id] += count[k][id];
            }
        }
        MergeSort(pivot_array[id][id],0,count[id][id]-1);
    }
    i = 0;
    printf("result:\n");
    for(int k=0; k<num_threads; k++)
    {
        for(int m=0; m<count[k][k]; m++)
        {
            printf("%d ",pivot_array[k][k][m]);
        }
        printf("\n");
    }
}

int main()
{
    int array[36] = { 16,2,17,24,33,28,30,1,0,27,9,25, 34,23,19,18,11,7,21,13,8,35,12,29 , 6,3,4,14,22,15,32,10,26,31,20,5 };
    double begin,end,time;
    begin = omp_get_wtime();
    PSRS(array, 36);
/*  MergeSort(list,0,35);
    for(int i=0; i<36; i++)
    {
        printf("%d ",list[i]);
    }*/
    end = omp_get_wtime();
    time = end - begin;
    printf("The running time is %lfs\n",time);
    return 0;
}