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

一. OpenMP简介

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

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

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

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

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

并行域结构:
一个能被多个线程执行的程序块,它是最基本的OpenMP并行结构,具体格式如下:
# pragma omp parallel [if (scalar_expression) | private (list) | shared (list) | defalut (shared | none) | firstprivate (list) | reduction (operator:list) | copyin (list)] newline

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支持
1

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.
2

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段;
全局交换:各处理器将其有序段按段号交换到对应的处理器中;
归并排序:各处理器对接收到的元素进行归并排序。
3

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;
}

4.运行结果:
4
显示排序结果正确,但该PSRS排序带来的额外开销使得所用时间并没有比归并排序算法快,至少在排序规模100以下,归并排序算法占优势。

已标记关键词 清除标记
©️2020 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页