8 线性方程组的迭代解法
在阶数较大、系数阵为稀疏阵的情况下，可以采用迭代法求解线性方程组。用迭代法(Iterative Method)求解线性方程组的优点是方法简单，便于编制计算机程序，但必须选取合适的迭代格式及初始向量，以使迭代过程尽快地收敛。迭代法根据迭代格式的不同分成雅可比(Jacobi)迭代、高斯-塞德尔(Gauss-Seidel)迭代和松弛(Relaxation)法等几种。在本节中，我们假设系数矩阵A的主对角线元素，且按行严格对角占优(Diagonal Donimant)，即：
雅可比迭代
雅可比迭代及其串行算法
雅可比迭代的原理是：对于求解n阶线性方程组Ax=b，将原方程组的每一个方程ai1x1+ ai2x2+…+ ainxn= bi改写为未知向量x的分量的形式：
然后使用第k-1步所计算的变量xi(k -1)来计算第k步的xi(k)的值：
这里，xi(k)为第k次迭代得到的近似解向量x(k)= (x1 (k), x2(k) , …, xn(k) )T的第i个分量。取适当初始解向量x(0)代入上述迭代格式中，可得到x(1)，再由x(1)得到x(2)，依次迭代下去得到近似解向量序列{x(k)}。若原方程组的系数矩阵按行严格对角占优，则{x(k)}收敛于原方程组的解x。实际计算中，一般认为，当相邻两次的迭代值xi(k +1)与xi(k) i=(1,2, …,n)很接近时，xi(k +1)与准确解x中的分量xi也很接近。因此，一般用判断迭代是否收敛。如果取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则下述雅可比迭代算法20.1的一轮计算时间为n2+n=O(n2)。
算法20.1 单处理器上求解线性方程组雅可比迭代算法
输入：系数矩阵An×n，常数向量b n×1，ε，初始解向量xn×1
输出：解向量xn×1
Begin
(1)for i=1 to n do
xi=bi/aii
end for
(2)diff=ε
(3)while (diff ≥ ε) do
(3.1)diff=0
(3.2)for i=1 to n do
(i)newxi= bi
(ii)for j= 1 to n do
if (j ≠ i) then
newxi= newxi- aij xj
end if
end for
(iii)newxi= newxi/ aii
end for
(3.3)for i=1 to n do
(i)diff=max {diff, |newxi- xi|}
(ii)xi=newxi
end for
end while
End
雅可比迭代并行算法
A是一个n阶系数矩阵，b是n维向量，在求解线性方程组Ax=b时，若处理器个数为p，则可对矩阵A按行划分以实现雅可比迭代的并行计算。设矩阵被划分为p块，每块含有连续的m行向量，这里，编号为i的处理器含有A的第im至第(i+1)m-1行数据，同时向量b中下标为im至(i+1)m-1的元素也被分配至编号为i的处理器(i=0,1, …,p-1)，初始解向量x被广播给所有处理器。
在迭代计算中，各处理器并行计算解向量x的各分量值，编号为i的处理器计算分量xim至x(i+1)m-1的新值。并求其分量中前后两次值的最大差localmax，然后通过归约操作的求最大值运算求得整个n维解向量中前后两次值的最大差max并广播给所有处理器。最后通过扩展收集操作将各处理器中的解向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下：
算法20.2 求解线性方程组的雅可比迭代并行算法
输入：系数矩阵An×n，常数向量b n×1，ε，初始解向量xn×1
输出：解向量xn×1
Begin
对所有处理器my_rank(my_rank=0,…, p-1)同时执行如下的算法:
while (max>ε) do
(1)for i=0 to m-1 do /*各个处理器由所存的行计算出解x的相应的分量*/
(1.1)sum=0.0
(1.2)for j=0 to n-1 do
if (j ≠ (my_rank*m+i)) then
sum=sum+a[i,j]*x[j]
end if
end for
(1.3)x1[i]=(b[i] - sum)/a[i,my_rank*m+i]
end for
(2)/*求出本处理器计算的x的相应的分量的新值与原值的差的最大值localmax */
localmax=│x1[0]-x[0]│
(3)for i=1 to m-1 do
if (│x1[i]-x[i] │>localmax) then
localmax =│x1[i]-x[i] │
end if
end for
(4)用Allgather操作将x的所有分量的新值广播到所有处理器中
(5)用Allreduce操作求出所有处理器中localmax值的最大值max并广播到所有处理器中
end while
End
若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则一轮迭代的计算时间为mn+m；另外，各处理器在迭代中做一次归约操作，通信量为1，一次扩展收集操作，通信量为m，需要的通信时间为，因此算法20.2的一轮并行计算时间为。
MPI源程序请参见所附光盘。
高斯-塞德尔迭代
高斯-塞德尔迭代及其串行算法
高斯-塞德尔迭代的基本思想与雅可比迭代相似。它们的区别在于，在雅可比迭代中，每次迭代时只用到前一次的迭代值，而在高斯-塞德尔迭代中，每次迭代时充分利用最新的迭代值。一旦一个分量的新值被求出，就立即用于后续分量的迭代计算，而不必等到所有分量的新值被求出以后。设方程组Ax=b的第i个方程为：
=
高斯-塞德尔迭代公式为：
取适当的x(0)作为初始向量，由上述迭代格式可得出近似解向量{x(k)}。若原方程组的系数矩阵是按行严格对角占优的，则{x(k)}收敛于方程组的解x，若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则下述高斯-塞德尔迭代算法20.3的一轮计算时间为n2+n=O(n2)。
算法20.3 单处理器上求解线性方程组的高斯-塞德尔迭代算法
输入：系数矩阵An×n，常数向量b n×1，ε，初始解向量xn×1
输出：解向量xn×1
Begin
(1)for i=1 to n do
xi=0
end for
(2)p=ε+1
(3)while (p ≥ ε) do
for i=1 to n do
(i) t = xi
(ii) s=0
(iii)for j= 1 to n do
if (j ≠ i) then
s= s+ aij xj
end if
end for
(iv) xi=(bi-s)/ aii
(v) if (│xi-t│>p) then p=│xi-t│end if
end for
end while
End
高斯-塞德尔迭代并行算法
在并行计算中，高斯-塞德尔迭代采用与雅可比迭代相同的数据划分。对于高斯-塞德尔迭代，计算xi的新值时，使用xi+1, …,xn-1的旧值和x0, …,xi-1的新值。计算过程中xi与x0, …,xi-1及xi+1, …,xn-1的新值会在不同的处理器中产生，因此可以考虑采用时间偏移的方法，使各个处理器对新值计算的开始和结束时间产生一定的偏差。编号为my_rank的处理器一旦计算出xi(my_rank×m ≤ i < (my_rank+1)×m)的新值，就立即广播给其余处理器，以供各处理器对x的其它分量计算有关xi的乘积项并求和。当它计算完x的所有分量后，它还要接收其它处理器发送的新的x分量，并对这些分量进行求和计算，为计算下一轮的xi作准备。计算开始时，所有处理器并行地对主对角元素右边的数据项进行求和，此时编号为0的处理器（简称为P0）计算出x0,然后广播给其余处理器，其余所有的处理器用x0的新值和其对应项进行求和计算,接着P0计算出x1,x2, …,当P0完成对xm-1的计算和广播后，P1计算出xm，并广播给其余处理器，其余所有的处理器用xm的新值求其对应项的乘积并作求和计算。然后P1计算出xm+1,xm+2, …,当P1完成对x2*m-1 的计算和广播后，P2计算出x2*m …，如此重复下去，直至xn-1在Pp-1中被计算出并广播至其余的处理器之后，P0计算出下一轮的新的x0，这样逐次迭代下去，直至收敛为止。具体算法框架描述如下：
算法20.4 求解线性方程组的高斯-塞德尔迭代并行算法
输入：系数矩阵An×n，常数向量b n×1，ε，初始解向量xn×1
输出：解向量xn×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 /*主行所在的处理器*/
temp= x[j], x[j]= (b[j]- sum[j])/a[j,j] /* 产生x(j)的新的值*/
if (│x[j]-temp│<ε) 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
sum[i]=sum[i]+a[i,j]*x[j]
end if
end for
else /*其它处理器*/
接收广播来的x[j]的新值
/*对所存各行计算x[j]所对于的内积项并累加*/
for i=my-rank* m to (my-rank+1)* m-1 do
sum[i]=sum[i]+a[i,j]*x[j]
end for
end if
end for
(2.3)用Allreduce操作求出所有处理器中iteration值的和total并广播到所有处理器中
end while
End
若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间。在算法开始阶段，各处理器并行计算主对角线右边元素相应的乘积项并求和,所需时间mn-(1+m)m/2，进入迭代计算后，虽然各个处理器所负责的x的分量在每一轮计算中的开始时间是不一样的，但一轮迭代中的计算量都是相等的，因此不妨取0号处理器为对象分析一轮迭代计算的时间，容易得出0号处理器计算时间为mn+m；另外它在一轮迭代中做广播操作n次，通信量为1，归约操作1次，通信量为1，所有的通信时间为，因此高斯-塞德尔迭代的一轮并行计算时间为。
MPI源程序请参见章末附录。
松弛法
松弛法及其串行算法
为了加快雅可比迭代与高斯-塞德尔迭代的收敛速度，可采用松弛法。以高斯-塞德尔迭代为例，首先，由高斯-塞德尔迭代格式求得第k+1次迭代值xi(k +1)，即：
然后计算第k+1次迭代值与第次迭代值之差；即：
i(k +1)-xi(k)
最后在第k次迭代值的基础上，加上这个差的一个倍数作为实际的第k+1次迭代值，即：
xi(k +1)= xi(k)+w(i(k +1)-xi(k))
其中w为一个常数。综合以上过程，可以得到如下迭代格式：
其中， w称为松弛因子，为保证收敛，要求松弛因子w满足：0<w<2。当w>1时，称之为超松弛法，此时，i(k +1)的比重被加大；当w=1时，它就成了一般的高斯-塞德尔迭代。实际应用时，可根据系数矩阵的性质及对其反复计算的经验来决定适合的松弛因子w。若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则下述算法20.5一轮计算时间计算为n2+n=O(n2)。
算法20.5 单处理器上松弛法求解线性方程组的算法
输入：系数矩阵An×n，常数向量b n×1，ε，初始解向量xn×1
输出：解向量xn×1
Begin
(1) for i=1 to n do
xi=0
end for
(2) p=ε+1
(3) while ( p ≥ ε) do
(3.1) for i=1 to n do
(i) s=0
(ii) for j= 1 to n do
if (j ≠ i) then
s= s+ aij xj
end if
end for
(iii) yi= (1-w) xi +w(bi-s)/aii
(iv) if (│yi - xi│>p) then p=│yi - xi│end if
(v) xi =yi
end for
end while
End
松弛法并行算法
松弛法并行算法和高斯-塞德尔迭代并行算法基本相同，只是在计算分量的时候，将对x分量的新值的计算改为“x[j]=(1-w)*temp+w* (b[j]- sum[j])/a[j,j]”，其中temp为x[j]的原有值。具体并行算法描述如下:
算法20.6 求解线性方程组的松弛迭代并行算法
输入：系数矩阵An×n，常数向量b n×1，ε，初始解向量xn×1
输出：解向量xn×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 /*主行所在的处理器*/
temp= x[j]
x[j]=(1-w)*temp+w* (b[j]- sum[j])/a[j,j] /* 产生x[j]的新的值*/
if (│x[j]-temp│<ε) 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
sum[i]=sum[i]+a[i,j]*x[j]
end if
end for
(iii) else /*其它处理器*/
接收广播来的x[j]的新值
/*对所存各行计算x[j]所对于的内积项并累加*/
for i=my-rank* m to (my-rank+1)* m-1 do
sum[I]=sum[i]+a[i,j]*x[j]
end for
end if
end for
(2.3) 用Allreduce操作求出所有处理器中iteration值的和total并广播到所有处理器中
end while
End
与并行高斯-塞德尔迭代相似,并行松弛迭代法的一轮并行计算时间为: 。
MPI源程序请参见所附光盘。
小结
本章主要讨论线性方程组的迭代解法，这种方法是一种逐步求精的近似求解过程，其优点是简单，易于计算机编程，但它存在着迭代是否收敛以及收敛速度快慢的问题。一般迭代过程由预先给定的精度要求来控制，但由于方程组的准确解一般是不知道的，因此判断某次迭代是否满足精度要求也是比较困难的，需要根据具体情况而定。文献[1]给出了稀疏线性方程组迭代解法的详尽描述，还包含了多重网格(Multigrid)法、共轭梯度(Conjugate Gradient)法，[2]综述了稀疏线性方程组的并行求解算法，[3]综述了在向量机和并行机上偏微分方程的求解方法，[4]讨论了超立方多处理机上的多重网格算法，[5]讨论了并行共轭梯度算法，[6]深入而全面地论述了SIMD和MIMD模型上的数值代数、离散变换和卷积、微分方程、计算数论和最优化计算的并行算法，对并行排序算法也作了介绍。此外，在[7]中第九章的参考文献注释中还列举了大量有关参考文献，进一步深入研究的读者可在这些文献中获得更多的资料。
参考文献
陈国良 编著．并行计算——结构·算法·编程．高等教育出版社,1999.10
Heath M T, Ng E and Peyton B W. Parallel Algorithm for Space Linear Systems. SIAM Review,1991,33:420-460
Ortega J M, Voigt R G. Solution of Partial Differential Equations on Vector and Parallel Computers. SIAM Review,1985,27(2):149-240
Chan T F,Saad Y. Multigrid Algorithms on the Hypercube Multiprocessor. IEEE-TC,1986,C-35(11):969-977
Chronopoulos A T, Gear C W. On the Efficient Implementation of Pre-condition S-step Conjugate Gradient Methods on Multiprocessors with Memory Hierarchy. Parallel Computing,1989,11:37-53
李晓梅，蒋增荣 等编著．并行算法(第五章) ．湖南科技出版社, 1992
Quinn M J. Parallel Computing-Theory and Practic(second edition)McGraw-Hill, Inc., 1994
附录 高斯-塞德尔迭代并行算法的MPI源程序
源程序 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。
