用部分選主元的高斯消去法并行求解線性方程組 精靈論文

1、用部分選主元的高斯消去法并行求解線性 方程組劉向嬌，劉佳梅（寧夏大學(xué)數(shù)學(xué)計算機(jī)學(xué)院，銀川 750021） 摘要：高斯消去法，又稱高斯消元法，實(shí)際上就是我們俗稱的加減消元法。數(shù)學(xué)上，高斯消 去法或稱高斯約當(dāng)消去法，由高斯和約當(dāng)?shù)妹?很多人將高斯消去作為完整的高斯約當(dāng) 消去的前半部分)，它是線性代數(shù)中的一個算法，用于決定線性方程組的解，決定矩陣的秩， 以及決定可逆方矩陣的逆。當(dāng)用于一個矩陣時，高斯消去產(chǎn)生“行消去梯形形式”。用高斯 消去法求解線性方程組的解是一種比較常見的解線性方程組的方法，這種方法尤其在利用計 算機(jī)求解線性方程組時是更是常用。但大多數(shù)情況下都是用串行的算法來解方程組，本文介 紹了

2、利用高斯消去法并行求解線性方程組的方法。 關(guān)鍵詞：高斯消去法；求解線性方程組；并行中圖分類號：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高斯消去法在求解線性方程組的算法中，有兩類最基本的算法：直接法和迭代法。在直接法中最主 要的就是高斯消去法，它可以分為消元和回代兩個過程，消元過程是將方程組轉(zhuǎn)換為一個等 價的三角方程組，而回代過程則是逆反求解這個三角方程組12。2部分選主元的高斯消去思想kk前述的消去過程中，未知量是按

7、其出現(xiàn)于方程組中的自然順序消去的，所以又叫順序消 去法。實(shí)際上已經(jīng)發(fā)現(xiàn)順序消去法有很大的缺點(diǎn)。設(shè)用作除數(shù)的 a ( k 1) 為主元素，首先，消元過程中可能出現(xiàn) a ( k 1) 為零的情況，此時消元過程亦無法進(jìn)行下去；其次如果主元素a ( k 1)kk kk作者簡介：劉向嬌（1984-），女，碩士研究生，主要研究方向：并行處理與高性能計算. e-mail:liuxiangjiao042很小，由于舍入誤差和有效位數(shù)消失等因素，其本身常常有較大的相對誤差，用其作除數(shù)3， 會導(dǎo)致其它元素數(shù)量級的嚴(yán)重增長和舍入誤差的擴(kuò)散，使得所求的解誤差過大，以致失真。我們來看一個例子： 例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 顯然，這不是解。 造成這個現(xiàn)象的原因是：第一步主元素太小，使得消元后所得的三角形方程組很不準(zhǔn)確所致。如果我們選第二個方程中 x1 的系數(shù) 1.00 為主元素來消去第一個方程中的 x1，則得出如 下方程式：1.00 x1 = 1.00x1 = 1.00 這是真解的三位

9、正確舍入值。 從上述例子中可以看出，在消元過程中適當(dāng)選取主元素是十分必要的。誤差分析的理論和計算實(shí)踐均表明：順序消元法在系數(shù)矩陣 a 為對稱正定時，可以保證此過程對舍入誤差 的數(shù)值穩(wěn)定性，對一般的矩陣則必須引入選取主元素的技巧，方能得到滿意的結(jié)果。高斯消去法的一個簡單變形部分選主元高斯消去法，可以產(chǎn)生可靠的結(jié)果。在部分 選主元的高斯消去法的第 i 步，我們在第 i 行到第 n-1 行中尋找第 i 列元素的絕對值最大的 行并將這一行與第 i 行交換（變成主元）1。在部分主元消去法中，未知數(shù)仍然是順序地消去的，只是選各方程中要消去的那個未知 數(shù)的系數(shù)的絕對值最大的作為主元素，然后用順序消去法的公式

10、求解。例：用部分選主元高斯消去法求解方程2 x1 x2 + 3x3 = 14 x1 + 2 x2 + 5x3 = 4 x1 + 2 x2 = 7由于解方程組取決于它的系數(shù)，因此可用這些系數(shù)（包括右端項(xiàng)）所構(gòu)成的“增廣矩陣” 作為方程組的一種簡化形式。對這種增廣矩陣施行消元手續(xù)： 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 消元過程的結(jié)果歸結(jié)到下列三角形方程組： 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 列中位于未標(biāo)記行中的元素消為 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沒有獨(dú)立的數(shù)組來存儲向量 b，而是用 b 鄰接 a 創(chuàng)建一個 n 行 n+1 列的增廣矩陣。 因此，在這個算

13、法中 a 代表這個增廣矩陣。2算法每次迭代中間接地訪問主元行，而不是直接地將主元行和第 i 行相交換。數(shù)組 元素 loci記錄了第 i 次迭代中所用的主元行的下標(biāo)4。4并行部分選主元高斯消去算法4.1 算法思想在上面的串行算法中位于最內(nèi)層下標(biāo)為 k 的 for 循環(huán)和位于中間層的下標(biāo)為 j 的 for 循 環(huán)都可以拿來并行執(zhí)行，換句話說，一旦找到主元行，對所有沒被標(biāo)記的行的修改就可以同 時進(jìn)行。在每一行內(nèi)，一旦系數(shù) alocj，i/aloci，i被計算出來了，那么這一行內(nèi)從位置 i+1 開始的 n-1 個元素的修改工作也可以同時進(jìn)行。在并行化的過程中，由于高斯消去法是利用主元行 i 對其余各行

14、 j(ji)，作初等行變換， 各行計算之間沒有數(shù)據(jù)相關(guān)關(guān)系，因此可以對矩陣 a 按行劃分?？紤]到在計算過程中處理 器之間的負(fù)載均衡，對 a 采用行交叉劃分。設(shè)處理器個數(shù)為 p，矩陣 a 的階數(shù)為 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 個元素。矩陣劃分好以后，接下來的主要工作是尋找主元行。我們發(fā)現(xiàn)為了決定主元行，需要對 分布在不同處理器中的第 i 列的各個元素進(jìn)行規(guī)約。我們可以把這一過程稱為巡回賽，因?yàn)榕c主元行的第 i 列中存儲的具體數(shù)值

15、(勝利者的得分)相比，我們對主元行的位置(勝利者的身份)更感興趣。在第 i 次迭代進(jìn)行的過程中，主元行的確定需要兩個步驟。首先，每個進(jìn)程在它所負(fù)責(zé) 的未標(biāo)記的行中尋找在第 i 列數(shù)值最大的那一行。第二，進(jìn)程加入尋找主元行的巡回賽5。 在每次迭代過程中，這一步都是首先要完成的，之后還有一個涉及通信的任務(wù)，如圖 1所示。為了計算 aj，k的新值，其對應(yīng)的任務(wù)需要訪問 aj，i、apicked，i和 apicked，k 的值。每個任務(wù)至少分配到 a 的一行，所以這個任務(wù)既然控制了 aj，k，那么也就控制著 aj，i。但是 apicked，i和 apicked，k的值可能為另外一個任務(wù)所控制。因此，我

16、們還需 要做一次廣播3。圖 1 迭代時通信任務(wù)圖4.2 算法實(shí)現(xiàn)輸入：系數(shù)矩陣 ann，常數(shù)向量 bn1輸出：解向量 xn1begin對所有處理器 my_rank(my_rank=0，p-1)同時執(zhí)行如下的算法：/*消去過程*/(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 接收廣播來的主行元素存與數(shù)組 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 doif(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) 接收