开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 并行计算 > 查看源码
cjacobi.c
资源名称：mpisourceprogram.rar　[点击查看]
上传用户：jsljixie
上传日期：2013-08-15
资源大小：827k
文件大小：22k
源码类别：

并行计算

开发平台：

Visual C++

cjacobi.c：源码内容
  1. #include "stdio.h"
  2. #include "stdlib.h"
  3. #include "mpi.h"
  4. #include "math.h"
  5. #include "string.h"
  6. #define E 0.00001
  7. #define min -1
  8. #define intsize sizeof(int)
  9. #define charsize sizeof(char)
  10. #define floatsize sizeof(float)
  11. #define A(x,y) A[x*N+y]
  12. #define I(x,y) I[x*N+y]
  13. #define a(x,y) a[x*N+y]
  14. #define e(x,y) e[x*N+y]
  15. #define b(x) b[x]
  16. #define buffer(x,y) buffer[x*N+y]
  17. #define buffee(x,y) buffee[x*N+y]
  18. int M,N;
  19. int *b;
  20. int m,p;
  21. int myid,group_size;
  22. float *A,*I;
  23. float *a,*e;
  24. float max;
  25. float sum;
  26. double starttime,time1,time2,time3,endtime;
  27. MPI_Status status;
  29. float sgn(float v)
  30. {
  31.     float f;
  32.     if (v>=0) f=1;
  33.     else f=-1;
  34.     return f;
  35. }
  38. void read_fileA()
  39. {
  40.     int i,j;
  42.     FILE *fdA;
  43.     time1=MPI_Wtime();
  45.     fdA=fopen("dataIn.txt","r");
  46.     fscanf(fdA,"%d %d", &M, &N);
  48.     if(M != N)
  49.     {
  50.         puts("The input is error!");
  51.         exit(0);
  52.     }
  54.     m=N/p; if (N%p!=0) m++;
  55.     A=(float*)malloc(floatsize*N*m*p);
  56.     I=(float*)malloc(floatsize*N*m*p);
  58.     for(i = 0; i < M; i ++)
  59.     {
  60.         for(j = 0; j < M; j ++) fscanf(fdA, "%f", A+i*M+j);
  61.     }
  62.     fclose(fdA);
  64.     printf("Input of file "dataIn.txt"n");
  65.     printf("%dt %dn",M, N);
  66.     for(i=0;i<M;i++)
  67.     {
  68.         for(j=0;j<N;j++) printf("%ft",A(i,j));
  69.         printf("n");
  70.     }
  72.     for(i=0;i<N;i++)
  73.     {
  74.         for(j=0;j<N;j++)
  75.         {
  76.             if (i==j) I(i,j)=1;
  77.             else I(i,j)=0;
  78.         }
  79.     }
  80. }
  83. void send_a()
  84. {
  85.     int i;
  86.     for(i=1;i<p;i++)
  87.     {
  88.         MPI_Send(&(A(m*i,0)),m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  89.         MPI_Send(&(I(m*i,0)),m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD);
  90.     }
  91.     free(A);
  92.     free(I);
  93. }
  96. void get_a()
  97. {
  98.     int i,j;
  99.     if (myid==0)
  100.     {
  101.         for(i=0;i<m;i++)
  102.             for(j=0;j<N;j++)
  103.         {
  104.             a(i,j)=A(i,j);
  105.             e(i,j)=I(i,j);
  106.         }
  107.     }
  108.     else
  109.     {
  110.         MPI_Recv(a,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
  111.         MPI_Recv(e,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
  112.     }
  113. }
  116. int main(int argc,char **argv)
  117. {
  118.     float *c;
  119.     int k;
  120.     int loop;
  121.     int i,j,v,z,r,t,y;
  122.     int i1,j1;
  123.     float f,g,h;
  124.     float sin1,sin2,cos1;
  125.     float s1,c1;
  126.     float *br,*bt,*bi,*bj,*zi,*zj;
  127.     float *temp1,*temp2,*buffer,*buffee;
  128.     int counter,current;
  129.     int *buf;
  130.     int mml,mpl;
  131.     float bpp,bqq,bpq,bqp;
  132.     float lmax;
  134.     MPI_Init(&argc,&argv);
  135.     MPI_Comm_size(MPI_COMM_WORLD,&group_size);
  136.     MPI_Comm_rank(MPI_COMM_WORLD,&myid);
  137.     if (myid==0) starttime=MPI_Wtime();
  138.     p=group_size;
  140.     max=10.0;
  141.     loop=0;
  143.     if (myid==0) read_fileA();
  145.     MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);
  147.     m=N/p;
  148.     if (N%p!=0) m++;
  150.     a=(float*)malloc(floatsize*m*N);
  151.     e=(float*)malloc(floatsize*m*N);
  152.     b=(int*)malloc(intsize*m);
  153.     br=(float*)malloc(floatsize*N);
  154.     bt=(float*)malloc(floatsize*N);
  155.     bi=(float*)malloc(floatsize*m);
  156.     bj=(float*)malloc(floatsize*m);
  157.     zi=(float*)malloc(floatsize*m);
  158.     zj=(float*)malloc(floatsize*m);
  159.     temp1=(float*)malloc(floatsize*4);
  160.     temp2=(float*)malloc(floatsize*4*p);
  162.     if ((myid==p-1)&&(myid%2!=0))
  163.     {
  164.         buffer=(float*)malloc(floatsize*m/2*N);
  165.         buffee=(float*)malloc(floatsize*m/2*N);
  166.         buf=(int*)malloc(intsize*m/2);            /* store i,j */
  167.     }
  169.     if ((myid%2!=0)&&(myid!=p-1))
  170.     {
  171.         buffer=(float*)malloc(floatsize*m*N);
  172.         buffee=(float*)malloc(floatsize*m*N);
  173.         buf=(int*)malloc(intsize*m);
  174.     }
  176.     if (myid==0)
  177.     {
  178.         get_a();
  179.         send_a();
  180.     }
  181.     else
  182.     {
  183.         get_a();
  184.     }
  186.     if (myid==0)                                  /* start computing now */
  187.         time2=MPI_Wtime();
  189.     for(i=0;i<m;i++)
  190.         b(i)=myid*m+i;
  192.     while (fabs(max)>E)
  193.     {
  194.         loop++;
  195.         for(i=myid*m;i<(myid+1)*m-1;i++)
  196.             for(j=i+1;j<=(myid+1)*m-1;j++)
  197.         {
  198.             r=i%m; t=j%m;                         /* to make a(r,j) zero use a(r,i) a(t,j) */
  199.             if(a(r,j)!=0)
  200.             {
  201.                 f=(-a(r,j));
  202.                 g=(a(t,j)-a(r,i))/2;
  203.                 h=sgn(g)*f/sqrt(f*f+g*g);
  204.                 sin2=h;
  205.                 sin1=h/sqrt(2*(1+sqrt(1-h*h)));
  206.                 cos1=sqrt(1-sin1*sin1);
  208.                 bpp=a(r,i)*cos1*cos1+a(t,j)*sin1*sin1+a(r,j)*sin2;
  209.                 bqq=a(r,i)*sin1*sin1+a(t,j)*cos1*cos1-a(r,j)*sin2;
  210.                 bpq=0; bqp=0;
  212.                 for(v=0;v<N;v++)                  /* compute row */
  213.                     if ((v!=i)&&(v!=j))
  214.                 {
  215.                     br[v]=a(r,v)*cos1+a(t,v)*sin1;
  216.                     a(t,v)=(-a(r,v))*sin1+a(t,v)*cos1;
  217.                 }
  219.                 for(v=0;v<N;v++)                  /* row */
  220.                     if ((v!=i)&&(v!=j))
  221.                 {
  222.                     a(r,v)=br[v];
  223.                 }
  225.                 for(v=0;v<m;v++)
  226.                     br[v]=e(v,i)*cos1+e(v,j)*sin1;
  228.                 for(v=0;v<m;v++)
  229.                     e(v,j)=e(v,i)*(-sin1)+e(v,j)*cos1;
  231.                 for(v=0;v<m;v++)
  232.                     e(v,i)=br[v];
  234.                 for(v=0;v<m;v++)                  /* compute col */
  235.                     if ((v!=r)&&(v!=t))
  236.                 {
  237.                     bi[v]=a(v,i)*cos1+a(v,j)*sin1;
  238.                     a(v,j)=(-a(v,i))*sin1+a(v,j)*cos1;
  239.                 }
  241.                 for(v=0;v<m;v++)                  /* col */
  242.                     if ((v!=r)&&(v!=t))
  243.                         a(v,i)=bi[v];
  245.                 a(r,i)=bpp;
  246.                 a(t,j)=bqq;
  247.                 a(r,j)=bpq;
  248.                 a(t,i)=bqp;
  250.                 temp1[0]=sin1;
  251.                 temp1[1]=cos1;
  252.                 temp1[2]=(float)i;
  253.                 temp1[3]=(float)j;
  254.             }
  255.             else
  256.             {
  257.                 temp1[0]=0.0;
  258.                 temp1[1]=0.0;
  259.                 temp1[2]=0.0;
  260.                 temp1[3]=0.0;
  261.             }                                     /* end of if */
  263.             MPI_Allgather(temp1,4,MPI_FLOAT,temp2,4,MPI_FLOAT,MPI_COMM_WORLD);
  264.             current=0;
  266.             for(v=1;v<=p;v++)
  267.             {
  268.                 s1=temp2[(v-1)*4+0];
  269.                 c1=temp2[(v-1)*4+1];
  270.                 i1=(int)temp2[(v-1)*4+2];
  271.                 j1=(int)temp2[(v-1)*4+3];
  273.                 if ((s1!=0.0)||(c1!=0.0)||(i1!=0)||(j1!=0))
  274.                 {
  275.                     if (myid!=current)
  276.                     {
  278.                         for(z=0;z<m;z++)
  279.                         {
  280.                             zi[z]=a(z,i1)*c1 + a(z,j1)*s1;
  281.                             a(z,j1)=-a(z,i1)*s1 + a(z,j1)*c1;
  282.                         }
  284.                         for(z=0;z<m;z++)
  285.                             a(z,i1)=zi[z];
  287.                         for(z=0;z<m;z++)
  288.                             zi[z]=e(z,i1)*c1+e(z,j1)*s1;
  290.                         for(z=0;z<m;z++)
  291.                             e(z,j1)=-e(z,i1)*s1+e(z,j1)*c1;
  293.                         for(z=0;z<m;z++)
  294.                             e(z,i1)=zi[z];
  296.                     }                             /* if myid!=current */
  298.                 }                                 /* if */
  300.                 current=current+1;
  301.             }                                     /* for v */
  302.         }                                         /* for i,j */
  304.         for(counter=1;counter<=2*p-2;counter++)
  305.         {
  307.             if (myid==0)
  308.             {
  309.                 MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  310.                 MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  311.                 MPI_Send(&b(m/2),m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);
  313.                 MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  314.                 MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  315.                 MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
  316.             }
  318.             if ((myid==p-1)&&(myid%2!=0))
  319.             {
  320.                 for(i=m/2;i<m;i++)
  321.                     for(j=0;j<N;j++)
  322.                         buffer((i-m/2),j)=a(i,j);
  324.                 for(i=m/2;i<m;i++)
  325.                     for(j=0;j<N;j++)
  326.                         buffee((i-m/2),j)=e(i,j);
  328.                 for(i=m/2;i<m;i++)
  329.                     buf[i-m/2]=b(i);
  331.                 for(i=0;i<m/2;i++)
  332.                     for(j=0;j<N;j++)
  333.                         a((i+m/2),j)=a(i,j);
  335.                 for(i=0;i<m/2;i++)
  336.                     for(j=0;j<N;j++)
  337.                         e((i+m/2),j)=e(i,j);
  339.                 for(i=0;i<m/2;i++)
  340.                     b(m/2+i)=b(i);
  342.                 MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  343.                 MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  344.                 MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  345.                 MPI_Send(buffer,m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  346.                 MPI_Send(buffee,m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  347.                 MPI_Send(buf,m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  348.             }
  350.             if ((myid==p-1)&&(myid%2==0))
  351.             {
  352.                 MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  353.                 MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  354.                 MPI_Send(&b(m/2),m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  356.                 for(i=0;i<m/2;i++)
  357.                     for(j=0;j<N;j++)
  358.                         a((i+m/2),j)=a(i,j);
  360.                 for(i=0;i<m/2;i++)
  361.                     for(j=0;j<N;j++)
  362.                         e((i+m/2),j)=e(i,j);
  364.                 for(i=0;i<m/2;i++)
  365.                     b(i+m/2)=b(i);
  367.                 MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  368.                 MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  369.                 MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  370.             }
  372.             if ((myid!=0)&&(myid!=p-1))
  373.             {
  374.                 if(myid%2==0)
  375.                 {
  376.                     MPI_Send(&(a(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  377.                     MPI_Send(&(e(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  378.                     MPI_Send(&b(0),m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);
  379.                     MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  380.                     MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  381.                     MPI_Send(&b(m/2),m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  383.                     MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  384.                     MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  385.                     MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  386.                     MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  387.                     MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  388.                     MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
  389.                 }
  391.                 if(myid%2!=0)
  392.                 {
  393.                     for(i=0;i<m;i++)
  394.                         for(j=0;j<N;j++)
  395.                             buffer(i,j)=a(i,j);
  397.                     for(i=0;i<m;i++)
  398.                         for(j=0;j<N;j++)
  399.                             buffee(i,j)=e(i,j);
  401.                     for(i=0;i<m;i++)
  402.                         buf[i]=b(i);
  404.                     MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  405.                     MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  406.                     MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  407.                     MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  408.                     MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  409.                     MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
  411.                     MPI_Send(&(buffer(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  412.                     MPI_Send(&(buffee(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  413.                     MPI_Send(&buf[0],m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);
  414.                     MPI_Send(&(buffer(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  415.                     MPI_Send(&(buffee(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  416.                     MPI_Send(&buf[m/2],m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  417.                 }
  418.             }
  420.             for(i=0;i<m/2;i++)
  421.                 for(j=m/2;j<m;j++)
  422.             {
  423.                 if (a(i,b(j))!=0)
  424.                 {
  425.                     f=-a(i,b(j));
  426.                     g=(a(j,b(j))-a(i,b(i)))/2;
  427.                     h=sgn(g)*f/sqrt(f*f+g*g);
  428.                     sin2=h;
  429.                     sin1=h/sqrt(2*(1+sqrt(1-h*h)));
  430.                     cos1=sqrt(1-sin1*sin1);
  432.                     bpp=a(i,b(i))*cos1*cos1+a(j,b(j))*sin1*sin1+a(i,b(j))*sin2;
  433.                     bqq=a(i,b(i))*sin1*sin1+a(j,b(j))*cos1*cos1-a(i,b(j))*sin2;
  434.                     bpq=0; bqp=0;
  436.                     for(v=0;v<N;v++)              /* compute row */
  437.                         if ((v!=b(i))&&(v!=b(j)))
  438.                     {
  439.                         br[v]=a(i,v)*cos1+a(j,v)*sin1;
  440.                         a(j,v)=-a(i,v)*sin1+a(j,v)*cos1;
  441.                     }
  443.                     for(v=0;v<N;v++)              /* row */
  444.                         if ((v!=b(i))&&(v!=b(j)))
  445.                             a(i,v)=br[v];
  447.                     for(v=0;v<m;v++)
  448.                         br[v]=e(v,b(i))*cos1+e(v,b(j))*sin1;
  450.                     for(v=0;v<m;v++)
  451.                         e(v,b(j))=e(v,b(i))*(-sin1)+e(v,b(j))*cos1;
  453.                     for(v=0;v<m;v++)
  454.                         e(v,b(i))=br[v];
  456.                     for(v=0;v<m;v++)              /* compute col */
  457.                         if ((v!=i)&&(v!=j))
  458.                     {
  459.                         bi[v]=a(v,b(i))*cos1+a(v,b(j))*sin1;
  460.                         a(v,b(j))=-a(v,b(i))*sin1+a(v,b(j))*cos1;
  461.                     }
  463.                     for(v=0;v<m;v++)              /* col */
  464.                         if ((v!=i)&&(v!=j))
  465.                             a(v,b(i))=bi[v];
  467.                     a(i,b(i))=bpp;
  468.                     a(j,b(j))=bqq;
  469.                     a(i,b(j))=bpq;                /* 0 */
  470.                     a(j,b(i))=bqp;                /* 0 */
  472.                     temp1[0]=sin1;
  473.                     temp1[1]=cos1;
  474.                     temp1[2]=(float)b(i);
  475.                     temp1[3]=(float)b(j);
  476.                 }
  478.                 else
  479.                 {
  480.                     temp1[0]=0.0;
  481.                     temp1[1]=0.0;
  482.                     temp1[2]=0.0;
  483.                     temp1[3]=0.0;
  484.                 }
  486.                 MPI_Allgather(temp1,4,MPI_FLOAT,temp2,4,MPI_FLOAT,MPI_COMM_WORLD);
  487.                 current=0;
  489.                 for(v=1;v<=p;v++)
  490.                 {
  491.                     s1=temp2[(v-1)*4+0];
  492.                     c1=temp2[(v-1)*4+1];
  493.                     i1=(int)temp2[(v-1)*4+2];
  494.                     j1=(int)temp2[(v-1)*4+3];
  496.                     if ((s1!=0.0)||(c1!=0.0)||(i1!=0)||(j1!=0))
  497.                     {
  498.                         if (myid!=current)
  499.                         {
  500.                             for(z=0;z<m;z++)
  501.                             {
  502.                                 zi[z]=a(z,i1)*c1 + a(z,j1)*s1;
  503.                                 a(z,j1)=-a(z,i1)*s1 + a(z,j1)*c1;
  504.                             }
  505.                             for(z=0;z<m;z++)
  506.                                 a(z,i1)=zi[z];
  508.                             for(z=0;z<m;z++)
  509.                                 zi[z]=e(z,i1)*c1+e(z,j1)*s1;
  511.                             for(z=0;z<m;z++)
  512.                                 e(z,j1)=-e(z,i1)*s1+e(z,j1)*c1;
  514.                             for(z=0;z<m;z++)
  515.                                 e(z,i1)=zi[z];
  517.                         }                         /* if myid!=current */
  518.                     }                             /* if */
  519.                     current=current+1;
  520.                 }                                 /* for v */
  521.             }                                     /* for i,j */
  523.         }                                         /*  counter */
  525.         if (myid==0)
  526.         {
  527.             MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  529.             MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  530.             MPI_Send(&b(m/2),m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);
  532.             MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  533.             MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  534.             MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
  535.         }
  537.         if ((myid==p-1)&&(myid%2!=0))
  538.         {
  539.             for(i=m/2;i<m;i++)
  540.                 for(j=0;j<N;j++)
  541.                     buffer((i-m/2),j)=a(i,j);
  543.             for(i=m/2;i<m;i++)
  544.                 for(j=0;j<N;j++)
  545.                     buffee((i-m/2),j)=e(i,j);
  547.             for(i=m/2;i<m;i++)
  548.                 buf[i-m/2]=b(i);
  550.             for(i=0;i<m/2;i++)
  551.                 for(j=0;j<N;j++)
  552.                     a((i+m/2),j)=a(i,j);
  554.             for(i=0;i<m/2;i++)
  555.                 for(j=0;j<N;j++)
  556.                     e((i+m/2),j)=e(i,j);
  558.             for(i=0;i<m/2;i++)
  559.                 b(m/2+i)=b(i);
  561.             MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  562.             MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  563.             MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  564.             MPI_Send(buffer,m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  565.             MPI_Send(buffee,m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  566.             MPI_Send(buf,m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  567.         }
  569.         if ((myid==p-1)&&(myid%2==0))
  570.         {
  571.             MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  572.             MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  573.             MPI_Send(&b(m/2),m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  575.             for(i=0;i<m/2;i++)
  576.                 for(j=0;j<N;j++)
  577.                     a((i+m/2),j)=a(i,j);
  579.             for(i=0;i<m/2;i++)
  580.                 for(j=0;j<N;j++)
  581.                     e((i+m/2),j)=e(i,j);
  583.             for(i=0;i<m/2;i++)
  584.                 b(i+m/2)=b(i);
  586.             MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  587.             MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  588.             MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  589.         }
  591.         if ((myid!=0)&&(myid!=p-1))
  592.         {
  593.             if(myid%2==0)
  594.             {
  595.                 MPI_Send(&(a(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  596.                 MPI_Send(&(e(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  597.                 MPI_Send(&b(0),m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);
  598.                 MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  599.                 MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  600.                 MPI_Send(&b(m/2),m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  602.                 MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  603.                 MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  604.                 MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  605.                 MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  606.                 MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  607.                 MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
  608.             }
  610.             if(myid%2!=0)
  611.             {
  612.                 for(i=0;i<m;i++)
  613.                     for(j=0;j<N;j++)
  614.                         buffer(i,j)=a(i,j);
  616.                 for(i=0;i<m;i++)
  617.                     for(j=0;j<N;j++)
  618.                         buffee(i,j)=e(i,j);
  620.                 for(i=0;i<m;i++)
  621.                     buf[i]=b(i);
  623.                 MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  624.                 MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
  625.                 MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
  626.                 MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  627.                 MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
  628.                 MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
  630.                 MPI_Send(&(buffer(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  631.                 MPI_Send(&(buffee(0,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
  632.                 MPI_Send(&buf[0],m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);
  633.                 MPI_Send(&(buffer(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  634.                 MPI_Send(&(buffee(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
  635.                 MPI_Send(&buf[m/2],m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
  636.             }
  637.         }
  639.         lmax=min;
  641.         for(i=0;i<m;i++)
  642.             for(j=0;j<N;j++)
  643.                 if ((m*myid+i)!=j)
  644.                 {
  645.                     if (fabs(a(i,j))>lmax)
  646.                         lmax=fabs(a(i,j));
  647.                 }
  649.         MPI_Allreduce(&lmax,&max,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
  651.     }                                             /* while */
  653.     if (myid==0)
  654.     {
  655.         time3=MPI_Wtime();
  656.         A=(float*)malloc(floatsize*N*m*p);
  657.         I=(float*)malloc(floatsize*N*m*p);
  659.         for(i=0;i<m;i++)
  660.             for(j=0;j<N;j++)
  661.         {
  662.             A(i,j)=a(i,j);
  663.             I(i,j)=e(i,j);
  664.         }
  665.     }
  667.     if (myid!=0)
  668.     {
  669.         MPI_Send(a,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD);
  670.         MPI_Send(e,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD);
  671.     }
  672.     else
  673.         for(i=1;i<p;i++)
  674.     {
  675.         MPI_Recv(a,m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD,&status);
  676.         MPI_Recv(e,m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD,&status);
  678.         for(j=0;j<m;j++)
  679.             for(k=0;k<N;k++)
  680.                 A((i*m+j),k)=a(j,k);
  682.         for(j=0;j<m;j++)
  683.             for(k=0;k<N;k++)
  684.                 I((i*m+j),k)=e(j,k);
  685.     }
  687.     if (myid==0)
  688.     {
  689.         for(i=0;i<N;i++)
  690.             printf("the %dst envalue:%fn",i,A(i,i));
  692.         endtime=MPI_Wtime();
  694.         printf("n");
  695.         printf("Iteration num = %dn",loop);
  696.         printf("Whole running time    = %f secondsn",endtime-starttime);
  697.         printf("Distribute data time  = %f secondsn",time2-time1);
  698.         printf("Parallel compute time = %f secondsn",time3-time2);
  700.     }
  702.     MPI_Finalize();
  703.     free(a);
  704.     free(b);
  705.     free(c);
  706.     free(br);
  707.     free(bt);
  708.     free(bi);
  709.     free(bj);
  710.     free(zi);
  711.     free(zj);
  712.     free(buffer);
  713.     free(buf);
  714.     free(buffee);
  715.     free(A);
  716.     free(I);
  717.     free(temp1);
  718.     free(temp2);
  719.     return(0);
  720. }