OpenMP编程实例（向量点积、矩阵乘法和快排）

之前的一篇文章给出了用OpenMP编程求PI的四种方法，算是OpenMP编程的初级实例，本文将继续探讨用OpenMP并行编程解决更多的问题。

关于求PI的实现方法可见我的前一篇博客：/article/11836426.html

1.使用临界区，并行规约方式来实现两个向量点积的运算

首先是临界区的方法：

#include <stdio.h>
#include <omp.h>
int main(){
int i,n;
float a[100],b[100],result;
//some initializations
n=100;
result=0.0;
for(i=0;i<n;i++){
a[i]=i*1.0;
b[i]=i*2.0;
}
#pragma omp parallel for
for(i=0;i<n;i++){
#pragma omp critical
result=result+(a[i]*b[i]);
}

printf("%f\n",result);
}

这里很简单的给两个数组设置了初始值，也可以用随机数发生器来产生两个数组。

或者采用并行规约的方式：

#include <stdio.h>
#include <omp.h>
//reduction
int main(){
int i,n,chunk;
float a[100],b[100],result;
//some initializations
n=100;
chunk=10;
result=0.0;
for(i=0;i<n;i++){
a[i]=i*1.0;
b[i]=i*2.0;
}
#pragma omp parallel for default(shared)private(i)schedule(static,chunk)reduction(+:result)

for(i=0;i<n;i++)
result=result+(a[i]*b[i]);

printf("%f\n",result);
}

可以在两段代码里加入计时函数，然后比较两种方法的性能。

2.矩阵乘法

同样是简单赋值初始化两个N×N的矩阵，OpenMP实现如下：

#include <stdio.h>
#include <omp.h>
#define N 10
#define n 2
int main()
{
int a

;
int b

;
int c

;
int i,j,k;
//some initializations
for(i=0;i<N;i++){
for(j=0;j<N;j++){
a[i][j]=1;
b[i][j]=1;
c[i][j]=0;
}
}
omp_set_num_threads(n);
#pragma omp parallel shared(a,b,c) private(i,j,k)
{
#pragma omp for schedule(dynamic)
for(i=0;i<N;i++){
for(j=0;j<N;j++){
for(k=0;k<N;k++)
c[i][j]+=a[i][k]*b[k][j];
}
}
}

for(i=0;i<N;i++){
for(j=0;j<N;j++){
printf("%d ",c[i][j]);
}
printf("\n");
}
}

这里是采用打印的方式查看结果，可同样加入计时函数，分析加速比，分析时间与线程数的关系。

3.用OpenMP实现快速排序算法

回顾一下快排的串行算法，主要是以下一个分治的过程：

QUICKSORT(A,p,r)
if p<r
then q=PARTITION(A,p,r)
QUICKSORT(A,p,q-1)
QUICKSORT(A,q+1,r)

快排的关键是PATITION过程，它对子数组A[p..r]进行就地重排：

PARTITION(A,p,r)
x=A[r]
i=p-1
for j=p to r-1
do if A[j]<=x
then i=i+1
exchange A[i] and A[j]
exchange A[i+1] and A[r]
return i+1

基于这样的串行算法，我们可以给出快排的OpenMP实现：

#include<stdio.h>
#include<stdlib.h>
#include<omp.h>
#include<time.h>
#define NUM_THREADS 2

int Partition(int *data,int p,int r)
{
int x,i,j,tmp;
//即A[r]
x=data[r];
i=p-1;
for(j=p;j<r;j++)
if(data[j]<=x)
{
i++;
//exchange A[i] and A[j]
tmp=data[i];
data[i]=data[j];
data[j]=tmp;
}
//exchange A[i+1] and A[r]
tmp=data[i+1];
data[i+1]=data[r];
data[r]=tmp;

return i+1;
}
int* QuickSort_parallel(int* Data,int p,int r)
{
int q;
if(p<r)
{
q=Partition(Data,p,r);
#pragma omp parallel
{
#pragma omp sections nowait
{
#pragma omp section
QuickSort_parallel(Data,p,q-1);
#pragma omp section
QuickSort_parallel(Data,q+1,r);
}
}
}
return Data;
}
void main()
{
int N;
printf("输入产生多少个数进行排序：");
scanf("%d",&N);
int *Data;
Data=new int
;
for(int i=0;i<N;i++)
//0~32765
Data[i]=rand();
clock_t timep =clock();
omp_set_num_threads(NUM_THREADS);
Data=QuickSort_parallel(Data,0,N-1);
clock_t timeFinish =clock();
for (int i=0;i<N;i++)
{
printf("%d\n",Data[i]);
}
printf("并行快速排序所花时间：%d ms\n",timeFinish-timep);
}

4.总结

本文以向量点积运算、矩阵乘法和快速排序为例，进一步对OpenMP编程进行了实战，可以发现，OpenMP编程还是较容易上手和掌握的。