姬 進(jìn)，艾 莉，郭警濤，趙 君，閆 穩(wěn)

(航空工業(yè)西安航空計(jì)算技術(shù)研究所，陜西 西安 710065)

并行計(jì)算機(jī)不僅僅是一種獲得高性能的手段，它同時(shí)也具有將計(jì)算能力從單個(gè)處理器無(wú)縫地?cái)U(kuò)展到無(wú)限多個(gè)處理器的潛力。這種潛力在幾十年前就已被人們所熟知，但是直到20世紀(jì)80年代末，它才真正顯現(xiàn)出來(lái)。

本文先對(duì)并行計(jì)算的一些基本概念作個(gè)簡(jiǎn)單介紹，然后介紹線性方程組的兩種迭代法的并行算法的原理和實(shí)現(xiàn)，最后對(duì)這兩種算法作出相應(yīng)的測(cè)試和性能分析。

1 線性方程組的并行計(jì)算簡(jiǎn)述

并行計(jì)算就是在并行計(jì)算機(jī)上，將一個(gè)應(yīng)用分解成多個(gè)子任務(wù)，分配給不同的處理器，各個(gè)處理器之間相互協(xié)同，并行地執(zhí)行子任務(wù)，從而達(dá)到加快求解速度，或者提高求解應(yīng)用問(wèn)題規(guī)模的目的[1-2]。

線性方程組是指矩陣含有大量零元素的線性代數(shù)方程組。矩陣多種多樣，主要可以分為具有規(guī)則結(jié)構(gòu)的結(jié)構(gòu)和無(wú)規(guī)則結(jié)構(gòu)的矩陣基本類型。

求解線性方程組的算法中，有兩類最基本的算法，一類是直接法，即以消去為基礎(chǔ)的解法，如果不考慮誤差的影響，從理論上講，它可以在固定步數(shù)內(nèi)求得方程組的準(zhǔn)確解；另一類是迭代解法，它是一個(gè)逐步迭代求得近似解的過(guò)程，這種方法便于編制程序，但存在著迭代是否收斂及收斂速度快慢的問(wèn)題；在迭代過(guò)程中，由于極限過(guò)程一般不可能進(jìn)行到底，因此只能得到滿足一定精度要求的近似解。

2 迭代法介紹

在科學(xué)工程領(lǐng)域中，經(jīng)常需要求解形如Ax=b的線性方程組，其中A為n×m的矩陣，稱為系數(shù)矩陣；b為n維向量，稱為列向量；x為待求解的m維列向量，常簡(jiǎn)稱解向量。在本論文中僅考慮m=n的情形。

迭代法的一般形式化可以描述為：

x(k)=φk(x(k-1),…,x(k-l)),k=l,l+1,…, .

(1)

其中φk稱為迭代算子，x(0),x(1),…，x(l-1)稱為初值。通常稱上式為l步迭代法，當(dāng)l=1時(shí)稱為單步迭代法。如果算子φk與k無(wú)關(guān)，則稱為定常迭代，否則為非定常迭代。如果φk為線性算子，則稱為線性迭代，否則為非線性迭代。

在求解線性方程組的迭代法中，單步定常線性迭代扮演著重要的角色，該迭代可以描述為：

x(k)Gx(k-1)+c,k=1,2,3,…，.

(2)

其中x(k)與c為n維向量，G為n階矩陣。如果存在x，使得對(duì)任意初值x(0)，由上式產(chǎn)生的序列{x(k)}都收斂于x，則稱該迭代法收斂，且記1imx(k)=x，否則發(fā)散。

3 迭代法并行算法設(shè)計(jì)

3.1 雅克比迭代法

對(duì)于方程組Ax=b，記D,-L,-U分別是A的對(duì)角、嚴(yán)格下三角、嚴(yán)格上三角部分構(gòu)成的矩陣，即A=D-L-U。這時(shí)方程組(1)式可以變成：

Dx=b+(L+U)x.

(3)

如果方程組右邊的x已知，由于D是對(duì)角矩陣，可以很容易求得左邊的x，這就是雅克比(Jacobi)迭代的出發(fā)點(diǎn)。因此，對(duì)于給定的初值x(0)，雅克比迭代法如下：

x(k)=D-1(L+U)x(k-1)+D-1b.

(4)

記G=D-1(L+U)=I-D-1A，c=D-1b.

