jacobi.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 v1(x) v1[x]
  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;                                            /* the size of work on eath processor */
  18. float *B;                                         /* store matrix B */
  19. float *A;                                         /* store matrix A */
  20. float *V;                                         /* store the value x */
  21. double starttime;
  22. double time1;
  23. double time2;
  24. int my_rank;                                      /* store the identifier of each processor */
  25. int p;                                            /* store the number of processor */
  26. MPI_Status status;
  27. FILE *fdA,*fdB;
  28. int i,j,group_size;
  29. float sum;
  30. float *b;
  31. float *v;
  32. float *a;
  33. float *v1;
  34. float lmax;                                       /* the differnce between v1 and v */
  35. float max=1.0;                                    /* max of lmax */
  36. int loop=0;                                       /* store the number of iteration */
  37. /*name   Environment_Finalize
  38.  *input  a,b,v,v1
  39.  *usage  free a,b,v,v1
  40.  */
  41. void Environment_Finalize(float *a,float *b,float *v,float *v1)
  42. {
  43.     free(a);
  44.     free(b);
  45.     free(v);
  46.     free(v1);
  47.     if(my_rank==0)
  48.     {
  49.         free(A);
  50.         free(B);
  51.         free(V);
  52.     }
  53. }
  54. /*name   Input
  55.  *usage  read size,N,A,B,V from input file on the root
  56.  */
  57. void Input()
  58. {
  59.     if(my_rank==0)
  60.     {
  61.         starttime=MPI_Wtime();
  62.         fdA=fopen("dataIn.txt","r");              /* Read from file "dataIn.txt" */
  63.         fscanf(fdA,"%d %d", &size, &N);
  64.         if (size != N-1)
  65.         {
  66.             printf("the input is wrongn");
  67.             exit(1);
  68.         }
  69.         A=(float *)malloc(floatsize*size*size);
  70.         B=(float *)malloc(floatsize*size);
  71.         V=(float *)malloc(floatsize*size);
  72.         for(i = 0; i < size; i++)
  73.         {
  74.             for(j = 0; j < size; j++)
  75.             {
  76.                 fscanf(fdA,"%f", A+i*size+j);
  77.             }
  78.             fscanf(fdA,"%f", B+i);
  79.         }
  80.         for(i = 0; i < size; i ++)
  81.         {
  82.             fscanf(fdA,"%f", V+i);
  83.         }
  84.         fclose(fdA);
  85.     }
  86. }
  87. /*name   nodInit
  88.  *usage  allocate spaces for v,a,b,v1 on each processor
  89.  *       and store v value on the root
  90.  */
  91. void nodInit()
  92. {
  93.     m=size/p;if (size%p!=0) m++;
  94.     v=(float *)malloc(floatsize*size);
  95.     a=(float *)malloc(floatsize*m*size);
  96.     b=(float *)malloc(floatsize*m);
  97.     v1=(float *)malloc(floatsize*m);              /* store v value */
  98.     if (a==NULL||b==NULL||v==NULL||v1==NULL)
  99.         printf("allocate space fail!");
  100.     if (my_rank==0)
  101.     {
  102.         for(i=0;i<size;i++)
  103.             v(i)=V(i);
  104.     }
  105. }
  106. /*name   Data_send
  107.  *usage  initialize the value of a,b on root and broadcast
  108.  */
  109. void Data_send()
  110. {
  111.     if (my_rank==0)
  112.     {
  113.         for(i=0;i<m;i++)
  114.             for(j=0;j<size;j++)
  115.                 a(i,j)=A(i,j);
  116.         for(i=0;i<m;i++)
  117.             b(i)=B(i);
  118.     }
  119.     /* root sends all other nodes the corresponding data */
  120.     if (my_rank==0)
  121.     {
  122.         for(i=1;i<p;i++)
  123.         {
  124.             MPI_Send(&(A(m*i,0)),m*size,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  125.             MPI_Send(&(B(m*i)),m,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  126.         }
  127.         /* free(A); free(B); free(V); */
  128.     }
  129.     else
  130.     {
  131.         MPI_Recv(a,m*size,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
  132.         MPI_Recv(b,m,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status);
  133.     }
  134. }
  135. /*name   Compute
  136.  *usage  compute v1(i) =( b(i)-∑a(i,j) * v(j) ) / a(i,my_rank * m + i)
  137.  *       lmax = abs(v1(i)-v(my_rank * m + i))
  138.  *      collect the max lmax on each processor,if it is greater than E
  139.  *      then continue else end while
  140.  */
  141. void Compute()
  142. {
  143.     time1=MPI_Wtime();
  144.     /* computing start */
  145.     while (max>E)                                 /* The precision requirement */
  146.     {
  147.         lmax=0.0;
  148.         for(i=0;i<m;i++)
  149.         {
  150.             if(my_rank*m+i<size)
  151.             {
  152.                 sum=0.0;
  153.                 for(j=0;j<size;j++)
  154.                     if (j!=(my_rank*m+i))
  155.                         sum=sum+a(i,j)*v(j);
  156.                                                   /* computes the new elements */
  157.                 v1(i)=(b(i) - sum)/a(i,my_rank * m + i);
  158.                 if (fabs(v1(i)-v(i))>lmax)
  159.                     lmax=fabs(v1(i)-v(my_rank * m + i));
  160.             }
  161.         }
  162.         /*Find the max element in the vector*/
  163.         MPI_Allreduce(&lmax,&max,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
  164.         /*Gather all the elements of the vector from all nodes*/
  165.         MPI_Allgather(v1,m,MPI_FLOAT,v,m,MPI_FLOAT,MPI_COMM_WORLD);
  166.         loop++;
  167.     }                                             /* while */
  168.     time2=MPI_Wtime();
  169. }
  170. /*name   Output
  171.  *usage  print messages
  172.  */
  173. void Output()
  174. {
  175.     if (my_rank==0)
  176.     {
  177.         printf("Input of file "dataIn.txt"n");
  178.         printf("%dt%dn", size, N);
  179.         for(i = 0; i < size; i ++)
  180.         {
  181.             for(j = 0; j < size; j ++) printf("%ft",A(i,j));
  182.             printf("%fn",B(i));
  183.         }
  184.         printf("n");
  185.         for(i = 0; i < size; i ++)
  186.         {
  187.             printf("%ft", V(i));
  188.         }
  189.         printf("nn");
  190.         printf("nOutput of solutionn");
  191.         for(i = 0; i < size; i ++) printf("x[%d] = %fn",i,v(i));
  192.         printf("n");
  193.         printf("Iteration num = %dn",loop);
  194.         printf("Whole running time    = %f secondsn",time2-starttime);
  195.         printf("Distribute data time  = %f secondsn",time1-starttime);
  196.         printf("Parallel compute time = %f secondsn",time2-time1);
  197.         fdB=fopen("dataOut.txt","w");
  198.         fprintf(fdB,"Input of file "dataIn.txt"n");
  199.         fprintf(fdB,"%dt%dn", size, N);
  200.         for(i = 0; i < size; i ++)
  201.         {
  202.             for(j = 0; j < size; j ++) fprintf(fdB,"%ft",A(i,j));
  203.             fprintf(fdB,"%fn",B(i));
  204.         }
  205.         fprintf(fdB,"n");
  206.         for(i = 0; i < size; i ++)
  207.         {
  208.             fprintf(fdB,"%ft", V(i));
  209.         }
  210.         fprintf(fdB,"nn");
  211.         fprintf(fdB,"nOutput of solutionn");
  212.         for(i = 0; i < size; i ++) fprintf(fdB,"x[%d] = %fn",i,v(i));
  213.         fprintf(fdB,"n");
  214.         fprintf(fdB,"Iteration num = %dn",loop);
  215.         fprintf(fdB,"Whole running time    = %f secondsn",time2-starttime);
  216.         fprintf(fdB,"Distribute data time  = %f secondsn",time1-starttime);
  217.         fprintf(fdB,"Parallel compute time = %f secondsn",time2-time1);
  218.     }
  219. }
  220. /*name   main
  221.  *input  argc:parameter number of command lines
  222.  *       argv:array of each command lines' parameter
  223.  *output return 0 if program is running correctly
  224.  */
  225. int main(int argc, char **argv)
  226. {
  227.     MPI_Init(&argc,&argv);
  228.     MPI_Comm_size(MPI_COMM_WORLD,&group_size);
  229.     MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
  230.     p=group_size;
  231.     Input();
  232.     /*Broadcast the size of  factor Matrix*/
  233.     MPI_Bcast(&size,1,MPI_INT,0,MPI_COMM_WORLD);
  234.     nodInit();
  235.     /*Broadcast the initial vector*/
  236.     MPI_Bcast(v,size,MPI_FLOAT,0,MPI_COMM_WORLD);
  237.     Data_send();
  238.     Compute();
  239.     Output();
  240.     MPI_Barrier(MPI_COMM_WORLD);
  241.     MPI_Finalize();
  242.     Environment_Finalize(a,b,v,v1);
  243.     return (0);
  244. }