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

并行计算

开发平台:

Unix_Linux

  1. #include "stdio.h"
  2. #include "stdlib.h"
  3. #include "mpi.h"
  4. #include "math.h"
  5. #include "string.h"
  6. #define E 0.0001
  7. #define intsize sizeof(int)
  8. #define floatsize sizeof(float)
  9. #define A(x,y) A[x*col+y]
  10. #define I(x,y) I[x*col+y]
  11. #define B(x) B[x]
  12. #define a(x,y) a[x*col+y]
  13. #define e(x,y) e[x*col+y]
  14. int col, N;
  15. int m,p;
  16. int myid,group_size;
  17. float *A,*I,*B;
  18. float *a,*e;
  19. MPI_Status status;
  20. float starttime,endtime,time1,time2;
  21. /*读输入的矩阵A*/
  22. void read_fileA()
  23. {
  24.     int i,j;
  25.     FILE *fdA;
  26.     time1=MPI_Wtime();
  27.     fdA=fopen("dataIn.txt","r");
  28.     fscanf(fdA,"%d %d", &col, &N);
  29.     if (col != N)
  30.     {
  31.         printf("the input is wrongn");
  32.         exit(1);
  33.     }
  34.     A=(float *)malloc(floatsize*col*col);
  35.     I=(float *)malloc(floatsize*col*col);
  36.     for(i = 0; i < col; i++)
  37.     {
  38.         for(j = 0; j < col; j++)
  39.         {
  40.             fscanf(fdA,"%f", A+i*col+j);
  41.         }
  42.     }
  43.     fclose(fdA);
  44.     printf("Input of file "dataIn.txt"n");
  45.     printf("%dt%dn", col, N);
  46.     for(i = 0; i < col; i ++)
  47.     {
  48.         for(j = 0; j < col; j ++) printf("%ft",A(i,j));
  49.         printf("n");
  50.     }
  51.     for(i=0;i<col;i++)
  52.     {
  53.         for(j=0;j<col;j++)
  54.         {
  55.             if (i==j) I(i,j)=1.0;
  56.             else I(i,j)=0.0;
  57.         }
  58.     }
  59. }
  60. int main(int argc,char **argv)
  61. {
  62.     int loop;
  63.     int i,j,v;
  64.     int p,group_size,myid;
  65.     int k;
  66.     float aa,bb,rr,c,s;
  67.     float su;
  68.     float *ss;
  69.     float *sum;
  70.     float *temp,*temp1;
  71.     MPI_Init(&argc,&argv);
  72.     MPI_Comm_size(MPI_COMM_WORLD,&group_size);
  73.     MPI_Comm_rank(MPI_COMM_WORLD,&myid);
  74.     if (myid==0) starttime=MPI_Wtime();
  75.     p=group_size;
  76.     loop=0;
  77.     k=0;
  78.     if (myid==0)
  79.         read_fileA();
  80.     MPI_Bcast(&col,1,MPI_INT,0,MPI_COMM_WORLD);
  81.     m=col/p; if (col%p!=0) m++;
  82.     a=(float*)malloc(floatsize*m*col);
  83.     e=(float*)malloc(floatsize*m*col);
  84.     temp=(float*)malloc(floatsize*m);
  85.     temp1=(float*)malloc(floatsize*m);
  86.     ss=(float*)malloc(floatsize*3);
  87.     sum=(float*)malloc(floatsize*3);
  88.     if (myid==0)
  89.         B=(float*)malloc(floatsize*col);
  90.     if (myid==0)
  91.     {
  92.         for(i=0;i<m;i++)
  93.             for(j=0;j<col;j++)
  94.                 a(i,j)=A(i,j);
  95.         for(i=0;i<m;i++)
  96.             for(j=0;j<col;j++)
  97.                 e(i,j)=I(i,j);
  98.         for(i=1;i<p;i++)
  99.         {
  100.             MPI_Send(&(A(m*i,0)),m*col,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  101.             MPI_Send(&(I(m*i,0)),m*col,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  102.         }
  103.     }
  104.     else
  105.     {
  106.         MPI_Recv(a,m*col,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
  107.         MPI_Recv(e,m*col,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
  108.     }
  109.     if (myid==0)                                  /*开始计时*/
  110.         time2=MPI_Wtime();
  111.     while (k<=col*(col-1)/2)
  112.     {
  113.         for(i=0;i<col;i++)
  114.             for(j=i+1;j<col;j++)
  115.         {
  116.             ss[0]=0; ss[1]=0; ss[2]=0;
  117.             sum[0]=0.0; sum[1]=0.0; sum[2]=0.0;
  118.             for(v=0;v<m;v++)
  119.                 ss[0]=ss[0]+a(v,i)*a(v,j);        /* 计算本进程所存的列子向量a[*, i], a[*, j]的内积赋值给ss[0] */
  120.             for(v=0;v<m;v++)
  121.                 ss[1]=ss[1]+a(v,i)*a(v,i);        /* 计算本进程所存的列子向量a[*, i], a[*, i]的内积赋值给ss[1] */
  122.             for(v=0;v<m;v++)
  123.                 ss[2]=ss[2]+a(v,j)*a(v,j);        /* 计算本进程所存的列子向量a[*, j], a[*, j]的内积赋值给ss[2] */
  124.             /* 计算所有进程中ss[i]的和赋给sum [i] (i = 0, 1, 2) */
  125.             MPI_Allreduce(ss,sum,3,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
  126.             if (fabs(sum[0])>E)                   /* 各个进程并行进行对i和j列正交化 */
  127.             {
  128.                 aa=2*sum[0];
  129.                 bb=sum[1]-sum[2];
  130.                 rr=sqrt(aa*aa+bb*bb);
  131.                 if (bb>=0)
  132.                 {
  133.                     c=sqrt((bb+rr)/(2*rr));
  134.                     s=aa/(2*rr*c);
  135.                 }
  136.                 if (bb<0)
  137.                 {
  138.                     s=sqrt((rr-bb)/(2*rr));
  139.                     c=aa/(2*rr*s);
  140.                 }
  141.                 for(v=0;v<m;v++)
  142.                 {
  143.                     temp[v]=c*a(v,i)+s*a(v,j);
  144.                     a(v,j)=(-s)*a(v,i)+c*a(v,j);
  145.                     temp1[v]=c*e(v,i)+s*e(v,j);
  146.                     e(v,j)=(-s)*e(v,i)+c*e(v,j);
  147.                 }
  148.                 for(v=0;v<m;v++)
  149.                 {
  150.                     a(v,i)=temp[v];
  151.                     e(v,i)=temp1[v];
  152.                 }
  153.             }
  154.             else
  155.                 k++;
  156.         }                                         /* for */
  157.         loop ++;
  158.     }                                             /* while */
  159.     if (myid==0)                                  /*0号进程从其余各进程中接收子矩阵a和e,得到经过变换的矩阵A和I*/
  160.     {
  161.         for(i=0;i<m;i++)
  162.             for(j=0;j<col;j++)
  163.                 A(i,j)=a(i,j);
  164.         for(i=0;i<m;i++)
  165.             for(j=0;j<col;j++)
  166.                 I(i,j)=e(i,j);
  167.     }
  168.     if (myid!=0)
  169.         MPI_Send(a,m*col,MPI_FLOAT,0,myid,MPI_COMM_WORLD);
  170.     else
  171.     {
  172.         for(j=1;j<p;j++)
  173.         {
  174.             MPI_Recv(a,m*col,MPI_FLOAT,j,j,MPI_COMM_WORLD,&status);
  175.             for(i=0;i<m;i++)
  176.                 for(k=0;k<col;k++)
  177.                     A((j*m+i),k)=a(i,k);
  178.         }
  179.     }
  180.     if (myid!=0)
  181.         MPI_Send(e,m*col,MPI_FLOAT,0,myid,MPI_COMM_WORLD);
  182.     else
  183.     {
  184.         for(j=1;j<p;j++)
  185.         {
  186.             MPI_Recv(e,m*col,MPI_FLOAT,j,j,MPI_COMM_WORLD,&status);
  187.             for(i=0;i<m;i++)
  188.                 for(k=0;k<col;k++)
  189.                     I((j*m+i),k)=e(i,k);
  190.         }
  191.     }
  192.     if (myid==0)                                  /*0号进程负责计算特征值*/
  193.     {
  194.         for(j=0;j<col;j++)
  195.         {
  196.             su=0.0;
  197.             for(i=0;i<col;i++)
  198.                 su=su+A(i,j)*A(i,j);
  199.             B(j)=sqrt(su);
  200.         }
  201.         printf("n");
  202.         for(j=0;j<col;j++)
  203.         {
  204.             if (I(0,j)*A(0,j)>0)                  /*判断矩阵A特征值的符号*/
  205.                 printf("the envalue of the matrix is %fn",B(j));
  206.             if (I(0,j)*A(0,j)<0)
  207.                 printf("the envalue of the matrix is %fn",(-1)*B(j));
  208.             if (I(0,j)*A(0,j)==0)
  209.             {
  210.                 if (I(1,j)*A(1,j)>=0)             /*判断矩阵A特征值的符号*/
  211.                     printf("the envalue of the matrix is %fn",B(j));
  212.                 if (I(1,j)*A(1,j)<0)
  213.                     printf("the envalue of the matrix is %fn",(-1)*B(j));
  214.             }
  215.         }
  216.     }
  217.     if (myid==0)
  218.     {
  219.         endtime=MPI_Wtime();
  220.         printf("n");
  221.         printf("Iteration num = %dn",loop);
  222.         printf("Whole running time    = %f sn",time2-starttime); /*输出整个运行时间*/
  223.         printf("Distribute data time  = %f sn",time1-starttime); /*输出划分时间*/
  224.         printf("Parallel compute time = %f sn",time2-time1);     /*输出并行计算时间*/
  225.     }
  226.     MPI_Finalize();
  227.     free(a);
  228.     free(e);
  229.     free(ss);
  230.     free(sum);
  231.     free(A);
  232.     free(I);
  233.     free(B);
  234.     free(temp);
  235.     free(temp1);
  236.     return(0);
  237. }