8线性方程组的迭代解法 在阶数较大、系数阵为稀疏阵的情况下,可以采用迭代法求解线性方程组。用迭代法 ( Iterative Method)求解线性方程组的优点是方法简单,便于编制计算机程序,但必须选取合 适的迭代格式及初始向量,以使迭代过程尽快地收敛。迭代法根据迭代格式的不同分成雅可 比 Jacobi)迭代、高斯塞德尔( Gauss-Seidel迭代和松弛 Relaxation)法等几种。在本节中,我 们假设系数矩阵A的主对角线元素an≠0,且按行严格对角占优 Diagonal Donimant),即: la>∑nl 1.1雅可比迭代 111雅可比选代及其串行算法 雅可比迭代的原理是:对于求解n阶线性方程组A=b,将原方程组的每一个方程anx1+ x2+…+ aint=b改写为未知向量x的分量的形式: x=(b-∑可x)/an(1≤i≤n) j=l,j≠ 然后使用第k-1步所计算的变量x1)来计算第k步的x的值: 4)=(b-∑ax-)/an(≤;k≤n) 这里,x为第k次迭代得到的近似解向量x=(x(,x2,…,x)的第i个分量。取 适当初始解向量x0代入上述迭代格式中,可得到x1),再由x)得到x2),依次迭代下去得到 近似解向量序列{x)}。若原方程组的系数矩阵按行严格对角占优,则{x)}收敛于原方程组 的解x。实际计算中,一般认为,当相邻两次的迭代值x+)与x=(1,2,…n)很接近时, x*+)与准确解x中的分量x也很接近。因此,一般用max+x,判断迭代是否收敛 如果取一次乘法和加法运算时间或一次比较运算时间为一个单位时间,则下述雅可比迭代算 法20.1的一轮计算时间为m2+n=O(n2) 算法201单处理器上求解线性方程组雅可比选代算法 输入:系数矩阵Anxn,常数向量bnx1,ε,初始解向量xx1 输出:解向量xn× (I)for i=l to n do end for
1. 8 线性方程组的迭代解法 在阶数较大、系数阵为稀疏阵的情况下，可以采用迭代法求解线性方程组。用迭代法 (Iterative Method)求解线性方程组的优点是方法简单，便于编制计算机程序，但必须选取合 适的迭代格式及初始向量，以使迭代过程尽快地收敛。迭代法根据迭代格式的不同分成雅可 比(Jacobi)迭代、高斯-塞德尔(Gauss-Seidel)迭代和松弛(Relaxation)法等几种。在本节中，我 们假设系数矩阵 A 的主对角线元素 aii 0 ，且按行严格对角占优(Diagonal Donimant)，即： ( 1,2,..., ) 1 a a i n n j j i ii ij = = 1.1 雅可比迭代 1.1.1 雅可比迭代及其串行算法 雅可比迭代的原理是：对于求解 n 阶线性方程组 Ax=b，将原方程组的每一个方程 ai1x1+ ai2x2+…+ ainxn= bi 改写为未知向量 x 的分量的形式： ( ) / (1 ) 1, x b a x aii i n n j j i i = i − ij j = 然后使用第 k-1 步所计算的变量 xi (k -1)来计算第 k 步的 xi (k)的值： ( ) / (1 , ) 1, ( ) ( 1) x b a x aii i k n n j j i k i ij k i j = − = − 这里，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 也很接近。因此，一般用 (k) i (k ) i x -x 1 1 i n max + 判断迭代是否收敛。 如果取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则下述雅可比迭代算 法 20.1 的一轮计算时间为 n 2+n=O(n 2 )。 算法 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=c ()while(diff> a)do (3.1)dif=0 (3.2)for il to n do (iifor=I to n do if(≠D)th end if (.3for il to (i)diff-max diff, Inewx-x? (il)x=newx nd for end while 112雅可比迭代并行算法 A是一个n阶系数矩阵,b是n维向量,在求解线性方程组A=b时,若处理器个数为p, 则可对矩阵A按行划分以实现雅可比迭代的并行计算。设矩阵被划分为p块,每块含有连 续的m行向量,这里m=「m/pl,编号为i的处理器含有A的第m至第(计+1)m-1行数据,同 时向量b中下标为m至(计+1)m-1的元素也被分配至编号为i的处理器(=0,1,…p1),初始 解向量x被广播给所有处理器。 在迭代计算中,各处理器并行计算解向量x的各分量值,编号为i的处理器计算分量 xm至x+1)m1的新值。并求其分量中前后两次值的最大差 localmax,然后通过归约操作的求 最大值运算求得整个n维解向量中前后两次值的最大差max并广播给所有处理器。最后通 过扩展收集操作将各处理器中的解向量按处理器编号连接起来并广播给所有处理器,以进行 下一次迭代计算,直至收敛。具体算法框架描述如下 算法202求解线性方程组的雅可比迭代并行算法 输入:系数矩阵Axn,常数向量bn×1,ε,初始解向量x×1 输出:解向量xnx1 Begin 对所有处理器 my rank( my rank=0,…,p1)同时执行如下的算法 while(max>e)d (1)for=0tom-1do/各个处理器由所存的行计算出解x的相应的分量* (1.1)sm=0.0 (1.2)for j=0 to n-l do if(≠( my_ rank*m+)then sm=sm+a[*x团 nd if 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 1.1.2 雅可比迭代并行算法 A 是一个 n 阶系数矩阵，b 是 n 维向量，在求解线性方程组 Ax=b 时，若处理器个数为 p， 则可对矩阵 A 按行划分以实现雅可比迭代的并行计算。设矩阵被划分为 p 块，每块含有连 续的 m 行向量，这里 m = n/ p ，编号为 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)xl[=(b0-sum)/a[i, my rank*m+i] end for (2)/求出本处理器计算的x的相应的分量的新值与原值的差的最大值 localmax localmax= x1[0]-x[o] (3)for i=l to m-I do if( xl[]-x[d I>localmax)then localmax=|x1(小-xl end if (4)用 Allgather操作将x的所有分量的新值广播到所有处理器中 (5)用 Allreduce操作求出所有处理器中 locaImax值的最大值max并广播到所有处 理器中 end while 若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间,则一轮迭代的计算 时间为m+m;另外,各处理器在迭代中做一次归约操作,通信量为1,一次扩展收集操作, 通信量为m,需要的通信时间为4,(√p-1)+(m+1)n(p-1),因此算法202的一轮并行计算 时间为Tn=4(√P-)+(m+1n(P-1)+m+m。 MPI源程序请参见所附光盘 12高斯-塞德尔迭代 121高斯-塞德尔迭代及其串行算法 高斯-塞德尔迭代的基本思想与雅可比迭代相似。它们的区别在于,在雅可比迭代中 每次迭代时只用到前一次的迭代值,而在高斯塞德尔迭代中,每次迭代时充分利用最新的 迭代值。一旦一个分量的新值被求出,就立即用于后续分量的迭代计算,而不必等到所有分 量的新值被求出以后。设方程组Ax=b的第i个方程为 b,(i=1,2…n) 高斯-塞德尔迭代公式为 (b1-∑4y 取适当的x0作为初始向量,由上述迭代格式可得出近似解向量{x}。若原方程组的系 数矩阵是按行严格对角占优的,则{x}收敛于方程组的解x,若取一次乘法和加法运算时间 或一次比较运算时间为一个单位时间,则下述高斯塞德尔迭代算法203的一轮计算时间为 算法20.3单处理器上求解线性方程组的高斯塞德尔迭代算法 输入:系数矩阵Anxn,常数向量bnx1,ε,初始解向量xx1 输出:解向量xn×1
(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，需要的通信时间为 4t ( p −1) + (m +1)tw( p −1) s ，因此算法 20.2 的一轮并行计算 时间为 Tp = 4t s ( p −1) + (m +1)tw ( p −1) + mn + m。 MPI 源程序请参见所附光盘。 1.2 高斯-塞德尔迭代 1.2.1 高斯-塞德尔迭代及其串行算法 高斯-塞德尔迭代的基本思想与雅可比迭代相似。它们的区别在于，在雅可比迭代中， 每次迭代时只用到前一次的迭代值，而在高斯-塞德尔迭代中，每次迭代时充分利用最新的 迭代值。一旦一个分量的新值被求出，就立即用于后续分量的迭代计算，而不必等到所有分 量的新值被求出以后。设方程组 Ax=b 的第 i 个方程为： = n j 1 ij a j x = i b (i =1,2, ,n) 高斯-塞德尔迭代公式为： ( ) 1 1 ( ) 1 1 ( 1) ( 1) = − − = + − = + + n j i k ij j i j k i ij j ii k i b a x a x a x (i =1,2, , n) 取适当的 x (0)作为初始向量，由上述迭代格式可得出近似解向量{x (k)}。若原方程组的系 数矩阵是按行严格对角占优的，则{x (k)}收敛于方程组的解 x，若取一次乘法和加法运算时间 或一次比较运算时间为一个单位时间，则下述高斯-塞德尔迭代算法 20.3 的一轮计算时间为 n 2+n=O(n 2 )。 算法 20.3 单处理器上求解线性方程组的高斯-塞德尔迭代算法 输入：系数矩阵 An×n，常数向量 b n×1，ε，初始解向量 xn×1 输出：解向量 xn×1 Begin
d fo (2)p=a 3) while(p≥E)do (1)t=x (i)s= (iiiforj=I to n do if(≠)then end for (v)if(x-t>)then p=x-tend if end for end while 122高斯塞德尔迭代并行算法 在并行计算中,高斯-塞德尔迭代采用与雅可比迭代相同的数据划分。对于高斯塞德尔 迭代,计算x的新值时,使用x+,…,x-1的旧值和xn…x的新值。计算过程中x与x0,…x:1 及x+1,…,xn1的新值会在不同的处理器中产生,因此可以考虑采用时间偏移的方法,使各 个处理器对新值计算的开始和结束时间产生一定的偏差。编号为 my rank的处理器一旦计 算出x( my rank×m≤i<( my rank+1)×m)的新值,就立即广播给其余处理器,以供各处理 器对x的其它分量计算有关x的乘积项并求和。当它计算完x的所有分量后,它还要接收其 它处理器发送的新的x分量,并对这些分量进行求和计算,为计算下一轮的x作准备。计算 开始时,所有处理器并行地对主对角元素右边的数据项进行求和,此时编号为0的处理器(简 称为Po)计算出x然后广播给其余处理器,其余所有的处理器用x0的新值和其对应项进行 求和计算,接着P计算出x,x2…当P完成对xm-1的计算和广播后,P1计算出xm,并广播给 其余处理器,其余所有的处理器用xm的新值求其对应项的乘积并作求和计算。然后P1计算 出xm+1,xm+2,…,当P1完成对xm1的计算和广播后,P2计算出x·m…,如此重复下去,直 至xn在P1中被计算出并广播至其余的处理器之后,P0计算出下一轮的新的x0,这样逐次 迭代下去,直至收敛为止。具体算法框架描述如下 算法20.4求解线性方程组的高斯塞德尔迭代并行算法 输入:系数矩阵Axn,常数向量bn×1,ε,初始解向量x×1 输出:解向量x×1 对所有处理器 my rank( my rank=0,…,p-1)同时执行如下的算法 (I)for i=my-rank m to(my-rank+1)m-I do /*所有处理器并行地对主对角元素右边的数据求和* (1.1)stm{d=0.0 (1.2)forj=i+l to n-l d sumi=sum(i+aixi
(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 1.2.2 高斯-塞德尔迭代并行算法 在并行计算中，高斯-塞德尔迭代采用与雅可比迭代相同的数据划分。对于高斯-塞德尔 迭代，计算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]
d for end for (2) while(toal<n)do/ *total为新旧值之差小于ε的x的分量个数 / iteration为本处理器中新旧值之差小于e的x的分量个数* (2.2)forj=0ton-ldo/*依次以第0,1,…,n1行为主行* (i)i/m (ini( my rank=q)then倖*主行所在的处理器* lemp=x,xU=(b小- suml/al/产生x的新的值* if(xUJ-temp <e)then iteration=iteration+1 end if 将x的新值广播到其它所有处理器中 /*对其余行计算x所对于的内积项并累加* suml=0 for i=my-rank* m to(my-rank+1)*m-I do if(≠ i)then sm[=m{小+a[4小*xU end fo else/*其它处理器* 接收广播来的x的新值 /*对所存各行计算x所对于的内积项并累加* for i=my-rank* m to(my-rank+1)"m-I do sumli=sumi+a[ixlI end for end for (2.3)用 Allreduce操作求出所有处理器中 iteration值的和toal并广播到所有处 理器中 end while 若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间。在算法开始阶段, 各处理器并行计算主对角线右边元素相应的乘积项并求和所需时间m-(1+m)m/2,进入迭代 计算后,虽然各个处理器所负责的x的分量在每一轮计算中的开始时间是不一样的,但一轮 迭代中的计算量都是相等的,因此不妨取0号处理器为对象分析一轮迭代计算的时间,容易 得出0号处理器计算时间为mn+m;另外它在一轮迭代中做广播操作n次,通信量为1,归 约操作1次,通信量为1,所有的通信时间为m1+ln)gp+2(P-1)+tp(p-1),因此高 斯塞德尔迭代的一轮并行计算时间为Tp=m+m+m(+1m)bogp+2,(√P-1+1n(P-1)。 MPI源程序请参见章末附录
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，所有的通信时间为 n(t + t )log p + 2t ( p −1) + tw( p −1) s w s ，因此高 斯-塞德尔迭代的一轮并行计算时间为 T = mn + m + n(t + t )log p + 2t ( p −1) + tw(p −1) p s w s 。 MPI 源程序请参见章末附录