seidel.c
资源名称:mpi.zip [点击查看]
上传用户:hhyinxing
上传日期:2013-09-10
资源大小:833k
文件大小:6k
源码类别:
并行计算
开发平台:
Unix_Linux
- #include "stdio.h"
- #include "stdlib.h"
- #include "mpi.h"
- #include "math.h"
- #define E 0.0001
- #define a(x,y) a[x*size+y]
- #define b(x) b[x]
- #define v(x) v[x]
- /*A为size*size的系数矩阵*/
- #define A(x,y) A[x*size+y]
- #define B(x) B[x]
- #define V(x) V[x]
- #define intsize sizeof(int)
- #define floatsize sizeof(float)
- #define charsize sizeof(char)
- int size,N;
- int m;
- float *B;
- float *A;
- float *V;
- int my_rank;
- int p;
- MPI_Status status;
- FILE *fdA,*fdB,*fdB1;
- int i,j,my_rank,group_size;
- int k;
- float *sum;
- float *b;
- float *v;
- float *a;
- float *differ;
- float temp;
- int iteration,total,loop;
- int r,q;
- /*
- * 函数名:Environment_Finalize
- * 功能:释放各个进程的a,b,v
- */
- void Environment_Finalize(float *a,float *b,float *v)
- {
- free(a);
- free(b);
- free(v);
- }
- /*
- * 函数名: read_data()
- * 功能:读取数据
- */
- void read_data()
- {
- if(my_rank==0)
- {
- /*打开数据文件*/
- fdA=fopen("dataIn.txt","r");
- fscanf(fdA,"%d %d", &size, &N);
- if (size != N-1)
- {
- printf("the input is wrongn");
- exit(1);
- }
- /*分配空间*/
- A=(float *)malloc(floatsize*size*size);
- B=(float *)malloc(floatsize*size);
- V=(float *)malloc(floatsize*size);
- /*初始化矩阵A和向量B*/
- for(i = 0; i < size; i++)
- {
- for(j = 0; j < size; j++)
- fscanf(fdA,"%f", A+i*size+j);
- fscanf(fdA,"%f", B+i);
- }
- /*初始化向量V*/
- for(i = 0; i < size; i ++)
- fscanf(fdA,"%f", V+i);
- /*关闭数据文件*/
- fclose(fdA);
- /*输出数据文件信息*/
- printf("Input of file "dataIn.txt"n");
- printf("%dt%dn", size, N);
- /*输出矩阵A和向量B的值*/
- for(i = 0; i < size; i ++)
- {
- for(j = 0; j < size; j ++)
- printf("%ft",A(i,j));
- printf("%fn",B(i));
- }
- printf("n");
- /*输出解向量V的值*/
- for(i = 0; i < size; i ++)
- printf("%ft", V(i));
- printf("nn");
- printf("nOutput of resultn");
- }
- }
- /*
- * 函数名: broadcast_data()
- * 功能:广播数据
- */
- void broadcast_data()
- {
- /*0号进程将size广播给所有进程*/
- MPI_Bcast(&size,1,MPI_INT,0,MPI_COMM_WORLD);
- m=size/p;
- if (size%p!=0) m++;
- /*为每个进程分配需要的空间*/
- v=(float *)malloc(floatsize*size);
- a=(float *)malloc(floatsize*m*size);
- b=(float *)malloc(floatsize*m);
- sum=(float *)malloc(floatsize*m);
- /*判断是否分配成功*/
- if (a==NULL||b==NULL||v==NULL)
- {
- printf("allocate space fail!");
- exit(1);
- }
- /*如果是0号进程,就将本进程的解向量初始化为问题的解向量*/
- if (my_rank==0)
- {
- for(i=0;i<size;i++)
- v(i)=V(i);
- }
- /*初始解向量v被广播给所有进程*/
- MPI_Bcast(v,size,MPI_FLOAT,0,MPI_COMM_WORLD);
- /*0号进程采用行块划分将矩阵A划分为大小为m*size的p块子矩阵,
- 将B划分为大小为m的p块子向量,依次发送给1至p-1号进程*/
- if (my_rank==0)
- {
- /*初始化0号进程的a和b*/
- for(i=0;i<m;i++)
- for(j=0;j<size;j++)
- a(i,j)=A(i,j);
- for(i=0;i<m;i++)
- b(i)=B(i);
- /*将各自的a和b发送到1至p-1号进程*/
- for(i=1;i<p;i++)
- {
- MPI_Send(&(A(m*i,0)), m*size,MPI_FLOAT,i,i,MPI_COMM_WORLD);
- MPI_Send(&(B(m*i)),m, MPI_FLOAT,i,i,MPI_COMM_WORLD);
- }
- /*释放空间*/
- free(A); free(B); free(V);
- }
- else
- {
- /*如果不是0号进程,就接收从0号进程传来的子矩阵a和子向量b*/
- MPI_Recv(a,m*size,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD, &status);
- MPI_Recv(b,m,MPI_FLOAT,0,my_rank, MPI_COMM_WORLD,&status);
- }
- }
- /*
- * 函数名: cal_cross_data()
- * 功能:主对角元素右边的数据求和
- */
- void cal_cross_data()
- {
- /*所有进程并行地对主对角元素右边的数据求和*/
- for(i=0;i<m;i++)
- {
- sum[i]=0.0;
- for(j=0;j<size;j++)
- if (j>(my_rank*m+i))
- sum[i]=sum[i]+a(i,j)*v(j);
- }
- }
- /*
- * 函数名: die_dai_cal()
- * 功能:各个进程并行迭代求解
- */
- void die_dai_cal()
- {
- /*total为新旧值之差小于E的x的分量的个数*/
- while (total<size)
- {
- iteration=0;
- /*iteration为本进程中新旧值之差小于E的x的分量个数*/
- total=0;
- for(j=0;j<size;j++)
- {
- r=j%m;
- q=j/m;
- /*编号为q的进程负责计算解向量并广播给所有进程*/
- if (my_rank==q) /*主行所在的进程*/
- {
- temp=v(my_rank*m+r);
- for(i=0;i<r;i++)
- sum[r]=sum[r]+a(r,my_rank*m+i)*v(my_rank*m+i);
- /*计算出的解向量*/
- v(my_rank*m+r)=(b(r)-sum[r]) / a(r,my_rank*m +r);
- if (fabs(v(my_rank*m+r)-temp)<E)
- iteration++;
- /*广播解向量*/
- MPI_Bcast(&v(my_rank*m+r), 1,MPI_FLOAT,my_rank,MPI_COMM_WORLD);
- sum[r]=0.0;
- for(i=0;i<r;i++)
- sum[i]=sum[i]+a(i,my_rank*m+r)*v(my_rank*m+r);
- }
- else /*其他进程接受广播来的x[j]的新值*/
- /*各进程对解向量的其它分量计算有关乘积项并求和*/
- {
- MPI_Bcast(&v(q*m+r),1, MPI_FLOAT,q,MPI_COMM_WORLD);
- for(i=0;i<m;i++)
- sum[i]=sum[i]+a(i,q*m +r)*v(q*m+r);
- }
- }
- /*通过归约操作的求和运算以决定是否进行下一次迭代*/
- /*allreduce过程求出所有进程iteration和total,并广播到所有进程*/
- MPI_Allreduce(&iteration,&total,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);
- loop++;
- if (my_rank==0)
- printf("in the %d times total vaule = %dn",loop,total);
- }
- }
- /*
- * 函数名: put_data
- * 功能:输出数据
- */
- void put_data()
- {
- if (my_rank==0)
- {
- for(i = 0; i < size; i ++)
- printf("x[%d] = %fn",i,v(i));
- printf("n");
- /*输出迭代次数*/
- printf("Iteration num = %dn",loop);
- }
- }
- /*
- * 函数名: main
- * 功能:主程序
- * 输入:argc为命令行参数个数;
- * argv为每个命令行参数组成的字符串数组。
- * 输出:返回0代表程序正常结束;
- */
- int main(int argc, char **argv)
- {
- loop = total = 0;
- /*MPI程序初始化*/
- MPI_Init(&argc,&argv);
- MPI_Comm_size(MPI_COMM_WORLD, &group_size);
- MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
- p=group_size;
- /*0号进程负责读取数据*/
- read_data();
- /*广播数据*/
- broadcast_data();
- /*所有进程并行地对主对角元素右边的数据求和*/
- cal_cross_data();
- /*迭代求解*/
- die_dai_cal();
- /*0号进程将结果输出*/
- put_data();
- MPI_Barrier(MPI_COMM_WORLD);
- MPI_Finalize();
- Environment_Finalize(a,b,v);
- return (0);
- }