用部分选主元的高斯消去法并行求解线性方程组 精灵论文.doc

《用部分选主元的高斯消去法并行求解线性方程组 精灵论文.doc》由会员分享，可在线阅读，更多相关《用部分选主元的高斯消去法并行求解线性方程组 精灵论文.doc（7页珍藏版）》请在三一办公上搜索。

1、用部分选主元的高斯消去法并行求解线性 方程组刘向娇，刘佳梅（宁夏大学数学计算机学院，银川 750021） 摘要：高斯消去法，又称高斯消元法，实际上就是我们俗称的加减消元法。数学上，高斯消 去法或称高斯约当消去法，由高斯和约当得名(很多人将高斯消去作为完整的高斯约当 消去的前半部分)，它是线性代数中的一个算法，用于决定线性方程组的解，决定矩阵的秩， 以及决定可逆方矩阵的逆。当用于一个矩阵时，高斯消去产生“行消去梯形形式”。用高斯 消去法求解线性方程组的解是一种比较常见的解线性方程组的方法，这种方法尤其在利用计 算机求解线性方程组时是更是常用。但大多数情况下都是用串行的算法来解方程组，本文介 绍了

2、利用高斯消去法并行求解线性方程组的方法。 关键词：高斯消去法；求解线性方程组；并行中图分类号：TP3Portion of pivoting in parallel Gaussian elimination method for solving linear equationsLiu Xiangjiao, Liu Jiamei(Ningxia University,YinChuan 750021)Abstract: Gaussian elimination, also calls Gauss the elimination, in fact is the elimination by addn

3、which we are named. In mathematics, the gaussian elimination or called Gauss - approximately works as elimination, by the Gauss peace treaty, when acquires fame (many people to eliminate Gauss takes complete Gauss - when approximately eliminates the first half part), it is in a linear algebra algori

4、thm, uses in deciding that system of linear equationss solution, decides the matrix the order, as well as decision reversible side matrix going against. When uses in a matrix, Gauss eliminates the production “the line to eliminate the trapezoidal form”. With the gaussian elimination solution system

5、of linearequationss solution is one quite common solution system of linear equationss method, this methodwhen solves the system of linear equations especially using the computer is commonly used. But in the majority situations was solves the equation set with the serial algorithm, this article intro

6、duced used the gaussian elimination parallel solution system of linear equationss method.Keywords:Gaussian elimination; system of linear equations; parallel1高斯消去法在求解线性方程组的算法中，有两类最基本的算法：直接法和迭代法。在直接法中最主 要的就是高斯消去法，它可以分为消元和回代两个过程，消元过程是将方程组转换为一个等 价的三角方程组，而回代过程则是逆反求解这个三角方程组12。2部分选主元的高斯消去思想kk前述的消去过程中，未知量是按

7、其出现于方程组中的自然顺序消去的，所以又叫顺序消 去法。实际上已经发现顺序消去法有很大的缺点。设用作除数的 a ( k 1) 为主元素，首先，消元过程中可能出现 a ( k 1) 为零的情况，此时消元过程亦无法进行下去；其次如果主元素a ( k 1)kk kk作者简介：刘向娇（1984-），女，硕士研究生，主要研究方向：并行处理与高性能计算. E-mail:liuxiangjiao042很小，由于舍入误差和有效位数消失等因素，其本身常常有较大的相对误差，用其作除数3， 会导致其它元素数量级的严重增长和舍入误差的扩散，使得所求的解误差过大，以致失真。我们来看一个例子： 例0.0001x1 + 1

8、.00 x2 = 1.001.00 x1 + 1.00 x2 = 2.00它的精确解为：x = 10000 1.000101 9999x = 9998 0.999902 9999用顺序消去法，第一步以 0.0001 为主元，从第二个方程中消 x1 后可得： 10000 x2= 10000x2 = 1.00回代可得x1 = 0.00 显然，这不是解。 造成这个现象的原因是：第一步主元素太小，使得消元后所得的三角形方程组很不准确所致。如果我们选第二个方程中 x1 的系数 1.00 为主元素来消去第一个方程中的 x1，则得出如 下方程式：1.00 x1 = 1.00x1 = 1.00 这是真解的三位

