1、1. 8 线性方程组的迭代解法 在阶数较大、系数阵为稀疏阵的情况下，可以采用迭代法求解线性方程组。用 迭代法(Iterative Method)求解线性方程组的优点是方法简单，便于编制计算机程序，但必须选取合适的迭代格式及初始向量，以使迭代过程尽快地收敛。迭代法根据迭代格式的不同分成 雅可比 (Jacobi)迭代、 高斯 -塞德尔 (Gauss-Seidel)迭代和 松弛 (Relaxation)法等几种。在本节中，我们假设系数矩阵 A 的主对角线元素 0iia ，且按行严格 对角占优 (Diagonal Donimant)，即： ), . . . ,2,1(1niaa nj ijijii 1
2、.1 雅可比迭代 1.1.1 雅可比迭代 及其串行算法 雅可比迭代的原理是：对于求解 n 阶线性方程组 Ax=b，将原方程组的每一个方程 ai1x1+ ai2x2+ + ainxn= bi 改写为未知向量 x 的分量的形式： )1(/)( ,1 niaxabx iin ijj jijii 然后使用第 k-1 步所计算的变量 xi(k -1)来计算第 k 步的 xi(k)的值： ),1(/)( ,1 )1()( nkiaxabx iin ijj kijik ji 这里， xi(k)为第 k 次迭代得到的近似解向量 x(k)= (x1 (k), x2(k) , , xn(k) )T 的第 i 个分
3、量。取适当初始解向量 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 -xx 1ni1max 判断迭代是否收敛。如果取一次乘法和加法运 算时间或一次比较运算时间为一个单位时间，则下述雅可比迭代算法 20.1 的一轮计算时间为 n2+n=O(n2)。 算法 20
4、.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)f
5、or 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 行向量，这里 pnm / ，编号为 i 的处理器含有 A 的第 im 至第 (i+1)m-1 行数据，同时向量 b 中下标为 im 至 (i+1)m-1 的元素也被分配至编号为 i 的处理器 (i=0,1,
6、 ,p-1)，初始解向量 x 被广播给所有处理器。 在迭代计算中，各处理器并行计算解向量 x 的各分量值，编号为 i 的处理器计算分量xim 至 x(i+1)m-1 的新值。并求其分量中前后两次值的最大差 localmax，然后通过归约操作的求最大值运算求得整个 n 维解向量中前后 两次值的最大差 max 并广播给所有处理器。最后通过扩展收集操作将各处理器中的解向量按处理器编号连接起来并广播给所有处理器，以进行下一次迭代计算，直至收敛。具体算法框架描述如下： 算法 20.2 求解线性方程组的雅可比迭代并行算法 输入： 系数矩阵 An n，常数向量 b n 1， ，初始解向量 xn 1 输出：
7、解向量 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+ai,j*xj end if end for (1.3)x1i=(bi - sum)/ai,my_rank*m+i end for (2)/*求出本处理器计算的 x 的相应的分量的新值与原值的差的最大值 loca
8、lmax */ localmax=x10-x0 (3)for i=1 to m-1 do if (x1i-xi localmax) then localmax =x1i-xi end if end for (4)用 Allgather 操作将 x 的所有分量的新值广播到所有处理器中 (5)用 Allreduce 操作求出所有处理器中 localmax 值的最大值 max 并广播到所有处理器中 end while End 若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则一轮迭代的计算时间为 mn+m；另外，各处理器在迭代中做一次归约操作，通信量为 1，一次扩展收集操作，通信量为 m
9、，需要 的通信时间为 )1()1()1(4 ptmpt ws ，因此算法 20.2 的一轮并行计算时间为 mmnptmptT wsp )1()1()1(4 。 MPI 源程序请参见所附光盘。 1.2 高斯 -塞德尔迭代 1.2.1 高斯 -塞德尔迭代及其串行算法 高斯 -塞德尔迭代的基本思想与雅可比迭代相似。它们的区别在于，在雅可比迭代中，每次迭代时只用到前一次的迭代值，而在高斯 -塞德尔迭代中，每次迭代时充分利用最新的迭代值。一旦一个分量的新值被求出，就立即用于后续分量的迭代计算，而不必等到所有分量的新值被求出以后。设方程组 Ax=b 的第 i 个方程为： nj1 ija jx = ib )
10、,2,1( ni 高斯 -塞德尔迭代公式为： )(1 1 )(11 )1()1( nij kjijij kjijiiiki xaxabax ),2,1( ni 取适当的 x(0)作为初始向量，由上述迭代格式可得出近似解向量 x(k)。若原方程组的系数矩阵是按行严格对角占优的，则 x(k)收敛于方程组 的解 x，若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间，则下述高斯 -塞德尔迭代算法 20.3 的一轮计算时间为n2+n=O(n2)。 算法 20.3 单处理器上求解线性方程组的高斯 -塞德尔迭代算法 输入： 系数矩阵 An n，常数向量 b n 1， ，初始解向量 xn 1 输出：
11、 解向量 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-tp) then p=xi-tend if end for end while End 1.2.2 高斯 -塞德尔迭代并行算法 在并行计算中，高斯 -塞德尔迭代采用与雅可比迭代相同的数据划分。
12、对于高斯 -塞德尔迭代，计算 xi的新值时，使用 xi+1, ,xn-1的旧值和 x0, ,xi-1的新值。计算过程中 xi与 x0, ,xi-1及 xi+1, ,xn-1 的新值会在不同的处理器中产生，因此可以考虑采用时间偏移的方法，使各个处理器对新值计算的开始和结束时间产生一定的偏差。编号为 my_rank 的处理器一旦计算出 xi(my_rank m i 1 时，称之为超松弛法，此时， x i(k +1)的比重被加大；当 w=1 时，它就成了一般的高斯 -塞德尔迭代。实际应用时，可根据系数矩阵的性质及对其反复计算的经验来决定适合的松弛因子 w。若取一次乘法和加法运算时间或一次比较运算时间
13、为一个单位时间，则下述算法 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 +
14、w(bi-s)/aii (iv) if (yi - xip) then p=yi - xiend if (v) xi =yi end for end while End 1.3.2 松弛法并行算法 松弛法并行算法和高斯 -塞德尔迭代并行算法基 本相同，只是在计算分量的时候，将对 x分量的新值的计算改为“ xj=(1-w)*temp+w* (bj- sumj)/aj,j”，其中 temp 为 xj的原有值。具体并行算法描述如下 : 算法 20.6 求解线性方程组的松弛迭代并行算法 输入： 系数矩阵 An n，常数向量 b n 1， ，初始解向量 xn 1 输出： 解向量 xn 1 Begin 对
15、所有处理器 my_rank(my_rank=0, , p-1)同时执行如下的算法 : /*所有处理器并行地对主对角元素右边的数据求和 */ (1) for i=my-rank* m to (my-rank+1)*m-1 do (1.1)sumi=0.0 (1.2)for j=i+1 to n-1 do sumi=sumi+ai,j*xj end for end for (2) while (total(my_rank*m+i) sumi=sumi+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+) sumr=sumr+ a(r,my_rank*m+i)*v(my_rank*m+i); /*计算出的解向量 */ v(my_rank*m+r)= (b(r)-sumr)/ a(r,my_rank*m +r); if (fabs(v(my_rank*m+r) -temp)E) iteration+; /*广播解向量 */
