朱京喬,趙永華

1(中國(guó)科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100190)

2(中國(guó)科學(xué)院大學(xué),北京 100049)

對(duì)稱(chēng)矩陣的特征求解問(wèn)題在科學(xué)計(jì)算領(lǐng)域普遍存在,包括計(jì)算物理、計(jì)算化學(xué)、結(jié)構(gòu)力學(xué)、納米材料等前沿學(xué)科.通常的處理方法是先通過(guò)正交變換把原矩陣化為三對(duì)角矩陣,如使用Householder 或Givens方法,再求解三對(duì)角陣的特征值和特征向量,然后映射回原矩陣的特征值和特征向量.三對(duì)角陣的特征求解算法有很多種,常見(jiàn)的有QR 法、二分法、MRRR 方法和分而治之方法等等.由于QR 法和二分法求解特征向量時(shí)往往需要顯式的正交過(guò)程,時(shí)間復(fù)雜度為O(n3),常用于中小規(guī)模矩陣的特征求解.而分治法和MRRR 法能避免這一過(guò)程,分治法的平均時(shí)間復(fù)雜度僅為O(n2.3),其分而治之思想對(duì)于大規(guī)模矩陣特征問(wèn)題的并行求解求解有著天然的優(yōu)勢(shì).

分而治之求解對(duì)稱(chēng)三對(duì)角矩陣特征問(wèn)題的算法最早由Cuppen 提出[1],后M.Gu和S.C.Eisenstat 等人作出改進(jìn),提高了該算法求解的特征向量的正交性[2],使得該算法能夠在實(shí)際求解中應(yīng)用.該算法采用分而治之思想,將原矩陣劃分為若干子矩陣,先求出子矩陣的特征分解,再通過(guò)修正計(jì)算,逐級(jí)合并子矩陣的結(jié)果,回代得到原矩陣的特征分解.分治法具有很強(qiáng)的靈活性,適合大規(guī)模三對(duì)角矩陣特征問(wèn)題求解的并行化求解.

分治算法的并行化工作很早就開(kāi)始展開(kāi)[3-5],如基于分布式存儲(chǔ)的MPI 進(jìn)程級(jí)并行,基于共享內(nèi)存的openMP 線(xiàn)程級(jí)并行,還有基于SMP 的MPI/openMP混合并行等模型.并行化的實(shí)現(xiàn)主要基于數(shù)據(jù)并行的考量,通過(guò)劃分?jǐn)?shù)據(jù)到多個(gè)核或節(jié)點(diǎn)上計(jì)算,在負(fù)載均衡的條件下能取得較好的并行效果.但對(duì)于分治和遞歸類(lèi)的場(chǎng)景,當(dāng)數(shù)據(jù)依賴(lài)呈現(xiàn)復(fù)雜的網(wǎng)狀結(jié)構(gòu)時(shí),難以實(shí)現(xiàn)理想的負(fù)載均衡效果.

基于任務(wù)的并行編程模型按照功能將一個(gè)應(yīng)用劃分成多個(gè)能并行執(zhí)行的模塊,可以是細(xì)粒度的,也可以是粗粒度的.常見(jiàn)的任務(wù)并行模型有TBB,Cilk和SMPSS 等,使用這些模型時(shí)無(wú)需關(guān)心任務(wù)的具體分配和調(diào)度細(xì)節(jié).

Cilk 作為一個(gè)基于任務(wù)的并行編程模型,為共享內(nèi)存系統(tǒng)提供了一個(gè)高效的任務(wù)竊取調(diào)度器,適合用來(lái)為分治和遞歸類(lèi)算法實(shí)現(xiàn)高效的多核任務(wù)并行.對(duì)于并行編程,Cilk 為串行程序的并行化提供了向量化、鎖、視圖等機(jī)制[6],通過(guò)關(guān)鍵字來(lái)聲明操作,改造后的代碼具有語(yǔ)義化、可讀性強(qiáng)的特點(diǎn).

1 分治算法求解的基本原理

1.1 矩陣劃分

對(duì)任意n階實(shí)對(duì)稱(chēng)三對(duì)角矩陣T,基于秩1 修正從第k行作如下劃分：

其中,T′1,T′2是三對(duì)角子矩陣塊,ek代表第k個(gè)元素為1 的單位向量,元素β是T的第k個(gè)副對(duì)角元素.為使修正矩陣元素一致,對(duì)T1的末位元素和T2的首元素都減去了相同的元素.

