8线性方程组的迭代解法 在阶数较大、系数阵为稀疏阵的情况下,可以采用迭代法求解线性方程组。用迭代法 ( Iterative Method)求解线性方程组的优点是方法简单,便于编制计算机程序,但必须选取合 适的迭代格式及初始向量,以使迭代过程尽快地收敛。迭代法根据迭代格式的不同分成雅可 比 Jacobi)迭代、高斯塞德尔( Gauss-Seidel迭代和松弛 Relaxation)法等几种。在本节中,我 们假设系数矩阵A的主对角线元素an≠0,且按行严格对角占优 Diagonal Donimant),即: la>∑nl(=12…,n) 1.1雅可比迭代 111雅可比选代及其串行算法 雅可比迭代的原理是:对于求解n阶线性方程组Ac=b,将原方程组的每一个方程anx1+ x2+…+ aint=b改写为未知向量x的分量的形式: x=(b-∑ajxj)/ai(≤i≤n) 然后使用第k-1步所计算的变量x1)来计算第k步的x的值: =(b 人9(k-)a(sk≤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,ε,初始解向量xx 输出:解向量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 (inews= b (iifor=I to n do if(≠D)then newx- newx而x d if (.3for il to (i)diff==max ( diff, Inewx-xilP (ii x=. 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求解线性方程组的雅可比迭代并行算法 输入:系数矩阵Anxn,常数向量bnx1,E,初始解向量xn×1 输出:解向量xnx1 Begin 对所有处理器 my rank( my rank=0,…,p1)同时执行如下的算法 hile(max>c)d (1)for=0tom-1do/各个处理器由所存的行计算出解x的相应的分量* (1.1)sm=0.0 (1.2)for j=0 to n-l do if(+(my rank*m+i))then sum=suma[, xll 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 end for (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源程序请参见所附光盘。 1.2高斯-塞德尔迭代 121高斯塞德尔迭代及其串行算法 高斯-塞德尔迭代的基本思想与雅可比迭代相似。它们的区别在于,在雅可比迭代中 每次迭代时只用到前一次的迭代值,而在高斯塞德尔迭代中,每次迭代时充分利用最新的 迭代值。一旦一个分量的新值被求出,就立即用于后续分量的迭代计算,而不必等到所有分 量的新值被求出以后。设方程组Ax=b的第i个方程为 ∑anx b,(i=1,2…n) 高斯-塞德尔迭代公式为 (b1-∑ax 取适当的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 (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+1,…,x-1的旧值和xn…x的新值。计算过程中x与x0,…x1 及x+1,…,xn1的新值会在不同的处理器中产生,因此可以考虑采用时间偏移的方法,使各 个处理器对新值计算的开始和结束时间产生一定的偏差。编号为 my rank的处理器一旦计 算出x( my rank×m≤i<( my rank+1)×m)的新值,就立即广播给其余处理器,以供各处理 器对x的其它分量计算有关x的乘积项并求和。当它计算完x的所有分量后,它还要接收其 它处理器发送的新的x分量,并对这些分量进行求和计算,为计算下一轮的x作准备。计算 开始时,所有处理器并行地对主对角元素右边的数据项进行求和,此时编号为0的处理器(简 称为Po)计算出x然后广播给其余处理器,其余所有的处理器用x0的新值和其对应项进行 求和计算,接着Po计算出x,x2…当P完成对xm-1的计算和广播后,P1计算出xm,并广播给 其余处理器,其余所有的处理器用xm的新值求其对应项的乘积并作求和计算。然后P1计算 出xm+1,xm+2,…,当P1完成对xm1的计算和广播后,P2计算出x2·m…,如此重复下去,直 至xn在P1中被计算出并广播至其余的处理器之后,P0计算出下一轮的新的x0,这样逐次 迭代下去,直至收敛为止。具体算法框架描述如下 算法20.4求解线性方程组的高斯塞德尔迭代并行算法 输入:系数矩阵Axn,常数向量bn×1,ε,初始解向量x×1 输出:解向量xn× 对所有处理器 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)for=i+l to n-l do sumi=sum(i+axl]
(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]
end for end for (2) while(toal<n)do/ *total为新旧值之差小于ε的x的分量个数 / iteration为本处理器中新旧值之差小于e的x的分量个数* (2.2)for戶=0ton-ldo/*依次以第O,1, 1,…,n-1行为主行* (i)i/m (ini( my rank=q)then倖*主行所在的处理器* lemp=x,xU=(b小-smUy/*产生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(≠D)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+t)gp+2(P-1)+tp(p-1),因此高 斯塞德尔迭代的一轮并行计算时间为Tp=m+m+m(1+1n)gP+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 源程序请参见章末附录
13松弛法 131松弛法及其串行算法 为了加快雅可比迭代与高斯塞德尔迭代的收敛速度,可采用松弛法。以高斯-塞德尔迭 代为例,首先,由高斯塞德尔迭代格式求得第k+1次迭代值x+,即 (b2 a r(+l) 然后计算第k+1次迭代值与第k次迭代值之差;即: 最后在第k次迭代值的基础上,加上这个差的一个倍数作为实际的第k1次迭代值,即: v(x{+1-x) 其中w为一个常数。综合以上过程,可以得到如下迭代格式: (1-)x1+w(b-∑ax 其中,w称为松弛因子,为保证收敛,要求松弛因子v满足:01时,称 之为超松弛法,此时,体+)的比重被加大:当=1时,它就成了一般的高斯塞德尔迭代。 实际应用时,可根据系数矩阵的性质及对其反复计算的经验来决定适合的松弛因子w。若取 次乘法和加法运算时间或一次比较运算时间为一个单位时间,则下述算法20.5一轮计算 时间计算为n2+m=O(m)。 算法20.5单处理器上松弛法求解线性方程组的算法 输入:系数矩阵Axn,常数向量bnx1,E,初始解向量x×1 输出:解向量xn×1 (1)for i=l to n de end for (2)p=+1 (3) while(p≥ε)do (3. 1)for i=l to n do (ii)for j= if(≠then S St aij xj (i)y=(1-)x+(b-s)ai (iv)if (ly-xi >p)then p=y-xi end if (v)xi=i end for end while
1.3 松弛法 1.3.1 松弛法及其串行算法 为了加快雅可比迭代与高斯-塞德尔迭代的收敛速度，可采用松弛法。以高斯-塞德尔迭 代为例，首先，由高斯-塞德尔迭代格式求得第 k+1 次迭代值 xi (k +1)，即： ( ) 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 然后计算第 k+1 次迭代值与第 k 次迭代值之差；即： x i (k +1) -xi (k) 最后在第 k 次迭代值的基础上，加上这个差的一个倍数作为实际的第 k+1 次迭代值，即： xi (k +1)= xi (k)+w( x i (k +1) -xi (k) ) 其中 w 为一个常数。综合以上过程，可以得到如下迭代格式： ii n j i k ij j i j k i ij j k i k i a x w x w b a x a x 1 (1 ) ( ) 1 ( ) 1 1 ( 1) ( ) ( 1) = − + − − = + − = + + (i =1,2, , n) 其中， w 称为松弛因子，为保证收敛，要求松弛因子 w 满足：01 时，称 之为超松弛法，此时， x i (k +1)的比重被加大；当 w=1 时，它就成了一般的高斯-塞德尔迭代。 实际应用时，可根据系数矩阵的性质及对其反复计算的经验来决定适合的松弛因子 w。若取 一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则下述算法 20.5 一轮计算 时间计算为 n 2+n=O(n 2 )。 算法 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 132松弛法并行算法 松弛法并行算法和高斯-塞德尔迭代并行算法基本相同,只是在计算分量的时候,将对x 分量的新值的计算改为“x=(1-)*lemp+w*(b小-m)”,其中temp为功的原有值。 具体并行算法描述如下 算法20.6求解线性方程组的松弛迭代并行算法 输入:系数矩阵A2xn,常数向量bnx1,E,初始解向量xnx1 输出:解向量xn×1 对所有处理器 my rank( my rank=0,…,p-l)同时执行如下的算法 /*所有处理器并行地对主对角元素右边的数据求和 (1) for i=my-rank*m to(my-rank+ 1)*m-l do (1.1)sm{=0.0 (1.2)for=计+lton-ldo sm[=sm{小+a/*xU end for end for (2) while(toak<n)do/ total为新旧值之差小于c的x的分量个数* (2.1) iteration=0/ IteratIon为本处理器新旧值之差小于E的x的分量个数* (2.2)for户=0ton-1do/依次以第0,1,…,n-1行为主行* (i)if( my rank=q)then/*主行所在的处理器* x=(1-)*temp+n*b小]m)a*产生x的新的值* if(xlJ-temp<)then iteration= iteration +l end if 将x)的新值广播到其它所有处理器中 *对其余行计算x所对于的内积项并累加* for i=my-rank* m to( my-rank+1)*m-1 do if(≠then end if (i)else/其它处理器* 接收广播来的x]的新值 *对所存各行计算x所对于的内积项并累加* for i=my-rank* m to( my-rank+1)*m-1 do sm{=m{]+a[* nd fo end if (2.3)用 Allreduce操作求出所有处理器中 iteration值的和 total并广播到所有
End 1.3.2 松弛法并行算法 松弛法并行算法和高斯-塞德尔迭代并行算法基本相同，只是在计算分量的时候，将对 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 与并行高斯-塞德尔迭代相似,并行松弛迭代法的一轮并行计算时间为 T,=mn+m+m( s +iw)log p+2, ( p-1)+,(p-I) MPI源程序请参见所附光盘 14小结 本章主要讨论线性方程组的迭代解法,这种方法是一种逐步求精的近似求解过程,其优 点是简单,易于计算机编程,但它存在着迭代是否收敛以及收敛速度快慢的问题。一般迭代 过程由预先给定的精度要求来控制,但由于方程组的准确解一般是不知道的,因此判断某次 迭代是否满足精度要求也是比较困难的,需要根据具体情况而定。文献给出了稀疏线性 方程组迭代解法的详尽描述,还包含了多重网格( Multigrid法、共轭梯度( Conjugate Gradient 法,[2]综述了稀疏线性方程组的并行求解算法,[3]综述了在向量机和并行机上偏微分方程 的求解方法,[4]讨论了超立方多处理机上的多重网格算法,[5]讨论了并行共轭梯度算法 6]深入而全面地论述了SMD和MMD模型上的数值代数、离散变换和卷积、微分方程、 计算数论和最优化计算的并行算法,对并行排序算法也作了介绍。此外,在[刁中第九章的 参考文献注释中还列举了大量有关参考文献,进一步深入研究的读者可在这些文献中获得更 多的资料 参考文献 [1]陈国良编著.并行计算——结构·算法·编程.高等教育出版社,199910 12. Heath M T, Ng e and Peyton B W. Parallel Algorithm for Space Linear Systems. SIAM Review,1991,33:420-460 [3]. Ortega J M, VoigtRG. Solution of Partial Differential Equations on Vector and Parallel Computers. SIAM Review, 1985, 27(2): 149-240 [4]. Chan T F, Saad Y. Multigrid Algorithms on the Hypercube Multiprocessor IEEE-TC,1986C-35(11):969977 15. 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 6]李晓梅,蒋増荣等编著.并行算法(第五章).湖南科技出版社,1992 [7]. Quinn MJ. Parallel Computing-Theory and Practic(second edition )McGraw-Hill, Inc, 1994 附录高斯-塞德尔迭代并行算法的MP源程序 1.源程序 seidel.c #include"mpi. h #include "math. h
处理器中 end while End 与 并 行 高 斯 - 塞 德 尔 迭 代 相 似 , 并 行 松 弛 迭 代 法 的 一 轮 并 行 计 算 时 间 为 : T = mn + m + n(t + t )log p + 2t ( p −1) + tw(p −1) p s w s 。 MPI 源程序请参见所附光盘。 1.4 小结 本章主要讨论线性方程组的迭代解法，这种方法是一种逐步求精的近似求解过程，其优 点是简单，易于计算机编程，但它存在着迭代是否收敛以及收敛速度快慢的问题。一般迭代 过程由预先给定的精度要求来控制，但由于方程组的准确解一般是不知道的，因此判断某次 迭代是否满足精度要求也是比较困难的，需要根据具体情况而定。文献[1]给出了稀疏线性 方程组迭代解法的详尽描述，还包含了多重网格(Multigrid)法、共轭梯度(Conjugate Gradient) 法，[2]综述了稀疏线性方程组的并行求解算法，[3]综述了在向量机和并行机上偏微分方程 的求解方法，[4]讨论了超立方多处理机上的多重网格算法，[5]讨论了并行共轭梯度算法， [6]深入而全面地论述了 SIMD 和 MIMD 模型上的数值代数、离散变换和卷积、微分方程、 计算数论和最优化计算的并行算法，对并行排序算法也作了介绍。此外，在[7]中第九章的 参考文献注释中还列举了大量有关参考文献，进一步深入研究的读者可在这些文献中获得更 多的资料。 参考文献 [1]. 陈国良 编著．并行计算——结构·算法·编程．高等教育出版社,1999.10 [2]. Heath M T, Ng E and Peyton B W. Parallel Algorithm for Space Linear Systems. SIAM Review,1991,33:420-460 [3]. Ortega J M, Voigt R G. Solution of Partial Differential Equations on Vector and Parallel Computers. SIAM Review,1985,27(2):149-240 [4]. Chan T F,Saad Y. Multigrid Algorithms on the Hypercube Multiprocessor. IEEE-TC,1986,C-35(11):969-977 [5]. 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 [6]. 李晓梅，蒋增荣 等编著．并行算法(第五章) ．湖南科技出版社, 1992 [7]. Quinn M J. Parallel Computing-Theory and Practic(second edition)McGraw-Hill, Inc., 1994 附录 高斯-塞德尔迭代并行算法的 MPI 源程序 1. 源程序 seidel.c #include "stdio.h" #include "stdlib.h" #include "mpi.h" #include "math.h
#define E 0.0001 total=o #define a(x,y)ax*sizetyl MPI Init( &argc, &argv) #define b(x)bx] MPI Comm size(MPI COMM WORld #define v(x)vE /A为size·sie的系数矩阵考 MPI Comm rank(MPI COMM WORLD, #define A(x y)Ax*sizetyl &my rank); #define B(x)Bx P-group size, #define V(x)VIx #define intsize sizeof(int) if(my_ rank==0) #define floatsize sizeof(float) #define charsize sizeof( char) A=fopen(datain. txt,r"); &size, &N) if(size =N-1) float *B printf("the input is wrong"); float *A exit(1); float *V. Int my ran A=(float")malloc(floatsize*size*size) malloc(floatsize*size) FILE fdA. *fdB. fdB 1 ze, i++) void Environment Finalize(float *a float *b, float v) for(=0; j< size, j+ scanf( fdA, %f, A+i*size+j) fscanf(fdA, %f, B+i); for( fscanf(fdA, " %f, V+i); int main(int argc, char **argv) printf( "Input of file I"datain. txt"ln") int i j, my rank, group size; printf("%d\%dn", size, N); int k sum for(i=0 float b for(j=0; j float*a: printf("%ft", A(i)) float "differ printf("%fn", B() float int iteration, total, loop: printf("n"); for(i=0; 1< size, 1++) printf("%ft", vO)
#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("nn"); printf("oUtput of result") se MPI Recv(am*size, MPI FLOAT, 件0号处理器将sie广播给所有处理器 k MPI COMM WORld MPI Bcast(&size, 1, MPI INT,O, &status ); MPI COMM WORLD) MPI Recv(b, m, MPI FLOAT, O, my_ rank, maize/p if (size%p!=0)m++; MPI COMM WORLD, &status); v=(float ")malloc( floatsize* size) 所有处理器并行地对主对角元素右边的数 um(float")malloc(floatsize"m) 据求和* if(aF=NULL b==NULLIV=NULL r(i=0;1<m;计+) printf("allocate space fail! ) if(my rank==0 sum[j=0.0, for(F0 j<size j++) for(F0 i <size: 1++) if ((my rank*m+i) sumO=sum 0+a(ij*v0: 陴初始解向量ⅴ被广播给所有处理器* MPI Bcast(v, size, MPI FLOAT,O, MPI COMM WORLD) iteration 0号处理器采用行块划分将矩阵A划分为大 小为m*sie的p块子矩阵,将B划分为 for(F0 j<size j++) 大小为m的p块子向量,依次发送给1 至p-1号处理器* f (my rank==0 编号为q的处理器负责计算解向 量并广播给所有处理器* if(my for(=0; j<size: j++) a(lfA(1j; temp=v(my rank*m+r) for(F0, i<r; i++ b(l=B(), sum[r=sum[rl+ a(r,my rank"m+i)* v( my rank’m+i; MPI Send(&(a(m*,0)), m*size, /计算出的解向量考 MPI FLOATLI (my rank'mtrF MPI COMM WORLD) MPI Send&(B(m·)m, a(r, my_ rank*m +r); MPI FLOAT. LI if(fabs( v(my rank*m+r) MPI COMM WORLD) temp<e) Iteration++ free(A); free(B); free(V) 广播解向量*
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(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++; /*广播解向量*/