對(duì)于方程組Ax=b式，方程組的每個(gè)方程ailx1+ai2x2+…+ainxn=bi可以很容易地寫(xiě)成如下的形式：

(5)

(6)

設(shè)有p個(gè)處理器，則將矩陣劃分為p塊，每塊含有連續(xù)的m行向量，這里m=[n/p]，編號(hào)為i的處理器含有A的第im至第(i+1)m-1行數(shù)據(jù)，同時(shí)向量b中下標(biāo)為im至(i+1)m-1的元素也被分配至編號(hào)為i的處理器(i=0,1,…,p-1),初始解向量被廣播給所有的處理器。

在迭代計(jì)算中，各處理器并行計(jì)算解向量x的各分向量值，編號(hào)i的處理器計(jì)算分量xim至x(i+1)m-1的新值。并求其分量中前后兩次值的最大差localmax，然后通過(guò)歸約操作的求最大值運(yùn)算，求得整個(gè)n維解向量中前后兩次值的最大差值max，并廣播給所有的處理器。最后通過(guò)擴(kuò)展收集操作將各處理器中的解向量按處理器編號(hào)連接起來(lái)，并廣播給所有的處理器，以進(jìn)行下一次迭代計(jì)算，直至收斂。具體算法描述見(jiàn)算法1。

算法1 線性方程組的雅克比迭代并行算法

max=ε+1；

While max>ε do

Fori=0 tom-1 do

sum=0.0；

Forj=0 ton-1 do

Ifj≠myid×m+ithen

sum=sum+aijxi；

End if

End for

End for

Fori=1 tom-1 do

End if

End for

將x的所有分量的新值廣播到所有處理器；

求出所有處理器中l(wèi)ocalmax值的最大值max,

并廣播給所有的處理器中；

End while

3.2 高斯-賽德?tīng)柕?/h3>

在并行計(jì)算中采用與雅克比迭代類似的數(shù)據(jù)劃分。對(duì)于高斯-賽德?tīng)柕?，?jì)算xi的新值時(shí)，使用xi+1,…，xn-1的舊值和x0,…,xi-1的新值。計(jì)算過(guò)程中xi與x0,…,xi-1及xi+1,…,xn-1的新值會(huì)在不同的處理器中產(chǎn)生，因此可以考慮采用時(shí)間偏移的方法，使各個(gè)處理器對(duì)新值計(jì)算的開(kāi)始和結(jié)束時(shí)間產(chǎn)生一定的偏差。編號(hào)為myid的處理器一旦計(jì)算出xi(myid×m≤i≤(myid+1)×m)的新值，就立即廣播給其余處理器，以供各處理器對(duì)x的其他分量計(jì)算有關(guān)xi的乘積項(xiàng)并求和。當(dāng)它計(jì)算完x的所有分量后，它還要接收其他處理器發(fā)送的新的x分量，并對(duì)這些分量進(jìn)行求和計(jì)算，為計(jì)算下一輪的xi作準(zhǔn)備。計(jì)算開(kāi)始時(shí)，所有處理器并行地對(duì)主對(duì)角元素右邊的數(shù)據(jù)項(xiàng)進(jìn)行求和，此時(shí)編號(hào)為0的處理器(簡(jiǎn)稱為P0)計(jì)算出x0，然后廣播給其余處理器，其余所有的處理器用x0的新值和其對(duì)應(yīng)項(xiàng)進(jìn)行求和計(jì)算，接著P0計(jì)算出x1,x2,…,當(dāng)P0完成對(duì)xm-1的計(jì)算和廣播后，P1計(jì)算出xm，并廣播給其余處理器，其余所有的處理器用xm的新值求其對(duì)應(yīng)項(xiàng)的乘積并作求和計(jì)算。然后P1計(jì)算出xm+1,xm+2,…,當(dāng)P1完成對(duì)x2×m-1的計(jì)算和廣播后，P2計(jì)算出x2×m…，如此重復(fù)下去，直至xn-1在PP-1中被計(jì)算出并廣播至其余的處理器之后，P0計(jì)算出下一輪新的x0，這樣逐次迭代下去，直至收斂為止。具體算法描述見(jiàn)算法2。

算法2 高斯-賽德?tīng)柕⑿兴惴?/p>

Fori=myid×mto (myid+1)×m-1 do

sumi=0.0；

Forj=i+1 ton-1 do