分治過(guò)程一般采用折半劃分,對(duì)于劃分后的矩陣T′,若規(guī)模未達(dá)到指定閾值,則繼續(xù)劃分.

1.2 分治法求解過(guò)程

若劃分后的子矩陣規(guī)模達(dá)到設(shè)定閾值,調(diào)用QR或者其他算法得到子矩陣T’1,T’2的特征值分解：

令修正向量

則矩陣T的特征分解如下：

令D+βzzT=W,因?yàn)閐iag(Q1,Q2)為正交矩陣,所以T和W相似有相同的特征值,求T的特征分解T=QΛQT,可以通過(guò)先求出W的特征分解：W=PΛPT,再令Q=diag(Q1,Q2)P,即可求得T的特征分解.

對(duì)于W的特征分解,需要根據(jù)劃分后矩陣的特征值和原矩陣特征值之間的間隔性[7],通過(guò)求解對(duì)應(yīng)的secular 方程(3)：

得到相應(yīng)的的特征值λ.文獻(xiàn)[6]討論了方程(3)根的分布特性和解法,在特征區(qū)間內(nèi)使用牛頓法或拉格朗日法迭代逼近區(qū)間內(nèi)的特征值,線(xiàn)性時(shí)間內(nèi)即可收斂.解出所有的特征值后,再通過(guò)式(D-λI)-1x[1]求得每個(gè)特征值對(duì)應(yīng)的特征向量,回代到式(1)中得到T的特征分解.

整個(gè)求解過(guò)程在邏輯上可看作樹(shù)型結(jié)構(gòu),在葉子結(jié)點(diǎn)上進(jìn)行特征求解,在非葉子結(jié)點(diǎn)上進(jìn)行特征合并.

1.3 注意事項(xiàng)

文獻(xiàn)[5]中指出合并過(guò)程中如果D的主對(duì)角線(xiàn)上出現(xiàn)了相同的元素或者z向量中出現(xiàn)0 元素時(shí),在迭代前進(jìn)行deflation 操作可以避免特征區(qū)間內(nèi)迭代不收斂的情況.通過(guò)Givens 正交旋轉(zhuǎn)變換,部分特征值和特征向量就能原地得到.如果求得的特征值前后過(guò)于接近,按原方法求得的特征向量會(huì)丟失正交性,文獻(xiàn)[8,9]討論了特征值間隔過(guò)小時(shí)的特征向量的求法,避免了顯式正交過(guò)程.文獻(xiàn)[2]根據(jù)特征問(wèn)題的反問(wèn)題給出了特征向量的改進(jìn)求法,保證了計(jì)算結(jié)果的正交性.

2 分治算法的并行化

2.1 節(jié)點(diǎn)內(nèi)并行實(shí)現(xiàn)

單節(jié)點(diǎn)內(nèi)求解三對(duì)角矩陣塊的特征分解主要包含兩個(gè)過(guò)程,分別是求出最小劃分矩陣的特征分解以及逐步合并特征值和特征向量.由于每個(gè)子矩陣的特征計(jì)算以及同級(jí)子矩陣之間的合并過(guò)程是相互獨(dú)立的,下面結(jié)合Cilk 模型給出遞歸實(shí)現(xiàn)的分治算法在單節(jié)點(diǎn)多核環(huán)境下的并行化算法.

算法1.分治算法基于Cilk 的節(jié)點(diǎn)內(nèi)并行1) 判斷輸入矩陣的規(guī)模,如果矩陣規(guī)模小于等于設(shè)定 閾值,執(zhí)行2)；否則執(zhí)行3)；2) 調(diào)用特征求解算法進(jìn)行特征分解,返回結(jié)果；3) 劃分矩陣T 得子矩陣T1,T2,通過(guò)Cilk 執(zhí)行任務(wù)劃分,并行地對(duì)T1,T2 遞歸調(diào)用步驟1),求出子矩陣的特征值和特征向量后,進(jìn)行排序和迭代逼近等操作合并子矩陣的特征值和特征向量,得到T 的特征分解,返回結(jié)果.

Cilk 在程序遞歸執(zhí)行過(guò)程中不斷劃分新的任務(wù),產(chǎn)生Cilk 線(xiàn)程去處理,并把其分配到空閑的核上,和父線(xiàn)程并行執(zhí)行.一個(gè)Cilk 任務(wù)是在同步、生成新線(xiàn)程或返回結(jié)果之前執(zhí)行的最長(zhǎng)序列化指令集或代碼段[10].求解過(guò)程中,遞歸調(diào)用T1之前的過(guò)程可看作任務(wù)A,遞歸調(diào)用T2的過(guò)程可看作任務(wù)B,同步及合并過(guò)程可看作任務(wù)C.現(xiàn)假設(shè)程序只進(jìn)行了四次嵌套調(diào)用,則全局任務(wù)劃分后生成的任務(wù)流如圖1所示.

