高斯-塞德尔迭代并行算法

高斯-塞德尔迭代并行算法在并行计算中，高斯-塞德尔迭代采用与雅可比迭代相同的数据划分。对于高斯-塞德尔迭代，计算ix的新值时，使用11,,inxx的旧值和01,,ixx的新值。计算过程中ix与01,,ixx及11,,inxx的新值会在不同的处理器中产生，因此可以考虑采用时间偏移的方法，使各个处理器对新值计算的开始和结束时间产生一定的偏差。编号为my_rank的处理器一旦计算出(_(_1))ixmyrankmimyrankm的新值，就立即广播给其余处理器，以供各处理器对x的其它分量计算有关ix的乘积项并求和。当它计算完x的所有分量后，它还要接收其它处理器发送的新的x分量，并对这些分量进行求和计算，为计算下一轮的ix作准备。计算开始时，所有处理器并行地对主对角元素右边的数据项进行求和，此时编号为0的处理器（简称为0p）计算出0x,然后广播给其余处理器，其余所有的处理器用0x的新值和其对应项进行求和计算,接着0p计算出12,,,xx当0p完成对1mx的计算和广播后，1p计算出mx，并广播给其余处理器，其余所有的处理器用mx的新值求其对应项的乘积并作求和计算。然后1p计算出12,,,mmxx当1p完成对2*1mx的计算和广播后，2p计算出2*mx，如此重复下去，直至1nx在1pp中被计算出并广播至其余的处理器之后，0p计算出下一轮的新的0x，这样逐次迭代下去，直至收敛为止。具体算法框架描述如下：算法1求解线性方程组的高斯-塞德尔迭代并行算法输入：系数矩阵nnA，常数向量nnb，ε，初始解向量nnx输出：解向量nnxBegin对所有处理器my_rank(my_rank=0,…,p-1)同时执行如下的算法:(1)fori=my_rank*mto(my_rank+1)*m-1do/*所有处理器并行地对主对角元素右边的数据求和*/(1.1)sum[i]=0.0(1.2)forj=i+1ton-1dosum[i]=sum[i]=+a[i,j]*x[j]endforendfor(2)while(totaln)do/*total为新旧值之差小于ε的x的分量个数*/(2.1)iteration=0/*iteration为本处理器中新旧值之差小于ε的x的分量个数*/(2.2)forj=0ton-1do/*依次以第0,1,…,n-1行为主行*/(i)q=j/m(ii)ifmy_rank=qthen/*主行所在的处理器*/tempxj,xjbjsumj/aj,j/*产生xj的新的值*/if(xjtemp││)theniteration=iteration+1endif将x[j]的新值广播到其它所有处理器中/*对其余行计算x[j]所对于的内积项并累加*/sum[j]=0fori=my-rank*mto(my-rank+1)*m-1doif(j≠i)thensumisumiai,j*xjendifendforelse/*其它处理器*/接收广播来的x[j]的新值/*对所存各行计算x[j]所对于的内积项并累加*/forimyrank*mtomyrank1*m1do sumisumiai,j*xjendforendifendfor(2.3)用Allreduce操作求出所有处理器中iteration值的和total并广播到所有处理器中endwhileend若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间。在算法开始阶段，各处理器并行计算主对角线右边元素相应的乘积项并求和,所需时间mn1mm/2，进入迭代计算后，虽然各个处理器所负责的x的分量在每一轮计算中的开始时间是不一样的，但一轮迭代中的计算量都是相等的，因此不妨取0号处理器为对象分析一轮迭代计算的时间，容易得出0号处理器计算时间为mnm；另外它在一轮迭代中做广播操作n次，通信量为1，归约操作1次，通信量为1，所有的通信时间为()log2(1)(1)swswnttptptp，因此高斯-塞德尔迭代的一轮并行计算时间为()log2(1)(1)pswswTmnmnttptptp。源程序如下：1.源程序seidel.c#includestdio.h#includestdlib.h#includempi.h#includemath.h#defineE0.0001#definea(x,y)a[x*size+y]#defineb(x)b[x]#definev(x)v[x]/*A为size*size的系数矩阵*/#defineA(x,y)A[x*size+y]#defineB(x)B[x]#defineV(x)V[x]#defineintsizesizeof(int)#definefloatsizesizeof(float)#definecharsizesizeof(char)intsize,N;intm;float*B;float*A;float*V;intmy_rank;intp;MPI_Statusstatus;FILE*fdA,*fdB,*fdB1;voidEnvironment_Finalize(float*a,float*b,float*v){free(a);free(b);free(v);}intmain(intargc,char**argv){inti,j,my_rank,group_size;intk;float*sum;float*b;float*v;float*a;float*differ;floattemp;intiteration,total,loop;intr,q;loop=0;total=0;MPI_Init(&argc,&argv);MPI_Comm_size(MPI_COMM_WORLD,&group_size);MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);p=group_size;if(my_rank==0){fdA=fopen(dataIn.txt,r);fscanf(fdA,%d%d,&size,&N);if(size!=N-1){printf(theinputiswrong\n);exit(1);}A=(float*)malloc(floatsize*size*size);B=(float*)malloc(floatsize*size);V=(float*)malloc(floatsize*size);for(i=0;isize;i++){for(j=0;jsize;j++)fscanf(fdA,%f,A+i*size+j);fscanf(fdA,%f,B+i);}for(i=0;isize;i++)fscanf(fdA,%f,V+i);fclose(fdA);printf(Inputoffile\dataIn.txt\\n);printf(%d\t%d\n,size,N);for(i=0;isize;i++){for(j=0;jsize;j++)printf(%f\t,A(i,j));printf(%f\n,B(i));}printf(\n);for(i=0;isize;i++)printf(%f\t,V(i));printf(\n\n);printf(\nOutputofresult\n);}/*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(allocatespacefail!);if(my_rank==0){for(i=0;isize;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){for(i=0;im;i++)for(j=0;jsize;j++)a(i,j)=A(i,j);for(i=0;im;i++)b(i)=B(i);for(i=1;ip;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{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);}/*所有处理器并行地对主对角元素右边的数据求和*/for(i=0;im;i++){sum[i]=0.0;for(j=0;jsize;j++)if(j(my_rank*m+i))sum[i]=sum[i]+a(i,j)*v(j);}while(totalsize){iteration=0;total=0;for(j=0;jsize;j++){r=j%m;q=j/m;/*编号为q的处理器负责计算解向量并广播给所有处理器*/if(my_rank==q){temp=v(my_rank*m+r);for(i=0;ir;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;ir;i++)sum[i]=sum[i]+a(i,my_rank*m+r)*v(my_rank*m+r);}else/*各处理器对解向量的其它分量计算有关乘积项并求和*/{MPI_Bcast(&v(q*m+r),1,MPI_FLOAT,q,MPI_COMM_WORLD);for(i=0;im;i++)sum[i]=sum[i]+a(i,q*m+r)*v(q*m+r);}}/*通过归约操作的求和运算以决定是否进行下一次迭代*/MPI_Allreduce(&iteration,&total,1,MPI_FLOAT,MPI_SUM,MPI_COMM_WORLD);loop++;if(my_rank==0)printf(inthe%dtimestotalvaule=%d\n,loop,total);}if(my_rank==0){for(i=0;isize;i++)printf(x[%d]=%f\n,i,v(i));printf(\n);}printf(Iterationnum=%d\n,loop);MPI_Barrier(MPI_COMM_WORLD);MPI_Finalize();Environment_Finalize(a,b,v);return(0);运行实例编译：mpicc–oseidelseidel.cc运行：可以使用命令mpiru