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

并行计算

开发平台:

Unix_Linux

  1. #include "stdio.h"
  2. #include "stdlib.h"
  3. #include "mpi.h"
  4. #include "math.h"
  5. #define E 0.0001
  6. #define a(x,y) a[x*size+y]
  7. #define b(x) b[x]
  8. #define v(x) v[x]
  9. /*A为size*size的系数矩阵*/
  10. #define A(x,y) A[x*size+y]
  11. #define B(x) B[x]
  12. #define V(x) V[x]
  13. #define intsize sizeof(int)
  14. #define floatsize sizeof(float)
  15. #define charsize sizeof(char)
  16. int size,N;
  17. int m;
  18. float *B;
  19. float *A;
  20. float *V;
  21. int my_rank;
  22. int p;
  23. MPI_Status status;
  24. FILE *fdA,*fdB,*fdB1;
  25. int i,j,my_rank,group_size;
  26. int k;
  27. float *sum;
  28. float *b;
  29. float *v;
  30. float *a;
  31. float *differ;
  32. float temp;
  33. int iteration,total,loop;
  34. int r,q;
  35. /*
  36.  * 函数名:Environment_Finalize
  37.  * 功能:释放各个进程的a,b,v
  38.  */
  39. void Environment_Finalize(float *a,float *b,float *v)
  40. {
  41.     free(a);
  42.     free(b);
  43.     free(v);
  44. }
  45. /*
  46.  * 函数名: read_data()
  47.  * 功能:读取数据
  48.  */
  49. void read_data()
  50. {
  51.     if(my_rank==0)
  52.     {
  53.         /*打开数据文件*/
  54.         fdA=fopen("dataIn.txt","r");
  55.         fscanf(fdA,"%d %d", &size, &N);
  56.         if (size != N-1)
  57.         {
  58.             printf("the input is wrongn");
  59.             exit(1);
  60.         }
  61.         /*分配空间*/
  62.         A=(float *)malloc(floatsize*size*size);
  63.         B=(float *)malloc(floatsize*size);
  64.         V=(float *)malloc(floatsize*size);
  65.         /*初始化矩阵A和向量B*/
  66.         for(i = 0; i < size; i++)
  67.         {
  68.             for(j = 0; j < size; j++)
  69.                 fscanf(fdA,"%f", A+i*size+j);
  70.             fscanf(fdA,"%f", B+i);
  71.         }
  72.         /*初始化向量V*/
  73.         for(i = 0; i < size; i ++)
  74.             fscanf(fdA,"%f", V+i);
  75.         /*关闭数据文件*/
  76.         fclose(fdA);
  77.         /*输出数据文件信息*/
  78.         printf("Input of file "dataIn.txt"n");
  79.         printf("%dt%dn", size, N);
  80.         /*输出矩阵A和向量B的值*/
  81.         for(i = 0; i < size; i ++)
  82.         {
  83.             for(j = 0; j < size; j ++)
  84.                 printf("%ft",A(i,j));
  85.             printf("%fn",B(i));
  86.         }
  87.         printf("n");
  88.         /*输出解向量V的值*/
  89.         for(i = 0; i < size; i ++)
  90.             printf("%ft", V(i));
  91.         printf("nn");
  92.         printf("nOutput of resultn");
  93.     }
  94. }
  95. /*
  96.  * 函数名: broadcast_data()
  97.  * 功能:广播数据
  98.  */
  99. void broadcast_data()
  100. {
  101.     /*0号进程将size广播给所有进程*/
  102.     MPI_Bcast(&size,1,MPI_INT,0,MPI_COMM_WORLD);
  103.     m=size/p;
  104.     if (size%p!=0) m++;
  105.     /*为每个进程分配需要的空间*/
  106.     v=(float *)malloc(floatsize*size);
  107.     a=(float *)malloc(floatsize*m*size);
  108.     b=(float *)malloc(floatsize*m);
  109.     sum=(float *)malloc(floatsize*m);
  110.     /*判断是否分配成功*/
  111.     if (a==NULL||b==NULL||v==NULL)
  112.     {
  113.         printf("allocate space  fail!");
  114.         exit(1);
  115.     }
  116.     /*如果是0号进程,就将本进程的解向量初始化为问题的解向量*/
  117.     if (my_rank==0)
  118.     {
  119.         for(i=0;i<size;i++)
  120.             v(i)=V(i);
  121.     }
  122.     /*初始解向量v被广播给所有进程*/
  123.     MPI_Bcast(v,size,MPI_FLOAT,0,MPI_COMM_WORLD);
  124.     /*0号进程采用行块划分将矩阵A划分为大小为m*size的p块子矩阵,
  125.       将B划分为大小为m的p块子向量,依次发送给1至p-1号进程*/
  126.     if (my_rank==0)
  127.     {
  128.         /*初始化0号进程的a和b*/
  129.         for(i=0;i<m;i++)
  130.             for(j=0;j<size;j++)
  131.                 a(i,j)=A(i,j);
  132.         for(i=0;i<m;i++)
  133.             b(i)=B(i);
  134.         /*将各自的a和b发送到1至p-1号进程*/
  135.         for(i=1;i<p;i++)
  136.         {
  137.             MPI_Send(&(A(m*i,0)), m*size,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  138.             MPI_Send(&(B(m*i)),m, MPI_FLOAT,i,i,MPI_COMM_WORLD);
  139.         }
  140.         /*释放空间*/
  141.         free(A); free(B); free(V);
  142.     }
  143.     else
  144.     {
  145.         /*如果不是0号进程,就接收从0号进程传来的子矩阵a和子向量b*/
  146.         MPI_Recv(a,m*size,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD, &status);
  147.         MPI_Recv(b,m,MPI_FLOAT,0,my_rank, MPI_COMM_WORLD,&status);
  148.     }
  149. }
  150. /*
  151.  * 函数名: cal_cross_data()
  152.  * 功能:主对角元素右边的数据求和
  153.  */
  154. void cal_cross_data()
  155. {
  156.     /*所有进程并行地对主对角元素右边的数据求和*/
  157.     for(i=0;i<m;i++)
  158.     {
  159.         sum[i]=0.0;
  160.         for(j=0;j<size;j++)
  161.             if (j>(my_rank*m+i))
  162.                 sum[i]=sum[i]+a(i,j)*v(j);
  163.     }
  164. }
  165. /*
  166.  * 函数名: die_dai_cal()
  167.  * 功能:各个进程并行迭代求解
  168.  */
  169. void die_dai_cal()
  170. {
  171.     /*total为新旧值之差小于E的x的分量的个数*/
  172.     while (total<size)
  173.     {
  174.         iteration=0;
  175.         /*iteration为本进程中新旧值之差小于E的x的分量个数*/
  176.         total=0;
  177.         for(j=0;j<size;j++)
  178.         {
  179.             r=j%m;
  180.             q=j/m;
  181.             /*编号为q的进程负责计算解向量并广播给所有进程*/
  182.             if (my_rank==q)                       /*主行所在的进程*/
  183.             {
  184.                 temp=v(my_rank*m+r);
  185.                 for(i=0;i<r;i++)
  186.                     sum[r]=sum[r]+a(r,my_rank*m+i)*v(my_rank*m+i);
  187.                 /*计算出的解向量*/
  188.                 v(my_rank*m+r)=(b(r)-sum[r]) / a(r,my_rank*m +r);
  189.                 if (fabs(v(my_rank*m+r)-temp)<E)
  190.                     iteration++;
  191.                 /*广播解向量*/
  192.                 MPI_Bcast(&v(my_rank*m+r), 1,MPI_FLOAT,my_rank,MPI_COMM_WORLD);
  193.                 sum[r]=0.0;
  194.                 for(i=0;i<r;i++)
  195.                     sum[i]=sum[i]+a(i,my_rank*m+r)*v(my_rank*m+r);
  196.             }
  197.             else                                  /*其他进程接受广播来的x[j]的新值*/
  198.             /*各进程对解向量的其它分量计算有关乘积项并求和*/
  199.             {
  200.                 MPI_Bcast(&v(q*m+r),1, MPI_FLOAT,q,MPI_COMM_WORLD);
  201.                 for(i=0;i<m;i++)
  202.                     sum[i]=sum[i]+a(i,q*m +r)*v(q*m+r);
  203.             }
  204.         }
  205.         /*通过归约操作的求和运算以决定是否进行下一次迭代*/
  206.         /*allreduce过程求出所有进程iteration和total,并广播到所有进程*/
  207.         MPI_Allreduce(&iteration,&total,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
  208.         loop++;
  209.         if (my_rank==0)
  210.             printf("in the %d times total vaule = %dn",loop,total);
  211.     }
  212. }
  213. /*
  214.  * 函数名: put_data
  215.  * 功能:输出数据
  216.  */
  217. void put_data()
  218. {
  219.     if (my_rank==0)
  220.     {
  221.         for(i = 0; i < size; i ++)
  222.             printf("x[%d] = %fn",i,v(i));
  223.         printf("n");
  224.         /*输出迭代次数*/
  225.         printf("Iteration num = %dn",loop);
  226.     }
  227. }
  228. /*
  229.  * 函数名: main
  230.  * 功能:主程序
  231.  * 输入:argc为命令行参数个数;
  232.  *       argv为每个命令行参数组成的字符串数组。
  233.  * 输出:返回0代表程序正常结束;
  234.  */
  235. int main(int argc, char **argv)
  236. {
  237.     loop = total = 0;
  238.     /*MPI程序初始化*/
  239.     MPI_Init(&argc,&argv);
  240.     MPI_Comm_size(MPI_COMM_WORLD, &group_size);
  241.     MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
  242.     p=group_size;
  243.     /*0号进程负责读取数据*/
  244.     read_data();
  245.     /*广播数据*/
  246.     broadcast_data();
  247.     /*所有进程并行地对主对角元素右边的数据求和*/
  248.     cal_cross_data();
  249.     /*迭代求解*/
  250.     die_dai_cal();
  251.     /*0号进程将结果输出*/
  252.     put_data();
  253.     MPI_Barrier(MPI_COMM_WORLD);
  254.     MPI_Finalize();
  255.     Environment_Finalize(a,b,v);
  256.     return (0);
  257. }