9、正确舍入值。 从上述例子中可以看出，在消元过程中适当选取主元素是十分必要的。误差分析的理论和计算实践均表明：顺序消元法在系数矩阵 A 为对称正定时，可以保证此过程对舍入误差 的数值稳定性，对一般的矩阵则必须引入选取主元素的技巧，方能得到满意的结果。高斯消去法的一个简单变形部分选主元高斯消去法，可以产生可靠的结果。在部分 选主元的高斯消去法的第 i 步，我们在第 i 行到第 n-1 行中寻找第 i 列元素的绝对值最大的 行并将这一行与第 i 行交换（变成主元）1。在部分主元消去法中，未知数仍然是顺序地消去的，只是选各方程中要消去的那个未知 数的系数的绝对值最大的作为主元素，然后用顺序消去法的公式

10、求解。例：用部分选主元高斯消去法求解方程2 x1 x2 + 3x3 = 14 x1 + 2 x2 + 5x3 = 4 x1 + 2 x2 = 7由于解方程组取决于它的系数，因此可用这些系数（包括右端项）所构成的“增广矩阵” 作为方程组的一种简化形式。对这种增广矩阵施行消元手续： 2 4*1 1 312 5 4 72 0 第一步将 4 选为主元素，并把主元素所在的行定为主元行，然后将主元行换到第一行得 到 42754 10.561.251 2 131 第一步消元 0 2*0.5 112010.51.25 01.51 1.2510.51.251 第二步消元 01 0.250.5 第三步消元 01

11、0.250.5 00 0.8755.25 001 6 消元过程的结果归结到下列三角形方程组： x1 + 0.5x2 + 1.25x3 = 1回代，得x2 0.25x3x3= 0.5= 6 x1 x2= 9= 1 x3= 63串行部分选主元高斯消去算法先部分选主元然后回代的串行高斯消去算法如下所示：for i0 to n-1 loci i endforfor i 0 to n-1找到主元行 pickedmagnitude 0 for j i to n-1if|alocj，i|magnitudemagnitude |alocj，i|picked j endif endfortmp lociloci

12、 locpickedlocpicked tmp将第 i 列中位于未标记行中的元素消为 0for j i+1 to n-1t alocj，i/aloci，ifor k i+1 to n+1alocj，k alocj，k-aloci，k t endforendfor endfor回代for in-1 down to 0xi aloci，n/aloci，ifor j0 to i-1 doalocj，n alocj，n-xi alocj，iendfor endfor该算法有两个值得注意的特性：1没有独立的数组来存储向量 b，而是用 b 邻接 A 创建一个 n 行 n+1 列的增广矩阵。 因此，在这个算

13、法中 a 代表这个增广矩阵。2算法每次迭代中间接地访问主元行，而不是直接地将主元行和第 i 行相交换。数组 元素 loci记录了第 i 次迭代中所用的主元行的下标4。4并行部分选主元高斯消去算法4.1 算法思想在上面的串行算法中位于最内层下标为 k 的 for 循环和位于中间层的下标为 j 的 for 循 环都可以拿来并行执行，换句话说，一旦找到主元行，对所有没被标记的行的修改就可以同 时进行。在每一行内，一旦系数 alocj，i/aloci，i被计算出来了，那么这一行内从位置 i+1 开始的 n-1 个元素的修改工作也可以同时进行。在并行化的过程中，由于高斯消去法是利用主元行 i 对其余各行

14、 j(ji)，作初等行变换， 各行计算之间没有数据相关关系，因此可以对矩阵 A 按行划分。考虑到在计算过程中处理 器之间的负载均衡，对 A 采用行交叉划分。设处理器个数为 p，矩阵 A 的阶数为 n，m= n / p ， 对矩阵 A 行交叉划分后，编号为 i(i=0，1，p-1)的处理器含有 A 的第 i，i+p，i+(m-1)p 行和向量 B 的第 i，i+p，i+(m-1)p 一共 m 个元素。矩阵划分好以后，接下来的主要工作是寻找主元行。我们发现为了决定主元行，需要对 分布在不同处理器中的第 i 列的各个元素进行规约。我们可以把这一过程称为巡回赛，因为与主元行的第 i 列中存储的具体数值

