一. 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支持
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;
}
4.运行结果:
显示排序结果正确,但该PSRS排序带来的额外开销使得所用时间并没有比归并排序算法快,至少在排序规模100以下,归并排序算法占优势。