首页 > 其他分享 >OpenMP与MPI混合做方阵向量乘法

OpenMP与MPI混合做方阵向量乘法

时间:2023-06-07 19:24:23浏览次数:38  
标签:NULL int MPI 进程 OpenMP 乘法 my 4950.000000

按行分配

  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]#

 

标签:NULL,int,MPI,进程,OpenMP,乘法,my,4950.000000
From: https://www.cnblogs.com/ybqjymy/p/17464321.html

相关文章

  • OpenMP 归约和reduction子句
    简述归约归约操作在MPI里也学过,不过那时候还不太熟悉这种操作。当时只知道MPI_Reduce可以把全局求和和集合通信封装起来,非常方便。实际上将相同的二元归约操作符重复地应用到一个序列上得到结果的计算过程都可以称为归约。python里那个难理解的reduce()函数也就是归约:1......
  • OpenMP 传统形式的方阵向量并行乘法
    按行分配思路和MPI基本类似,不过OpenMP是共享内存的,不必做分发和聚集,申请的矩阵空间就不必是完全连续的。1#include<stdio.h>2#include<omp.h>3#include<stdlib.h>45#defineN400//规模(方针的阶数)6inti,j;//通用游标7double**mat=NULL;//矩阵对象......
  • compiler expression pattern match
     编译器中经常需要用到patternmatch。那么如何实现呢?比较直观的方法是使用递归。以patternmatch:y=a*(b+c)为例。首先,将其解析成一个抽象语法树:*a+bc其次,递归match:match(y,pattern)=>match(y,'*a+bc')左边是待检测的string,右边的是pattern。只要......
  • Python正则表达式学习(3)——re.compile()
    re.compile(pattern,flags=0)将正则表达式pattern编译为正则表达式对象,可用于使用其match()和search()方法进行匹配。顺序:prog=re.compile(pattern)result=prog.match(string)等价于:result=re.match(pattern,string)但是当单个程序中的表达式被多次使用时,使用re.comp......
  • Codeforces 1588F - Jumping Through the Array
    显然无法用polylog的数据结构维护,序列分块也不行,考虑询问分块。每\(B\)个询问处理一次。将这个询问中\(2,3\)操作涉及到的点设为“关键点”,那么容易发现,环上每一段以关键点结尾的链在这块操作的过程中始终保持不变,也就是说我们可以把它们缩在一起。先预处理出每个块的增量......
  • java反编译工具jd-gui和插件jd-eclipse,还有插件Enhanced Class Decompiler 3.3.0
    JD-GUI和JD-ECLIPSE可以直接在下面的网址进行下载http://java-decompiler.github.io/ (1)注意:JD-GUI.exe单机版有很多版本,有些旧版本反编译出来的源码和高版本反编译出来的源码是区别的1.低版本的反编译可能和实际源码有出入2.1.6.6版本反编译的源码中有中文无法正常复制? ......
  • MPI-IS/mesh库的安装与使用
    MPI-IS/mesh库MeshProcessing Library是由德国马克斯·普朗克计算机科学研究所(MPI-IS)开发的一个开源网格处理库,用于处理三维网格数据。MPI-ISMeshProcessingLibrary提供了一系列的网格处理算法,包括网格滤波、网格重建、网格配准、网格切割、网格拓扑结构处理等。它支......
  • 【Exception】maven-compiler-plugin 编译失败集锦
    1JDK明明是1.8为什么说编译环境和运行环境不一致?Whatfuck?JDK明明1.8为什么编译环境变成1.5了?Whatfuck?原因分析:奇怪的是我的机器上只安装了JDK8,为什么还会说不支持diamond和lambda呢?在Google大神的指引下,在MavenCompiler插件介绍里面找到了答案:Alsonotethat......
  • java基础知识之 算法 九九乘法表
    /**auther:kevindate:20110710function:amultiplicationtableEditon:1rt*/importjava.util.Scanner;//progarmusesclassScannerpublicclassJiujiu{//mainmethodbeginsexecutionofJavaapplicationpublicstaticvoidmain(Stringargs[]){J......
  • 最小二乘法的矩阵正则化改进——“岭回归”和“LASSO回归”算法
    看代码过程中发现了一个很奇怪的概念,叫做“最小二乘法的矩阵正则化”,这个词汇十分的陌生,虽然最小二乘法是知道的,但是用了矩阵正则化的最小二乘法是个什么东西呢?  相关代码见:强化学习:连续控制问题中Actor-Critic算法的linearbaseline  后来在网上一通查才知道,原来“最小二乘法......