sumi=sumi+aijxi；

End for

End for

total=0；

Whiletotal

iteration=0；

Forj=0 ton-1 do

q=j/m；

Ifmyid==qthen

temp=xj,xj=(bi-sumi)/aii；

If |xi-temp|<ε then

iteration=iteration+1；

End if

將xj的新值廣播到其他所有處理器；

sumi=0；

Fori=myid×mto (myid+1)×m-1 do

Ifj≠ithen

sumi=sumi+aijxi；

End if

End for

Else

接受廣播來(lái)的xj的新值；

Fori=myid×mto (myid+1)×m-1 do

sumi=sumi+aijxi；

End for

End if

End for

求出所有處理器中iteration值的和total并廣播到所有處理器；

End while

4 算法測(cè)試與性能分析

假設(shè)某個(gè)串行應(yīng)用程序在某臺(tái)并行計(jì)算機(jī)上單處理器上的執(zhí)行時(shí)間為T(mén)s，而該程序并行化后，P個(gè)進(jìn)程在P個(gè)處理器并行執(zhí)行所需要的時(shí)間為T(mén)p，則該并行程序在該計(jì)算機(jī)上的加速比可定義為：

(7)

效率定義為：

(8)

加速比和效率是衡量一個(gè)并行程序性能的最基本評(píng)價(jià)方法，顯然執(zhí)行最慢的進(jìn)程將決定并行程序的性能。

4.1 雅克比迭代法性能分析

該算法的時(shí)間復(fù)雜度為O(mm)，其中m為每個(gè)進(jìn)程中矩陣的行數(shù)。在單進(jìn)程計(jì)算時(shí)，m=n，此時(shí)的時(shí)間復(fù)雜度為O(n2)，從表1中可以看到執(zhí)行時(shí)間隨著規(guī)模的增大也相應(yīng)增大了n2倍。

表1 單進(jìn)程雅克比迭代并行算法運(yùn)行時(shí)間

當(dāng)求解規(guī)模在10000時(shí)，其對(duì)應(yīng)的單進(jìn)程執(zhí)行時(shí)間為72.340110 s，兩個(gè)進(jìn)程并行執(zhí)行時(shí)間為25.584798 s，可以求得其加速比S2=2.82，效率E2=1.41，四個(gè)進(jìn)程并行執(zhí)行時(shí)間為13.505695 s，易求得加速比S4=5.50，效率E4=1.375。兩個(gè)或四個(gè)進(jìn)程并行后，迭代次數(shù)減少了(從表2和表3中可以看到)，所以其并行后計(jì)算時(shí)間明顯減少，從而導(dǎo)致加速比明顯增大，也導(dǎo)致該并行算法的效率很高。

表2 兩個(gè)進(jìn)程雅克比迭代并行算法運(yùn)行時(shí)間

表3 四個(gè)進(jìn)程雅克比迭代并行算法運(yùn)行時(shí)間

4.2 高斯-賽德?tīng)柕阅芊治?/h3>

通過(guò)表4和表5可看到：高斯-賽德?tīng)柕牟⑿行暂^差，當(dāng)增加進(jìn)程數(shù)量時(shí)，時(shí)間并沒(méi)有減少，這是由于高斯-賽德?tīng)柕看我们懊嬗?jì)算的新值，因此算法會(huì)有一定的等待時(shí)間，導(dǎo)致整體性能下降。

表4 兩進(jìn)程高斯-賽德?tīng)柕⑿杏?jì)算運(yùn)行時(shí)間

表5 四進(jìn)程高斯-賽德?tīng)柕⑿杏?jì)算運(yùn)行時(shí)間

5 結(jié)束語(yǔ)

線性方程組的迭代求解是一種逐步求得近似解的過(guò)程，這種方法的好處是方法簡(jiǎn)單，便于編制程序，但是迭代法必須符合合適的迭代格式及初始向量，以使得迭代過(guò)程盡快收斂。本文介紹了兩種經(jīng)典的迭代算法——雅克比迭代和高斯-賽德?tīng)柕⒃O(shè)計(jì)了其并行算法。雅克比迭代有著內(nèi)在的并行性，因此它的并行性較高斯-賽德?tīng)柕^好，但是高斯-賽德?tīng)柕蒙洗巫钚碌闹?，因此在并行性上有所欠缺，但是它的收斂速度卻提高了。