15、(胜利者的得分)相比，我们对主元行的位置(胜利者的身份)更感兴趣。在第 i 次迭代进行的过程中，主元行的确定需要两个步骤。首先，每个进程在它所负责 的未标记的行中寻找在第 i 列数值最大的那一行。第二，进程加入寻找主元行的巡回赛5。 在每次迭代过程中，这一步都是首先要完成的，之后还有一个涉及通信的任务，如图 1所示。为了计算 aj，k的新值，其对应的任务需要访问 aj，i、apicked，i和 apicked，k 的值。每个任务至少分配到 A 的一行，所以这个任务既然控制了 aj，k，那么也就控制着 aj，i。但是 apicked，i和 apicked，k的值可能为另外一个任务所控制。因此，我

16、们还需 要做一次广播3。图 1 迭代时通信任务图4.2 算法实现输入：系数矩阵 Ann，常数向量 bn1输出：解向量 xn1Begin对所有处理器 my_rank(my_rank=0，p-1)同时执行如下的算法：/*消去过程*/(1) for i=0 to m-1 do for j=0 to p-1 do(1.1) if(my_ranklmax0) thenlmax0=ak，t，lmax1=k，lmax2=t，lmax3=my_rankend if end forend forend if(1.2) if (my_rankj) then(i) v=i*p+j(ii) lmax0=ai，v(ii

17、i) for k=i to m-1 do for t=v to n-1 doif(|ak，t|)lmax0) thenlmax0=ak，t，lmax1=k，lmax2=t，lmax3=my_rankend if end forend forend if(1.3) 将每一个处理器中的 lmax 元素广播到其他所有处理器中(1.4) maxvalve=getpivort(max)，maxrow=getrow(max)，maxrank=getrank(max) (1.5) if(my_rank=j) then(i) if(maxrank=j and imaxrow) theninnerexchang

18、erow()end if(ii) if (maxrankj) thenouterexchangerow()end if end ifend if(1.6) if(my_rank=j) then(i) for k=v+1 to n-1 doai，k=ai，k/ai，vend for(ii) bi=bi/ai，v，ai，v=1(iii) for k=v+1 to n-1 dofk=ai，kend for(iv) fn=bi(v) 广播主行到所有处理器else 接收广播来的主行元素存与数组 f 中end if(1.7) if(my_rankj) thenfor k=i+1 to m-1 do(i)

19、for w=v+1 to n-1 doak，w=ak，w-fw*ak，vend for(ii) bk=bk-fn*ak，vend for end if(1.8) if(my_rankj) thenfor k=i to m-1 do(i) for w=v+1 to n-1 doak，w=ak，w-fw*ak，vend for(ii) bk=bk-fn*ak，vend for end ifend for end for/*回代过程*/(2) for i=0 to m-1 dosumi=0.0end for(3) for i=m-1 downto 0 dofor j=p-1 downto 0 doi

20、f(my_rank=j) then(i) xi*p+j=(bi-sumi)/ai，i*p+j(ii) 将 xi*p+j广播到所有处理器中(iii) for k=0 to i-1 dosumk=sumk+ak，i*p+j* xi*p+jend forelse /*非主行所在的处理器*/(iv) 接收广播来的 xi*p+j的值(v) if(my_rankj) then for k=0 to i-1 dosumk=sumk+ak，i*p+j* xi*p+jend for end if(vi) if(my_rankj) thenfor k=0 to i dosumk=sumk+ak，i*p+j* xi

21、*p+jend for end ifEndend if end forend for4.3 算法分析在第 i 次迭代主元行确定的两个步骤中，每个进程在它所负责的未标记的行中寻找在第 i 列中最大的那一行的时间复杂度是 O(n/p)，而进程加入寻找主元行的巡回赛的时间复杂度 是 O(logp)。综合考虑通信的总开销，我们发现面向行的部分选主元的高斯消去算法的总消息延迟是O(nlogp)，而总的消息传输时间是 O(n2logp)6。参考文献 (References)1 陈国良 安虹 陈崚等.并行算法实践M.高等教育出版社.2004 年 1 月.第九章 P504-511 2 胡哲光.在 C+中实现带主元选择的高斯消去法求解线性方程J.大众科技.2006 年第八期 3 陈文光译.MPI 与 Open MP 并行程序设计(C 语言版)M.清华大学出版社.2004 年 10 月4 陈国良等.并行算法实践M.高等教育出版社.2004 年 1 月5 陆鑫达等译.并行程序设计(第 2 版)M.机械工业出版社.2005 年 5 月6 张云泉 陈英译.并行算法导论M.机械工业出版社.2004 年 2 月