2.2 多節(jié)點(diǎn)混合并行實(shí)現(xiàn)

實(shí)現(xiàn)分治算法的多節(jié)點(diǎn)并行,需要先把原始矩陣按邏輯樹(shù)結(jié)構(gòu)劃分成小的子矩陣分發(fā)到葉子結(jié)點(diǎn)上進(jìn)行初步計(jì)算,然后通過(guò)MPI 通信匯總每棵子樹(shù)的計(jì)算結(jié)果,把整理后的結(jié)果重新發(fā)送到子結(jié)點(diǎn)進(jìn)行合并操作,重復(fù)此過(guò)程,直到返回到樹(shù)的根結(jié)點(diǎn).為了減少進(jìn)程間通信,充分利用Cilk 任務(wù)并行機(jī)制,需要對(duì)原矩陣進(jìn)行較粗粒度的劃分,以榨取單節(jié)點(diǎn)處理器的計(jì)算性能.

圖1 分治算法的任務(wù)并行流程

對(duì)于n階對(duì)稱(chēng)三對(duì)角矩陣T在N個(gè)進(jìn)程上的特征求解,假設(shè)N=2k,n=2j,算法初始時(shí)把原矩陣劃分為N個(gè)子矩陣,每個(gè)子矩陣階數(shù)為2j-k,下面給出分治算法在多節(jié)點(diǎn)上的混合并行實(shí)現(xiàn).

算法2.分治算法基于MPI/Cilk 的節(jié)點(diǎn)間并行首先,對(duì)于劃分到N 個(gè)進(jìn)程上的每個(gè)初始子矩陣,調(diào)用算法1 求出各部分特征值和特征向量；接著,進(jìn)行p 輪循環(huán)(p 對(duì)應(yīng)邏輯樹(shù)的高度).每輪循環(huán)先把進(jìn)程分組并確定主控進(jìn)程.每個(gè)MPI 進(jìn)程根據(jù)自己的進(jìn)程類(lèi)別按序選擇執(zhí)行下列步驟：1) 組內(nèi)進(jìn)程通過(guò)互補(bǔ)通信交換求得的局部特征向量矩陣,接著向主控進(jìn)程發(fā)送本進(jìn)程求得的局部特征值和拼接的修正向量；2) 主控進(jìn)程收集所有組內(nèi)進(jìn)程求得的局部特征值進(jìn)行排序,并保留排序后的位置索引數(shù)組,并根據(jù)此數(shù)組將修正向量排序,將排好序的特征值和修正向量劃分到組內(nèi)各進(jìn)程；3) 組內(nèi)進(jìn)程接收主控進(jìn)程發(fā)送過(guò)來(lái)的排好序的部分特征值和修正向量,求解secular 方程和特征向量；4) 各組內(nèi)進(jìn)程按列執(zhí)行分塊矩陣乘法,更新當(dāng)前輪循環(huán)的部分特征向量矩陣,接著執(zhí)行下一輪循環(huán).

2.3 算法分析

圖1最上方左面對(duì)應(yīng)遞歸的最外層,每個(gè)出度大于1 的結(jié)點(diǎn)都會(huì)派生出新的Cilk 線(xiàn)程,執(zhí)行新的任務(wù).程序執(zhí)行過(guò)程產(chǎn)生的所有的Cilk 線(xiàn)程數(shù)為圖1中節(jié)點(diǎn)數(shù),關(guān)鍵路徑長(zhǎng)度為圖1中虛線(xiàn)部分路徑長(zhǎng)度.程序的并行性為Cilk 線(xiàn)程數(shù)與關(guān)鍵路徑長(zhǎng)度之比,執(zhí)行時(shí)間大于等于所有任務(wù)平均到每個(gè)核上運(yùn)行的時(shí)間和關(guān)鍵路徑上任務(wù)運(yùn)行花費(fèi)的總時(shí)間.

