按行分配
1 #include<stdio.h>
2 #include<mpi.h>
3 #include<stdlib.h>
4 #include<omp.h>
5
6 #define N 100
7
8 //time_t start,end;//开始和结束时间
9 double start,end;
10
11 int main(int argc,char* argv[])
12 {
13 //读入10进制表示的进程数,用于OpenMP
14 int thrdCnt=strtol(argv[1],NULL,10);
15
16 int *vec=NULL;//列向量
17 double *mat=NULL;//自己进程的那部分矩阵
18 int my_rank;//自己进程的进程号
19 int comm_sz;//总的进程数目
20 int my_row;//本进程处理的行数
21 int i,j;//通用游标
22 double *result=NULL;//用来存本进程计算出的结果向量
23 double *all_rst=NULL;//只给0号进程存总的结果
24
25 /**********************/
26 MPI_Init(NULL,NULL);
27 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
28 MPI_Comm_size(MPI_COMM_WORLD,&comm_sz);
29
30 //计时开始
31 if(my_rank==0)
32 {
33 //start=time(NULL);
34 start=MPI_Wtime();
35 }
36
37 //本进程处理的行数就是总阶数/进程数
38 my_row=N/comm_sz;
39
40 //为每个进程都申请空间
41 mat=malloc(N*my_row*sizeof(double)); //my_row行的小矩阵
42 vec=malloc(N*sizeof(int)); //每个进程各自读入列向量
43 result=malloc(my_row*sizeof(double));//每个进程各自的结果向量
44
45 //赋值(省去读入的步骤,实际做时是读入自己那部分)
46 //因为parallel for指令对外层j块选,故对vec[j]的操作无冲突
47 //此外,mat[]的每个元素只被一个线程操作一次,无冲突
48 # pragma omp parallel for num_threads(thrdCnt) private(i)
49 for(j=0;j<N;j++) {
50 for(i=0;i<my_row;i++)
51 mat[i*N+j]=j;
52 vec[j]=1;
53 }
54
55 //计算:j=0时做的是赋值
56 //因为parallel for指令对外层i块选,故对result[i]的操作无冲突
57 # pragma omp parallel for num_threads(thrdCnt) private(j)
58 for(i=0;i<my_row;i++) {
59 //j=0时做的是赋值
60 result[i]=mat[i*N]*vec[0];
61 //j!=0时做的是加赋值
62 for(j=1;j<N;j++)
63 result[i]+=mat[i*N+j]*vec[j];
64 }
65
66 //聚集给0号进程
67 if(my_rank==0)
68 {
69 //先开辟存储空间
70 all_rst=malloc(N*sizeof(double));
71
72 //聚集
73 MPI_Gather
74 (
75 result, /*发送内容的地址*/
76 my_row, /*发送的长度*/
77 MPI_DOUBLE, /*发送的数据类型*/
78 all_rst, /*接收内容的地址*/
79 my_row, /*接收的长度*/
80 MPI_DOUBLE, /*接收的数据类型*/
81 0, /*接收至哪个进程*/
82 MPI_COMM_WORLD /*通信域*/
83 );
84 }
85 else
86 {
87 //聚集
88 MPI_Gather
89 (
90 result, /*发送内容的地址*/
91 my_row, /*发送的长度*/
92 MPI_DOUBLE, /*发送的数据类型*/
93 all_rst, /*接收内容的地址*/
94 my_row, /*接收的长度*/
95 MPI_DOUBLE, /*接收的数据类型*/
96 0, /*接收至哪个进程*/
97 MPI_COMM_WORLD /*通信域*/
98 );
99 }
100
101 //0号进程负责输出
102 if(my_rank==0)
103 {
104 //计时结束
105 //end=time(NULL);
106 end=MPI_Wtime();
107 //计算时间
108 //printf("time=%f\n",difftime(end,start));
109 printf("time=%e\n",end-start);
110 printf("The Result is:\n");
111 //改变跨度可以采样获取结果,快速结束I/O
112 for(i=0;i<N;i+=N/11)
113 printf("%f\n",all_rst[i]);
114 }
115
116 MPI_Finalize();
117 /**********************/
118
119 //最终,free应无遗漏
120 free(all_rst);
121 free(mat);
122 free(vec);
123 free(result);
124
125 return 0;
126 }
输出
1 [root@hostlzh mpi_omp]# mpicc -fopenmp -o row.o row.c
2 [root@hostlzh mpi_omp]# mpiexec -n 4 ./row.o 7
3 time=7.662773e-03
4 The Result is:
5 4950.000000
6 4950.000000
7 4950.000000
8 4950.000000
9 4950.000000
10 4950.000000
11 4950.000000
12 4950.000000
13 4950.000000
14 4950.000000
15 4950.000000
16 4950.000000
17 [root@hostlzh mpi_omp]#
按列分配
1 #include<stdio.h>
2 #include<mpi.h>
3 #include<stdlib.h>
4 #include<omp.h>
5
6 #define N 100
7
8 //time_t start,end;//开始和结束时间
9 double start,end;
10
11 int main(int argc,char* argv[])
12 {
13 //读入10进制表示的进程数,用于OpenMP
14 int thrdCnt=strtol(argv[1],NULL,10);
15
16 int *vec=NULL;//列向量
17 double *mat=NULL;//自己进程的那部分矩阵
18 int my_rank;//自己进程的进程号
19 int comm_sz;//总的进程数目
20 int my_col;//本进程处理的列数
21 int i,j;//通用游标
22 double *result=NULL;//用来存本进程计算出的结果向量
23 double *all_rst=NULL;//只给0号进程存总的结果
24
25 /**********************/
26 MPI_Init(NULL,NULL);
27 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
28 MPI_Comm_size(MPI_COMM_WORLD,&comm_sz);
29
30 //计时开始
31 if(my_rank==0)
32 {
33 //start=time(NULL);
34 start=MPI_Wtime();
35 }
36
37 //本进程处理的列数就是总阶数/进程数
38 my_col=N/comm_sz;
39
40 //为每个进程都申请空间
41 mat=malloc(my_col*N*sizeof(double)); //my_col列的小矩阵
42 vec=malloc(my_col*sizeof(int)); //每个进程各自读入一段列向量
43 result=malloc(N*sizeof(double));//每个进程各自的结果向量
44
45 //赋值(省去读入的步骤,实际做时是读入自己那部分)
46 //因为parallel for指令对外层i块选,故对vec[i]的操作无冲突
47 //此外,mat[]的每个元素只被一个线程操作一次,无冲突
48 # pragma omp parallel for num_threads(thrdCnt) private(j)
49 for(i=0;i<my_col;i++) {
50 for(j=0;j<N;j++)
51 mat[j*my_col+i]=my_rank*my_col+i;//注意和my_rank有关
52 vec[i]=1;
53 }
54
55 //计算
56 //因为parallel for指令对外层i块选,故对result[i]的操作无冲突
57 # pragma omp parallel for num_threads(thrdCnt) private(j)
58 for(i=0;i<N;i++) {
59 //j=0时做的是赋值
60 result[i]=mat[i*my_col]*vec[0];
61 for(j=1;j<my_col;j++)
62 result[i]+=mat[i*my_col+j]*vec[j];
63 }
64
65 //归约给0号进程
66 if(my_rank==0)
67 {
68 //先开辟存储空间
69 all_rst=malloc(N*sizeof(double));
70 //将0号进程自己的result写入
71 for(i=0;i<N;i++)
72 all_rst[i]=result[i];
73
74 /*
75 一种改进的想法是,用0号进程自己的result数组,
76 作为MPI_Reduce归约的第2个参数.
77 书上指出这种方式可能会出现混乱,不建议我们这样做
78 */
79
80 //归约
81 MPI_Reduce
82 (
83 result, /*发送内容的地址*/
84 all_rst, /*接收内容的地址*/
85 N, /*归约操作的长度*/
86 MPI_DOUBLE, /*归约的数据类型*/
87 MPI_SUM, /*归约操作符*/
88 0, /*接收至哪个进程*/
89 MPI_COMM_WORLD /*通信域*/
90 );
91 }
92 else
93 {
94 //归约
95 MPI_Reduce
96 (
97 result, /*发送内容的地址*/
98 all_rst, /*接收内容的地址*/
99 N, /*归约操作的长度*/
100 MPI_DOUBLE, /*归约的数据类型*/
101 MPI_SUM, /*归约操作符*/
102 0, /*接收至哪个进程*/
103 MPI_COMM_WORLD /*通信域*/
104 );
105 }
106
107 //0号进程负责输出
108 if(my_rank==0)
109 {
110 //计时结束
111 //end=time(NULL);
112 end=MPI_Wtime();
113 //计算时间
114 //printf("time=%f\n",difftime(end,start));
115 printf("time=%e\n",end-start);
116 printf("The Result is:\n");
117 //改变跨度可以采样获取结果,快速结束I/O
118 for(i=0;i<N;i+=N/11)
119 printf("%f\n",all_rst[i]);
120 }
121
122 MPI_Finalize();
123 /**********************/
124
125 //最终,free应无遗漏
126 free(all_rst);
127 free(mat);
128 free(vec);
129 free(result);
130
131 return 0;
132 }
输出
1 [root@hostlzh mpi_omp]# mpicc -fopenmp -o col.o col.c
2 [root@hostlzh mpi_omp]# mpiexec -n 10 ./col.o 3
3 time=4.535699e-02
4 The Result is:
5 4950.000000
6 4950.000000
7 4950.000000
8 4950.000000
9 4950.000000
10 4950.000000
11 4950.000000
12 4950.000000
13 4950.000000
14 4950.000000
15 4950.000000
16 4950.000000
17 [root@hostlzh mpi_omp]#
按块分配
1 #include<stdio.h>
2 #include<mpi.h>
3 #include<stdlib.h>
4 #include<omp.h>
5
6 #define N 100
7
8 //time_t start,end;//开始和结束时间
9 double start,end;
10
11 //自己的求平方根的函数,因为math里的sqrt报错
12 int mysqrt(int k)
13 {
14 int i;
15 for(i=0;i*i<k;i++)
16 ;
17 if(i*i==k)
18 return i;
19 return -1;
20 }
21
22 int main(int argc,char* argv[])
23 {
24 //读入10进制表示的进程数,用于OpenMP
25 int thrdCnt=strtol(argv[1],NULL,10);
26
27 int *vec=NULL;//列向量
28 double *mat=NULL;//自己进程的那部分矩阵
29 int my_rank;//自己进程的进程号
30 int comm_sz;//总的进程数目
31 int i,j,k;//通用游标
32 double *result=NULL;//用来存本进程计算出的结果向量
33 double *all_rst=NULL;//只给0号进程存总的结果
34 int t;//每行(列)的子阵块数
35 int my_N;//每行(列)的子阵块维度
36
37 /**********************/
38 MPI_Init(NULL,NULL);
39 MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
40 MPI_Comm_size(MPI_COMM_WORLD,&comm_sz);
41
42 //计时开始
43 if(my_rank==0)
44 {
45 //start=time(NULL);
46 start=MPI_Wtime();
47 }
48
49 t=mysqrt(comm_sz);//一行(列)有多少个子阵块
50 if(t==-1)
51 {
52 if(my_rank==0)//只给0号进程输出的权力
53 printf("t is -1\n");
54 return 0;//但所有进程遇此情况都应结束
55 }
56 //本进程处理的块阶数=总阶数/sqrt(进程数)
57 my_N=N/t;
58
59 //为每个进程都申请空间
60 mat=malloc(my_N*my_N*sizeof(double)); //my_N阶的小方阵
61 vec=malloc(my_N*sizeof(int)); //每个进程各自读入自己需要的那段列向量
62 result=malloc(my_N*sizeof(double));//每个进程各自的结果向量
63
64 //赋值(省去读入的步骤,实际做时是读入自己那部分)
65 //因为parallel for指令对外层i块选,故对vec[i]的操作无冲突
66 //此外,mat[]的每个元素只被一个线程操作一次,无冲突
67 # pragma omp parallel for num_threads(thrdCnt) private(j)
68 for(i=0;i<my_N;i++) {
69 for(j=0;j<my_N;j++)
70 mat[i*my_N+j]=my_rank%t*my_N+j;
71 vec[i]=1;
72 }
73
74 //计算
75 //因为parallel for指令对外层i块选,故对result[i]的操作无冲突
76 # pragma omp parallel for num_threads(thrdCnt) private(j)
77 for(i=0;i<my_N;i++) {
78 //j=0时做的是赋值
79 result[i]=mat[i*my_N]*vec[0];
80 for(j=1;j<my_N;j++)
81 result[i]+=mat[i*my_N+j]*vec[j];
82 }
83
84 /*聚集到0号进程上
85 实际上是开了comm_sz个my_N长度组成的总向量作聚集*/
86 if(my_rank==0)
87 {
88 //结果向量一共是comm_sz个长度为my_N的子向量
89 all_rst=malloc(my_N*comm_sz*sizeof(double));
90 //聚集
91 MPI_Gather(result,
92 my_N,
93 MPI_DOUBLE,
94 all_rst,
95 my_N,
96 MPI_DOUBLE,
97 0,
98 MPI_COMM_WORLD
99 );
100 }
101 else
102 //聚集
103 MPI_Gather(result,
104 my_N,
105 MPI_DOUBLE,
106 all_rst,
107 my_N,
108 MPI_DOUBLE,
109 0,
110 MPI_COMM_WORLD
111 );
112
113 /*从结果长向量中计算结果
114 把后面的分量加到0~my_N-1分量上去*/
115 if(my_rank==0)
116 {
117 //后面的分量加上来
118 //因为每个线程的i范围无重叠,且N<my_N
119 //因此对all_rst[]数组中的i*N+k下标所在内存空间不需各自保护
120 # pragma omp parallel for num_threads(thrdCnt) private(j,k)
121 for(i=0;i<t;i++) //一共有t行
122 {
123 for(j=1;j<t;j++) //每行有t个
124 {
125 for(k=0;k<my_N;k++)
126 {
127 all_rst[i*N+k]+=all_rst[i*N+j*my_N+k];
128 }
129 }
130 }
131
132 //补到前面
133 //和前面同样的原因,对all_rst数组不需保护
134 # pragma omp parallel for num_threads(thrdCnt) private(k)
135 for(i=1;i<t;i++)
136 {
137 for(k=0;k<my_N;k++)
138 {
139 all_rst[i*my_N+k]=all_rst[i*N+k];
140 }
141 }
142
143 //计时结束
144 //end=time(NULL);
145 end=MPI_Wtime();
146 //计算时间
147 //printf("time=%f\n",difftime(end,start));
148 printf("time=%e\n",end-start);
149 printf("The Result is:\n");
150 //改变跨度可以采样获取结果,快速结束I/O
151 for(i=0;i<N;i+=N/11)
152 printf("%f\n",all_rst[i]);
153
154 }
155
156 MPI_Finalize();
157 /**********************/
158
159 //最终,free应无遗漏
160 free(all_rst);
161 free(mat);
162 free(vec);
163 free(result);
164
165 return 0;
166 }
输出
1 [root@hostlzh mpi_omp]# mpicc -fopenmp -o block.o block.c
2 [root@hostlzh mpi_omp]# mpiexec -n 4 ./block.o 5
3 time=2.295089e-02
4 The Result is:
5 4950.000000
6 4950.000000
7 4950.000000
8 4950.000000
9 4950.000000
10 4950.000000
11 4950.000000
12 4950.000000
13 4950.000000
14 4950.000000
15 4950.000000
16 4950.000000
17 [root@hostlzh mpi_omp]#