closure.c
上传用户:hhyinxing
上传日期:2013-09-10
资源大小:833k
文件大小:5k
源码类别:

并行计算

开发平台:

Unix_Linux

  1. #include "math.h"
  2. #include "stdio.h"
  3. #include "stdlib.h"
  4. #include "mpi.h"
  5. #define INTSIZE sizeof(int)
  6. #define CHARSIZE sizeof(char)
  7. #define M(i,j) M[i*n*p+j]
  8. #define a(i,j) a[i*n*p+j]
  9. #define b(i,j) b[i*n+j]
  10. #define c(l,i,j) c[i*n*p+l*n+j]
  11. int *a,*b,*c,*d;
  12. /**输入的邻接矩阵**/
  13. int *M;
  14. /**处理器数目**/
  15. int p;
  16. /**顶点数目**/
  17. int N;
  18. /**分配该单个处理器的顶点数目**/
  19. int n;
  20. /**处理器编号**/
  21. int myid;
  22. MPI_Status status;
  23. /* 
  24.  * 函数名: readmatrix 
  25.  * 功能: 读入邻接矩阵
  26.  * 输入: 
  27.  * 输出:返回0代表程序正常结束;否则程序出错。
  28.  */ 
  29. int readmatrix()
  30. {
  31. int i,j;
  32. printf("Input the size of matrix:");
  33. scanf("%d",&N);
  34. n=N/p;
  35. if(N%p!=0)n++;
  36. /*给M矩阵分配空间*/
  37. M=(int*)malloc(INTSIZE*n*p*n*p);
  38. if(!M)
  39. {
  40. error("failed to allocate space for M");
  41. }/*if*/
  42. printf("Input the matrix:n");
  43. /*输入矩阵*/
  44. for(i=0;i<N;i++)
  45. {
  46. for(j=0;j<N;j++)
  47. {
  48. scanf("%d",&(M(i,j)));
  49. if (i==j) M(i,j)=1;
  50. }/*for*/
  51. for(j=N;j<n*p;j++) M(i,j)=0;
  52. }/*for*/
  53. for(i=N;i<n*p;i++)
  54. for(j=0;j<n*p;j++)
  55. if(i==j) M(i,j)=1;
  56. else M(i,j)=0;
  57. return(0);
  58. } /*readmatrix*/
  59. /* 
  60.  * 函数名: error 
  61.  * 功能: 提示错误并退出
  62.  * 输入:message为要提示的消息
  63.  * 输出:输出出错的处理器号和消息
  64.  *       返回0代表程序正常结束;否则程序出错。
  65.  */ 
  66. int error(message)
  67. char *message;
  68. {
  69. printf("processor %d:%sn",myid,message);
  70. /*退出*/
  71. exit(0);
  72. return(0);
  73. } /*error*/
  74. /* 
  75.  * 函数名: sendmatrix 
  76.  * 功能: 向各处理器发送邻接矩阵数据
  77.  * 输入: 
  78.  * 输出:返回0代表程序正常结束;否则程序出错。
  79.  */ 
  80. int sendmatrix()
  81. {
  82. int i,j;
  83. for(i=1;i<p;i++)
  84. {
  85. MPI_Send(&(M(n*i,0)),n*n*p,MPI_INT,i,i,MPI_COMM_WORLD);
  86. for(j=0;j<n*p;j++)
  87. MPI_Send(&(M(j,n*i)),n,MPI_INT,i,i,MPI_COMM_WORLD);
  88. } /*for*/
  89. return(0);
  90. } /*sendmatrix*/
  91. /* 
  92.  * 函数名: getmatrix 
  93.  * 功能: 各处理器接收邻接矩阵数据
  94.  * 输入: 
  95.  * 输出:返回0代表程序正常结束;否则程序出错。
  96.  */ 
  97. int getmatrix()
  98. {
  99. int i,j;
  100. if(myid==0)
  101. {
  102. for(i=0;i<n;i++)
  103. for(j=0;j<n*p;j++)
  104. a(i,j)=M(i,j);
  105. for(i=0;i<n*p;i++)
  106. for(j=0;j<n;j++)
  107. b(i,j)=M(i,j);
  108. }
  109. else
  110. {
  111. MPI_Recv(a,n*n*p,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
  112. for(j=0;j<n*p;j++)
  113. MPI_Recv(&(b(j,0)),n,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
  114. } /*if*/
  115. return(0);
  116. }
  117. /* 
  118.  * 函数名: writeback 
  119.  * 功能: 并行矩阵运算过程,将结果保存在c矩阵中
  120.  * 输入: 
  121.  * 输出:返回0代表程序正常结束;否则程序出错。
  122.  */ 
  123. int paramul()
  124. {
  125. int l,i,j,s;
  126. /**以下两变量表示前一个和后一个处理器标识**/
  127. int last,next;
  128. for(l=0;l<p;l++)
  129. {
  130. for(i=0;i<n;i++)
  131. for(j=0;j<n;j++)
  132. for(c(((l+myid)%p),i,j)=0,s=0;s<n*p;s++)
  133. if(a(i,s)&&b(s,j))
  134. {c(((l+myid)%p),i,j)=1;break;}
  135. last=(p+myid-1)%p;
  136. next=(myid+1)%p;
  137. if(l!=p-1)
  138. {
  139. if(myid%2==0)
  140. /**偶数号处理器直接接收数据b**/
  141. {
  142. MPI_Send(b,n*n*p,MPI_INT,last,last,MPI_COMM_WORLD);
  143. MPI_Recv(b,n*n*p,MPI_INT,next,myid,MPI_COMM_WORLD,&status);
  144. }
  145. else
  146. {
  147. /**奇数号处理器需要先将原数据保存,防止丢失**/
  148. for(i=0;i<n*n*p;i++)
  149. d[i]=b[i];
  150. MPI_Recv(b,n*n*p,MPI_INT,next,myid,MPI_COMM_WORLD,&status);
  151. MPI_Send(d,n*n*p,MPI_INT,last,last,MPI_COMM_WORLD);
  152. }
  153. }
  154. }/*for*/
  155. return(0);
  156. }
  157. /* 
  158.  * 函数名: writeback 
  159.  * 功能: 将c矩阵保存的矩阵相乘结果写回处理器0的数组M
  160.  * 输入: 
  161.  * 输出:返回0代表程序正常结束;否则程序出错。
  162.  */ 
  163. int writeback()
  164. int i;
  165. if(myid==0)
  166. {
  167. for(i=0;i<n*n*p;i++)
  168. M(0,i)=c(0,0,i);
  169. for(i=1;i<p;i++)
  170. {
  171. MPI_Recv(&(M(i*n,0)),n*n*p,MPI_INT,i,i,MPI_COMM_WORLD,&status);
  172. }/*for*/
  173. }/*for*/
  174. else
  175. MPI_Send(c,n*n*p,MPI_INT,0,myid,MPI_COMM_WORLD);
  176. return(0);
  177. }
  178. /* 
  179.  * 函数名: output 
  180.  * 功能: 输出结果
  181.  * 输入: 
  182.  * 输出:输出结果矩阵M
  183.  *       返回0代表程序正常结束;否则程序出错。
  184.  */ 
  185. int output()
  186. {
  187. int i,j;
  188. for(i=0;i<N;i++)
  189. {
  190. for(j=0;j<N;j++)
  191. printf("%d ",M(i,j));
  192. printf("n");
  193. }
  194. return(0);
  195. }
  196. /******************** main ********************/
  197. /* 
  198.  * 函数名: main 
  199.  * 功能: 启动MPI计算;
  200.  *       确定进程数和进程标识符;
  201.  *       调用主进程和从进程程序并行求解传递闭包问题。
  202.  * 输入:argc为命令行参数个数; 
  203.  *       argv为每个命令行参数组成的字符串数组。 
  204.  * 输出:返回0代表程序正常结束;否则程序出错。
  205.  */ 
  206. int main(argc,argv)
  207. int argc;
  208. char **argv;
  209. {
  210. int i,j;
  211. int group_size;
  212. int *temp;
  213.     /*初始化*/
  214. MPI_Init(&argc,&argv);
  215. /*确定工作组中的进程个数*/
  216. MPI_Comm_size(MPI_COMM_WORLD,&group_size);
  217. /*确定自己在工作组中的进程标识符*/
  218. MPI_Comm_rank(MPI_COMM_WORLD,&myid);
  219. p=group_size;
  220. /**处理器0读入矩阵并发送,步骤(1)**/
  221. if(myid==0)
  222. {
  223. /*输入邻接矩阵*/
  224. readmatrix();
  225. for(i=1;i<p;i++)
  226. MPI_Send(&N,1,MPI_INT,i,i,MPI_COMM_WORLD);
  227. }
  228. else
  229. MPI_Recv(&N,1,MPI_INT,0,myid,MPI_COMM_WORLD,&status);
  230. n=N/p;
  231. if(N%p!=0)n++;
  232. /*分配存储空间*/
  233. a=(int*)malloc(INTSIZE*n*n*p);
  234. b=(int*)malloc(INTSIZE*n*n*p);
  235. c=(int*)malloc(INTSIZE*n*n*p);
  236. d=(int*)malloc(INTSIZE*n*n*p);
  237. /*分配失败*/
  238. if(a==NULL||b==NULL||c==NULL)
  239. error("failed to allocate space for a,b and c");
  240. /**logN次矩阵自乘,步骤(3)**/
  241. for(i=0;i<=log(N)/log(2);i++)
  242. {
  243. if(myid==0)printf("loop %d:n",i);
  244. if(myid==0) output();
  245. /**处理器0向各处理器发送必要数据,步骤(3.1)(3.2)**/
  246. if(myid==0) sendmatrix();
  247. /**步骤(3.3)**/
  248. getmatrix();
  249. /**矩阵相乘,步骤(3.4)**/
  250. paramul();
  251. /**结果写回,步骤(3.5)**/
  252. writeback();
  253. };
  254. if(myid==0)printf("loop %d:n",i);
  255. if(myid==0)output();
  256. /*结束MPI计算*/
  257. MPI_Finalize();
  258. /*释放存储空间*/
  259. free(a);
  260. free(b);
  261. free(c);
  262. free(M);
  263. return(0);
  264. } /* main */