算法2 中在對(duì)特征值進(jìn)行排序后,同時(shí)記錄排序后各位置元素的原來(lái)位置的索引,避免對(duì)特征向量重復(fù)排序,方便下一輪特征向量的計(jì)算.在第p輪循環(huán)中,從0 號(hào)進(jìn)程開(kāi)始每2p+1個(gè)進(jìn)程分為一組,每隔2p+1個(gè)進(jìn)程被選為選為主控線(xiàn)程.更新局部特征向量矩陣時(shí),為了減少數(shù)據(jù)通信開(kāi)銷(xiāo),避免主控進(jìn)程進(jìn)行統(tǒng)一收發(fā),使組內(nèi)進(jìn)程間只進(jìn)行必要的局部特征向量的互補(bǔ)交換,然后按列執(zhí)行分塊矩陣乘法,分別計(jì)算更新的特征向量矩陣的部分,通信量為O(n2/N).

算法1和算法2 中的合并過(guò)程涉及大量計(jì)算密集型操作,占據(jù)程序執(zhí)行的主要時(shí)間,可以使用Cilk 提供的數(shù)據(jù)并行機(jī)制加速.Cilk 通過(guò)數(shù)據(jù)劃分形成一個(gè)個(gè)可供調(diào)度的線(xiàn)程級(jí)任務(wù)并行執(zhí)行,通過(guò)基于貪心策略的任務(wù)竊取方式執(zhí)行調(diào)度：當(dāng)一個(gè)核完成自己的任務(wù)而其它核還有很多任務(wù)未完成,此時(shí)會(huì)將忙的核上的任務(wù)重新分配給空閑的核.

由于在迭代計(jì)算不同特征值的近似值時(shí),所執(zhí)行的迭代次數(shù)可能不同,均勻分配任務(wù)可能導(dǎo)致線(xiàn)程間的負(fù)載不平衡.Cilk 通過(guò)工作密取機(jī)制從其他線(xiàn)程獲取一部分工作量,能夠避免單核上出現(xiàn)任務(wù)過(guò)載和饑餓等待的問(wèn)題,同時(shí)能將密取的次數(shù)控制在最低水平,減少密取帶來(lái)的性能開(kāi)銷(xiāo).若單節(jié)點(diǎn)內(nèi)劃分的數(shù)據(jù)粒度適當(dāng),就能利用Cilk 機(jī)制充分發(fā)揮多核處理器的計(jì)算和調(diào)度能力,得到較好的負(fù)載均衡.

3 實(shí)驗(yàn)結(jié)果與分析

我們?cè)谥锌圃骸霸背闵系腸puII 計(jì)算隊(duì)列進(jìn)行了數(shù)值實(shí)驗(yàn),節(jié)點(diǎn)上搭載的是Intel E5-2680 V3 芯片,每塊芯片包含12 顆主頻為2.5 GHz 的計(jì)算核心.MPI程序使用Intel mpiicpc 編譯并鏈接到Cilk Plus 運(yùn)行時(shí).實(shí)驗(yàn)對(duì)象是30 000 階主對(duì)角元素是4,副對(duì)角元素是1 實(shí)對(duì)稱(chēng)三對(duì)角矩陣.矩陣劃分求解的規(guī)模閾值設(shè)為200 階,實(shí)驗(yàn)最終求出全部特征值和特征向量(通過(guò)環(huán)境變量設(shè)置每個(gè)節(jié)點(diǎn)使用2～12 個(gè)Cilk 線(xiàn)程進(jìn)行實(shí)驗(yàn)).

表1是MPI 模型與MPI/Cilk 混合模型在4～128個(gè)核、2～12 個(gè)Cilk 線(xiàn)程參與計(jì)算下所用時(shí)間匯總.從表中可以看出,相同計(jì)算條件下MPI/Cilk 混合并行計(jì)算時(shí)間要少于純MPI 方法,而隨著參與計(jì)算的Cilk 線(xiàn)程數(shù)的增加,計(jì)算用時(shí)也越來(lái)越少,這表明當(dāng)計(jì)算任務(wù)密集時(shí),計(jì)算效率隨著可供調(diào)度的線(xiàn)程的增加而增加.當(dāng)參與計(jì)算的Cilk 線(xiàn)程數(shù)超過(guò)10 時(shí),對(duì)于當(dāng)前規(guī)模的矩陣計(jì)算獲得的加速提升效果逐漸變得有限,這是因?yàn)镃ilk 是基于貪心策略執(zhí)行調(diào)度的,當(dāng)任務(wù)粒度較小時(shí)不會(huì)頻繁的執(zhí)行任務(wù)竊取,從而減小系統(tǒng)開(kāi)銷(xiāo).而在較少節(jié)點(diǎn)參與計(jì)算時(shí),獲得的加速效果較為明顯,這說(shuō)明對(duì)數(shù)據(jù)進(jìn)行粗粒度劃分時(shí)能取得較好性能.

