按行分配
思路和MPI基本类似，不过OpenMP是共享内存的，不必做分发和聚集，申请的矩阵空间就不必是完全连续的。
1 #include<stdio.h>
2 #include<omp.h>
3 #include<stdlib.h>
4
5 #define N 400 //规模(方针的阶数)
6 int i,j;//通用游标
7 double **mat=NULL;//矩阵对象
8 int *vec=NULL;//列向量对象
9 double *result=NULL;//结果向量对象
10
11 int main(int argc,char* argv[])
12 {
13 //读入10进制表示的进程数
14 int thrdCnt=strtol(argv[1],NULL,10);
15 //矩阵一级挂载空间
16 mat=(double**)malloc(N*sizeof(double *));
17 //存向量的空间
18 vec=(int *)malloc(N*sizeof(int));
19 //存结果的空间
20 result=(double *)malloc(N*sizeof(double));
21
22 //并行for:申请存矩阵的空间
23 //因为OpenMP不需要分发聚集等,所以不必做连续空间申请
24 # pragma omp parallel for num_threads(thrdCnt)
25 for(i=0;i<N;i++)
26 mat[i]=(double*)malloc(N*sizeof(double));
27
28 //并行for:为矩阵和列向量赋值(实际做时是从文件读入)
29 # pragma omp parallel for num_threads(thrdCnt) private(j)
30 for(i=0;i<N;i++)
31 {
32 vec[i]=1;
33 for(j=0;j<N;j++)
34 mat[i][j]=j;
35 }
36
37 //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断
38 # pragma omp parallel for num_threads(thrdCnt)
39 for(i=0;i<N;i++)
40 result[i]=mat[i][0]*vec[0];
41
42 //并行for:计算结果
43 # pragma omp parallel for num_threads(thrdCnt) private(j)
44 for(i=0;i<N;i++)
45 for(j=1;j<N;j++)
46 result[i]+=mat[i][j]*vec[j];
47
48 //采样输出结果看一下
49 for(i=0;i<N;i+=N/11)
50 printf("%f\n",result[i]);
51
52 return 0;
53 }
输出
1 [lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_row.o omp_row.c
2 [lzh@hostlzh OpenMP]$ ./omp_row.o 7
3 79800.000000
4 79800.000000
5 79800.000000
6 79800.000000
7 79800.000000
8 79800.000000
9 79800.000000
10 79800.000000
11 79800.000000
12 79800.000000
13 79800.000000
14 79800.000000
15 [lzh@hostlzh OpenMP]$
按列分配
按列分配有很多细节要注意，如果不对result数组保护则可能会发生冲突，如果用critical或者atomic会导致计算过程实际是串行的，虽然atomic的加解锁过程是critical的7倍。
1 #include<stdio.h>
2 #include<omp.h>
3 #include<stdlib.h>
4
5 #define N 400 //规模(方针的阶数)
6 int i,j;//通用游标
7 double **mat=NULL;//矩阵对象
8 int *vec=NULL;//列向量对象
9 double *result=NULL;//结果向量对象
10
11 int main(int argc,char* argv[])
12 {
13 //读入10进制表示的进程数
14 int thrdCnt=strtol(argv[1],NULL,10);
15 //矩阵一级挂载空间
16 mat=(double**)malloc(N*sizeof(double *));
17 //存向量的空间
18 vec=(int *)malloc(N*sizeof(int));
19 //存结果的空间
20 result=(double *)malloc(N*sizeof(double));
21
22 //并行for:申请存矩阵的空间
23 //因为OpenMP不需要分发聚集等,所以不必做连续空间申请
24 # pragma omp parallel for num_threads(thrdCnt)
25 for(i=0;i<N;i++)
26 mat[i]=(double*)malloc(N*sizeof(double));
27
28 //并行for:为矩阵和列向量赋值(实际做时是从文件读入)
29 # pragma omp parallel for num_threads(thrdCnt) private(j)
30 for(i=0;i<N;i++)
31 {
32 vec[i]=1;
33 for(j=0;j<N;j++)
34 mat[i][j]=j;
35 }
36
37 //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断
38 # pragma omp parallel for num_threads(thrdCnt)
39 for(i=0;i<N;i++)
40 result[i]=mat[i][0]*vec[0];
41
42 //创建存放每个线程临时结果的数组,初始化为0
43 //这里calloc还能像前面那样改进,但为了程序可读性暂时不做改进
44 double* sum;
45
46 //OpenMP默认块划分给每个进程(除了最后一个进程)的循环次数
47 //是总的循环次数处以进程数向上取整,需要先计算出这个数
48 int needPlus;//记录是否向上+1
49 needPlus=((N-1)%thrdCnt==0)?0:1;
50 //每个进程(除了最后一个进程)的循环次数
51 int numThrdFor=(N-1)/thrdCnt+needPlus;
52
53 //并行for:计算结果,按列分配时这个大的外层for用j,且从1开始
54 # pragma omp parallel for num_threads(thrdCnt) \
55 private(i) firstprivate(sum)//对i只要每个线程私有,对sum数组还需初值
56 for(j=1;j<N;j++)
57 {
58 //本线程第一次运行时创建sum空间
59 if((j-1)%numThrdFor==0)
60 {
61 //创建存放自己线程临时结果的数组,初始化为0
62 //这里calloc还能像前面那样改进,但为了程序可读性暂时不做改进
63 sum=calloc(sizeof(double),N);
64 }
65
66 //对于这其中的每一短行
67 //(这个i被private保护,为每个线程单独创建了私有的i)
68 for(i=0;i<N;i++)
69 {
70 //如果用critical保护所有result[i]本质上这个计算完全是串行的
71 //想使用reduction子句,但是预编译只会做一次,而不会随for循环
72 //使用atomic使得这条语句是原子的,要执行就执行完而不拆分
73 //但用atomic这个计算也仍然是串行的,只是加解锁比critical快
74 //因此对每个线程拷贝分配一个求和数组,才能实现并行
75 sum[i]+=mat[i][j]*vec[j];//加到自己线程的sum数组上
76 }
77
78 //仅当本线程即将结束前,把算好的sum数组加上来
79 //判断条件:从1开始计数下能整除次数,或当前是最后一次循环
80 if(j%numThrdFor==0 || j==N-1)
81 {
82 for(i=0;i<N;i++)
83 # pragma omp atomic//在这里使用atomic保护result[i]
84 result[i]+=sum[i];
85 }
86 }
87
88 //采样输出结果看一下
89 for(i=0;i<N;i+=N/11)
90 printf("%f\n",result[i]);
91
92 return 0;
93 }
输出
1 [lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_col.o omp_col.c
2 [lzh@hostlzh OpenMP]$ ./omp_col.o 12
3 79800.000000
4 79800.000000
5 79800.000000
6 79800.000000
7 79800.000000
8 79800.000000
9 79800.000000
10 79800.000000
11 79800.000000
12 79800.000000
13 79800.000000
14 79800.000000
15 [lzh@hostlzh OpenMP]$
按块分配
每次循环都会重新开个线程，虽然操作了共享变量但是测试时没出现问题(为啥)。
1 #include<stdio.h>
2 #include<omp.h>
3 #include<stdlib.h>
4
5 #define N 400 //规模(方针的阶数)
6 int i,j;//通用游标
7 double **mat=NULL;//矩阵对象
8 int *vec=NULL;//列向量对象
9 double *result=NULL;//结果向量对象
10 int blkSqrt=-1;//块数的开算数平方
11
12 //自己的求平方根的函数,因为math里的sqrt报错
13 int mysqrt(int k)
14 {
15 int i;
16 for(i=0;i*i<k;i++)
17 ;
18 if(i*i==k)
19 return i;
20 return -1;
21 }
22
23 int main(int argc,char* argv[])
24 {
25 //读入10进制表示的进程数
26 int thrdCnt=strtol(argv[1],NULL,10);
27 //判断块数合法性
28 if(blkSqrt=mysqrt(thrdCnt)==-1)//块数不是完全平方数
29 {
30 printf("块数不是完全平方数!\n");
31 return 0;
32 }
33 //矩阵一级挂载空间
34 mat=(double**)malloc(N*sizeof(double *));
35 //存向量的空间
36 vec=(int *)malloc(N*sizeof(int));
37 //存结果的空间
38 result=(double *)malloc(N*sizeof(double));
39
40 //并行for:申请存矩阵的空间
41 //因为OpenMP不需要分发聚集等,所以不必做连续空间申请
42 # pragma omp parallel for num_threads(thrdCnt)
43 for(i=0;i<N;i++)
44 mat[i]=(double*)malloc(N*sizeof(double));
45
46 //并行for:为矩阵和列向量赋值(实际做时是从文件读入)
47 # pragma omp parallel for num_threads(thrdCnt) private(j)
48 for(i=0;i<N;i++)
49 {
50 vec[i]=1;
51 for(j=0;j<N;j++)
52 mat[i][j]=j;
53 }
54
55 //并行for:为结果向量赋初始值,这样能避免calloc时间或多余的if判断
56 # pragma omp parallel for num_threads(thrdCnt)
57 for(i=0;i<N;i++)
58 result[i]=mat[i][0]*vec[0];
59
60 //并行for:按块分配计算结果,即行列分别分开
61 # pragma omp parallel for num_threads(blkSqrt)
62 for(i=0;i<N;i++)
63 # pragma omp parallel for num_threads(blkSqrt)
64 for(j=1;j<N;j++)
65 result[i]+=mat[i][j]*vec[j];
66
67 //采样输出结果看一下
68 for(i=0;i<N;i+=N/11)
69 printf("%f\n",result[i]);
70
71 return 0;
72 }
输出
1 [lzh@hostlzh OpenMP]$ gcc -fopenmp -o omp_blk.o omp_blk.c
2 [lzh@hostlzh OpenMP]$ ./omp_blk.o 9
3 79800.000000
4 79800.000000
5 79800.000000
6 79800.000000
7 79800.000000
8 79800.000000
9 79800.000000
10 79800.000000
11 79800.000000
12 79800.000000
13 79800.000000
14 79800.000000
15 [lzh@hostlzh OpenMP]$
