relaxation.c
上传用户:jsljixie
上传日期:2013-08-15
资源大小:827k
文件大小:7k
源码类别:

并行计算

开发平台:

Visual C++

  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. #define A(x,y) A[x*size+y]
  10. #define B(x) B[x]
  11. #define V(x) V[x]
  12. #define intsize sizeof(int)
  13. #define floatsize sizeof(float)
  14. #define charsize sizeof(char)
  15. int size,N;                                       /*size 为系数矩阵的大小,N为size+1*/
  16. int m;
  17. float *B;
  18. float *A;
  19. float *V;                                         /* store the value x */
  20. double starttime;
  21. double time1;
  22. double time2;
  23. int my_rank;
  24. int p;
  25. MPI_Status status;
  26. FILE *fdA,*fdB,*fdB1;
  27. void Environment_Finalize(float *a,float *b,float *v)
  28. {
  29.     free(a);
  30.     free(b);
  31.     free(v);
  32. }
  33. /*
  34.  * 函数名: main
  35.  * 输入:argc为命令行参数个数;
  36.  *       argv为每个命令行参数组成的字符串数组。
  37.  * 输出:返回0代表程序正常结束
  38.  */
  39. int main(int argc, char **argv)
  40. {
  41.     int i,j,my_rank,group_size;
  42.     int k;
  43.     float *sum;
  44.     float *b;
  45.     float *v;
  46.     float *a;
  47.     float *differ;
  48.     float temp,w;
  49.     int iteration,total,loop;
  50.     int r,q;
  51.     loop=0;
  52.     total=0;
  53.     w=0.9;
  54.     MPI_Init(&argc,&argv);                        /* 启动计算 */
  55.     MPI_Comm_size(MPI_COMM_WORLD,&group_size);    /* 找进程数 */
  56.     MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);       /* 找自己的id */
  57.     p=group_size;
  58.     if(my_rank==0)
  59.     {
  60.         starttime=MPI_Wtime();                    /* 开始计时 */
  61.         fdA=fopen("dataIn.txt","r");              /* 读入输入的系数矩阵A,b */
  62.         fscanf(fdA,"%d %d", &size, &N);
  63.         /*size 要求 size=N+1*/
  64.         if (size != N-1)
  65.         {
  66.             printf("the input is wrongn");
  67.             exit(1);
  68.         }
  69.         /*分配系数矩阵A,与b,解x的空间*/
  70.         A=(float *)malloc(floatsize*size*size);
  71.         B=(float *)malloc(floatsize*size);
  72.         V=(float *)malloc(floatsize*size);
  73.         /*存放A,b,x*/
  74.         for(i = 0; i < size; i++)
  75.         {
  76.             for(j = 0; j < size; j++)
  77.             {
  78.                 fscanf(fdA,"%f", A+i*size+j);
  79.             }
  80.             fscanf(fdA,"%f", B+i);
  81.         }
  82.         for(i = 0; i < size; i ++)
  83.         {
  84.             fscanf(fdA,"%f", V+i);
  85.         }
  86.         fclose(fdA);
  87.         /*打印A,b,x*/
  88.         printf("Input of file "dataIn.txt"n");
  89.         printf("%dt%dn", size, N);
  90.         for(i = 0; i < size; i ++)
  91.         {
  92.             for(j = 0; j < size; j ++) printf("%ft",A(i,j));
  93.             printf("%fn",B(i));
  94.         }
  95.         printf("n");
  96.         for(i = 0; i < size; i ++)
  97.         {
  98.             printf("%ft", V(i));
  99.         }
  100.         printf("nn");
  101.         printf("nOutput of resultn");
  102.     }                                             /* end of if (my_rank==0) */
  103.     MPI_Bcast(&size,1,MPI_INT,0,MPI_COMM_WORLD);  /* 广播size到所有进程*/
  104.     m=size/p;if (size%p!=0) m++;                  /* 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.     if (a==NULL||b==NULL||v==NULL)
  111.         printf("allocate space  fail!");
  112.     if (my_rank==0)
  113.     {
  114.         for(i=0;i<size;i++)
  115.             v(i)=V(i);
  116.     }
  117.     MPI_Bcast(v,size,MPI_FLOAT,0,MPI_COMM_WORLD); /* 将初始x广播到所有进程 */
  118.     /* 给进程分配A,b,和x初值*/
  119.     if (my_rank==0)
  120.     {
  121.         for(i=0;i<m;i++)
  122.             for(j=0;j<size;j++)
  123.                 a(i,j)=A(i,j);
  124.         for(i=0;i<m;i++)
  125.             b(i)=B(i);
  126.     }
  127.     /* 将A,b和x初值分配到所有进程,每个进程分配m行 */
  128.     if (my_rank==0)                               /* 进程0将数据广播出去 */
  129.     {
  130.         for(i=1;i<p;i++)
  131.         {
  132.             MPI_Send(&(A(m*i,0)),m*size,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  133.             MPI_Send(&(B(m*i)),m,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  134.         }
  135.         free(A); free(B); free(V);
  136.     }
  137.     else                                          /* 其他进程接收数据A,b和x的初值*/
  138.     {
  139.         MPI_Recv(a,m*size,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
  140.         MPI_Recv(b,m,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
  141.     }
  142.     time1=MPI_Wtime();                            /* 取开始计算时间 */
  143.     /* 求每行主对角线右边的元素和 */
  144.     for(i=0;i<m;i++)                              /* 开始计算 */
  145.     {
  146.         sum[i]=0.0;
  147.         for(j=0;j<size;j++)
  148.         {
  149.             if (j>(my_rank*m+i))
  150.                 sum[i]=sum[i]+a(i,j)*v(j);
  151.         }
  152.     }
  153.     while (total<size)                            /* total为迭代差小于E的x分量个数,当所有解即所有x分量的迭代差均小于E时,停止计算 */
  154.     {
  155.         iteration=0;                              /* iteration为各个process中迭代差小于E的x分量个数 */
  156.         total=0;
  157.         for(j=0;j<size;j++)                       /* 依次以第0,1, ..., n-1行为主行 */
  158.         {
  159.             r=j%m; q=j/m;
  160.             if (my_rank==q)                       /* 主行所在的处理器 */
  161.             {
  162.                 temp=v(my_rank*m+r);
  163.                 /*sum[j]累加上本进程中已算出的新的x项,这些新项的行号小于当前处理的行号*/
  164.                 for(i=0;i<r;i++)
  165.                     sum[r]=sum[r]+a(r,my_rank*m+i)*v(my_rank*m+i);
  166.                 /*计算新的x值*/
  167.                 v(my_rank*m+r)=(1-w)*temp+w*(b(r)-sum[r])/a(r,my_rank*m+r);
  168.                 if (fabs(v(my_rank*m+r)-temp)<E)  /*新旧值差小于E,累加iteration*/
  169.                     iteration++;
  170.                 /*广播新的x值*/
  171.                 MPI_Bcast(&v(my_rank*m+r),1,MPI_FLOAT,my_rank,MPI_COMM_WORLD);
  172.                 sum[r]=0.0;                       /*当前行处理完毕后,sum清零,为下一次计算做准备*/
  173.                 for(i=0;i<r;i++)                  /*本process中前面行的sum累加新x值项,为下一次计算累加对角元素左边的项*/
  174.                     sum[i]=sum[i]+a(i,my_rank*m+r)*v(my_rank*m+r);
  175.             }
  176.             else                                  /*非主行所在的处理器*/
  177.             {
  178.                                                   /*接收广播来的新x值*/
  179.                 MPI_Bcast(&v(q*m+r),1,MPI_FLOAT,q,MPI_COMM_WORLD);
  180.                 for(i=0;i<m;i++)                  /*累加新x值对应的项*/
  181.                     sum[i]=sum[i]+a(i,q*m+r)*v(q*m+r);
  182.             }
  183.         }
  184.         /* 求所有处理器中iteration值的和total并广播到所有处理器中 */
  185.         MPI_Allreduce(&iteration,&total,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
  186.         loop++;
  187.         if (my_rank==0)                           /* 打印每次迭代的total值 */
  188.             printf("%-2d th total vaule = %dn",loop,total);
  189.     }                                             /* while */
  190.     time2=MPI_Wtime();                            /* 取结束时间 */
  191.     if (my_rank==0)                               /* 打印解x值以及运行、分配和计算时间 */
  192.     {
  193.         for(i = 0; i < size; i ++) printf("x[%d] = %fn",i,v(i));
  194.         printf("n");
  195.         printf("Iteration num = %dn",loop);
  196.         printf("Whole running time    = %f secondsn",time2-starttime);
  197.         printf("Distribute data time  = %f secondsn",time1-starttime);
  198.         printf("Parallel compute time = %f secondsn",time2-time1);
  199.     }
  200.     MPI_Barrier(MPI_COMM_WORLD);                  /* 同步所有进程 */
  201.     MPI_Finalize();                               /* 结束计算 */
  202.     Environment_Finalize(a,b,v);
  203.     return (0);
  204. }