圖2是MPI 模型和MPI/Cilk 模型在不同參數(shù)下的加速比變化趨勢(shì)曲線(xiàn)(已測(cè)出在單節(jié)點(diǎn)單線(xiàn)程環(huán)境下的平均串行時(shí)間為880 秒左右).可以看出,隨著參與計(jì)算的核數(shù)和Cilk 線(xiàn)程數(shù)增加,三條曲線(xiàn)的上升趨勢(shì)都逐漸變得平緩,出現(xiàn)這種現(xiàn)象的原因主要為計(jì)算任務(wù)粒度變細(xì)和進(jìn)程間通信開(kāi)銷(xiāo)增加導(dǎo)致的.MPI/Cilk方法的曲線(xiàn)在一開(kāi)始上升較快,隨著進(jìn)程數(shù)的增加,計(jì)算任務(wù)粒度逐漸減小,加速效果也隨之下降,但仍比純MPI 方法變緩得更慢.這表明MPI/Cilk 方法的并行效果要優(yōu)于純MPI 方法,擁有更好的可擴(kuò)展性.隨著可供調(diào)度的核數(shù)增多,使用更多的Cilk 線(xiàn)程數(shù)參與計(jì)算能獲得更好的計(jì)算性能,Cilk 動(dòng)態(tài)調(diào)度的優(yōu)勢(shì)也漸漸體現(xiàn)了出來(lái).

表2給出了MPI/Cilk 模型(使用8 個(gè)Cilk 線(xiàn)程)同ELPA、ScaLAPACK 軟件包實(shí)現(xiàn)的分治算法求解30000 階矩陣特征的計(jì)算時(shí)間比較.ELPA和ScaLAPACK 是求解線(xiàn)性代數(shù)問(wèn)題領(lǐng)域中應(yīng)用廣泛的并行軟件包.從表中數(shù)據(jù)可以得出,當(dāng)參與計(jì)算的核數(shù)較少時(shí),MPI/Cilk 方法與ELPA、Scalapack 中分治算法計(jì)算用時(shí)接近,隨著參與計(jì)算的核數(shù)增多,MPI/Cilk方法、ELPA 方法與ScaLAPACK 方法的計(jì)算性能逐漸拉開(kāi),但本文的MPI/Cilk 方法與ELPA 方法的擴(kuò)展性仍然存在著一定差距.這是由于ELPA 按照塊級(jí)循環(huán)分布矩陣,實(shí)現(xiàn)了對(duì)計(jì)算流程的更細(xì)粒度的控制.

表1 MPI 模型和MPI/Cilk 模型運(yùn)行時(shí)間

圖2 MPI 模型和MPI/Cilk 模型的加速比曲線(xiàn)

表2 MPI/Cilk 方法與ELPA、ScaLAPACK 求解器時(shí)間對(duì)比

4 結(jié)論與展望

本文針對(duì)對(duì)稱(chēng)三對(duì)角矩陣特征求解的分治算法,提出了一種改進(jìn)的基于MPI/Cilk 的混合并行實(shí)現(xiàn).在節(jié)點(diǎn)間,進(jìn)行粗粒度計(jì)算任務(wù)的劃分,更新局部特征矩陣時(shí)通過(guò)弱化主控進(jìn)程的中心作用,使組內(nèi)進(jìn)程之間只進(jìn)行必要的互補(bǔ)數(shù)據(jù)交換,避免了統(tǒng)一收發(fā)流程,減少了進(jìn)程間的通信開(kāi)銷(xiāo)；在節(jié)點(diǎn)內(nèi),利用Cilk 的任務(wù)并行機(jī)制提高了CPU 的利用率和特征求解過(guò)程中的負(fù)載均衡性.實(shí)驗(yàn)結(jié)果體現(xiàn)了該算法的性能.本文實(shí)現(xiàn)的分治算法與最新的ELPA 軟件包仍然存在一定的性能差距,還有可以提升的空間.此外,Cilk 線(xiàn)程啟動(dòng)數(shù)和數(shù)據(jù)劃分粒度對(duì)計(jì)算性能的影響以及同openMP 等其它任務(wù)并行模型的效果對(duì)比等等是值得進(jìn)行下一步研究的方向.