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

高斯-塞德尔迭代并行算法
在并行计算中，高斯-塞德尔迭代采用与雅可比迭代相同的数据划分。对于高斯-塞德尔迭代，计算的新值时，使用的旧值和的新值。计算过程中与及的新值会在不同的处理器中产生，因此可以考虑采用时间偏移的方法，使各个处理器对新值计算的开始和结束时间产生一定的偏差。编号为 my_rank 的处理器一旦计算出的新值，就立即广播给其余处理器，以供各处理器对x 的其它分量计算有关的乘积项并求和。当它计算完 x的所有分量后，它还要接收其它处理器发送的新的分量，并对这些分量进行求和计算，为计算下一轮的作准备。计算开始时，所有处理器并行地对主对角元素右边的数据项进行求和，此时编号为 0 的处理器 （简称为）计算出,然后广播给其余处理器，其余所有的处理器用的新值和其对应项进行求和计算,接着计算出当完成对的计算和广播后，计算出，并广播给其余处理器，其余所有的处理器用的新值求其对应项的乘积并作求和计算。然后计算出当完成对的计算和广播后，计算出，如此重复下去，直至在中被计算出并广播至其余的处理器之后，计算出下一轮的新的，这样逐次迭代下去，直至收敛为止。具体算法框架描述如下：
算法 1 求解线性方程组的高斯-塞德尔迭代并行算法
输入：系数矩阵，常数向量，ε，初始解向量
输出：解向量
Begin
对所有处理器 my_rank(my_rank=0,…, p-1)同时执行如下的算法:
(1) for i=my_rank*m to (my_rank+1)*m-1 do
/*所有处理器并行地对主对角元素右边的数据求和*/
(1.1) sum[i]=0.0
(1.2)for j=i+1 to n-1 do
sum[i]= sum[i]=+ a[i,j]*x[j]
end for
end for
(2) while (total<n) do /*total为新旧值之差小于 ε的x的分量个数*/
(2.1) iteration=0
/* iteration 为本处理器中新旧值之差小于 ε的x 的分量个数*/
(2.2)for j=0 to n-1 do /*依次以第 0,1, …, n-1行为主行*/
(i) q=j/m
(ii)if my_rank=q then /*主行所在的处理器*/
/* 产生的新的值*/
if () then iteration= iteration +1 end if
将x[j]的新值广播到其它所有处理器中
/*对其余行计算x[j]所对于的内积项并累加*/
sum[j]=0
for i=my-rank* m to (my-rank+1)*m-1 do
if (j ≠ i) then

end if
end for
else /*其它处理器*/
接收广播来的x[j]的新值
/*对所存各行计算x[j]所对于的内积项并累加*/
for do

end for
end if
end for
(2.3)用Allreduce操作求出所有处理器中iteration值的和total并广播到所有处
理器中
end while
end
若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间。在算法开始阶段，各处理器并行计算主对角线右边元素相应的乘积项并求和,所需时间， 进入迭代计算后，虽然各个处理器所负责的的分量在每一轮计算中的开始时间是不一样的，但一轮迭代中的计算量都是相等的，因此不妨取0 号处理器为对象分析一轮迭代计算的时间，容易得出0号处理器计算时间为；另外它在一轮迭代中做广播操作 n 次，通信量为 1，归约操作1 次，通信量为1，所有的通信时间为 ，因此高斯-塞德尔迭代的一轮并行计算时间为。
源程序如下：
1. 源程序 seidel.c
#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;

void Environment_Finalize(float *a,float *b,float *v)
{
free(a);
free(b);
free(v);
}

int main(int argc, char **argv)
{
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;
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("the input is wrong\n");
exit(1);
}
A=(float *)malloc(floatsize*size*size);
B=(float *)malloc(floatsize*size);
V=(float *)malloc(floatsize*size);

for(i = 0; i < size; i++)
{
for(j = 0; j < size; j++)
fscanf(fdA,"%f", A+i*size+j);
fscanf(fdA,"%f", B+i);
}
for(i = 0; i < size; i ++)
fscanf(fdA,"%f", V+i);
fclose(fdA);
printf("Input of file \"dataIn.txt\"\n");
printf("%d\t%d\n", size, N);
for(i = 0; i < size; i ++)
{
for(j = 0; j < size; j ++)
printf("%f\t",A(i,j));
printf("%f\n",B(i));
}
printf("\n");
for(i = 0; i < size; i ++)
printf("%f\t", V(i));
printf("\n\n");
printf("\nOutput of result\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("allocate space fail!");
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)
{
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);
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
{
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;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);
}
while (total<size)
{
iteration=0;
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
/*各处理器对解向量的其它分量计算有关乘积项并求和*/
{
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);
}
}
/*通过归约操作的求和运算以决定是否进行下一次迭代*/
MPI_Allreduce(&iteration,&total,1, MPI_FLOAT,MPI_SUM, MPI_COMM_WORLD);
loop++;
if (my_rank==0)
printf("in the %d times total vaule = %d\n",loop,total);
}
if (my_rank==0)
{
for(i = 0; i < size; i ++)
printf("x[%d] = %f\n",i,v(i));
printf("\n");
}
printf("Iteration num = %d\n",loop);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
Environment_Finalize(a,b,v);
return (0);
运行实例
编译：mpicc –o seidel seidel.cc
运行：可以使用命令mpirun–np ProcessSize seidel 来运行该程序，其中 ProcessSize是所使用的处理器个数,这里取为4。本实例中使用了mpirun –np 4 seidel
运行结果：
Input of file "dataIn.txt"
3 4
9.000000 -1.000000 -1.000000 7.000000
-1.000000 8.000000 0.000000 7.000000
-1.000000 0.000000 9.000000 8.000000
0.000000 0.000000 1.000000
Output of result
in the 1 times total vaule = 0
in the 2 times total vaule = 0
in the 3 times total vaule = 0
in the 4 times total vaule = 3
x[0] = 0.999998
x[1] = 1.000000
x[2] = 1.000000
Iteration num = 4
说明：该运行实例中，A为3×4的矩阵，B是长度为3的向量，它们的值存放于文档“dataIn.txt”
中，其中前3行为矩阵A，最后一行为向量 B，最后输出线性方程组 AX=B的解向量X。