朱京喬， 趙永華

(1. 中國(guó)科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心高性能部 北京 100190； 2. 中國(guó)科學(xué)院大學(xué) 計(jì)算與控制學(xué)院 北京 100049)

0 引言

科學(xué)計(jì)算和工程問題中經(jīng)常涉及稠密對(duì)稱矩陣的特征求解問題，通常的處理方法是先通過正交變換把原矩陣化為三對(duì)角矩陣，如Householder或Givens方法，再求解三對(duì)角陣的特征值和特征向量。三對(duì)角陣的特征求解算法有很多種，常用的有QR法、二分法、分而治之方法和多近似健壯表示(multiple relatively robust representations, MRRR)方法等等，本文工作基于分治法展開。

對(duì)稱三對(duì)角陣特征求解的分治算法，最早由文獻(xiàn)[1]提出,文獻(xiàn)[2]作出了改進(jìn)，提高了該算法的求解精度和特征向量的正交性。分治算法采用分而治之思想，將原矩陣劃分為若干小的子矩陣，先求出子矩陣的特征分解，接著通過修正計(jì)算，逐級(jí)合并子矩陣的結(jié)果，回代得到原矩陣的特征分解。分治法具有很強(qiáng)的靈活性，對(duì)于大規(guī)模三對(duì)角矩陣特征問題的并行求解是一種比較理想的算法。

現(xiàn)有很多模型都實(shí)現(xiàn)了針對(duì)分治算法的并行化，如基于分布式存儲(chǔ)的消息傳遞接口(message passing interface,MPI)進(jìn)程級(jí)并行[3]，基于共享內(nèi)存的openMP線程級(jí)并行，還有基于對(duì)稱多處理結(jié)構(gòu)(symmetrical multi procession,SMP)的MPI/openMP混合并行模型[4]等。以openMP為代表的線程級(jí)并行模型主要基于數(shù)據(jù)并行的考量，通過把數(shù)據(jù)劃分到多個(gè)核上實(shí)現(xiàn)并行，該模型在計(jì)算負(fù)載均衡的條件下能取得較好的效果。但對(duì)于分治和遞歸類的場(chǎng)景，當(dāng)數(shù)據(jù)依賴呈復(fù)雜網(wǎng)狀結(jié)構(gòu)時(shí)，則不能很好地實(shí)現(xiàn)負(fù)載均衡。

Cilk是一個(gè)基于任務(wù)的并行編程模型，為串行程序的并行化提供了向量化、視圖等機(jī)制(http:∥supertech.csail.mit.edu/cilk/manual-5.4.6.pdf)，通過關(guān)鍵字來(lái)聲明操作，改造后的代碼具有語(yǔ)義化、可讀性強(qiáng)的特點(diǎn)。Cilk為共享內(nèi)存系統(tǒng)提供了一個(gè)高效的任務(wù)竊取調(diào)度器，適合用來(lái)實(shí)現(xiàn)高效的多核任務(wù)并行。

1 三對(duì)角矩陣特征求解分治算法的基本原理

1.1 三對(duì)角矩陣的劃分

對(duì)任意n階實(shí)對(duì)稱三對(duì)角矩陣T，基于秩1修正作如下劃分。

其中：T1、T2是三對(duì)角子矩陣塊，β元素所在矩陣構(gòu)成原矩陣的修正矩陣。為使修正矩陣元素一致，T1的末位元素和T2的首元素都減去相同的β。 經(jīng)過劃分原矩陣變?yōu)?/p>

T=T′+βzzT，

其中Z=(0,…,0,1,1,0,…,0)T。 1的索引代表對(duì)原矩陣進(jìn)行劃分的位置。 對(duì)于劃分后的矩陣T′，若規(guī)模未達(dá)到指定閾值，則繼續(xù)劃分。

1.2 分治法求解過程

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

(1)

此時(shí)T可以表示為

(2)

對(duì)于D+βxxT的特征分解，需要根據(jù)劃分后矩陣的特征值和原矩陣特征值的間隔性[5]，通過求解對(duì)應(yīng)的secular方程

(3)

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

圖1 分治算法結(jié)構(gòu)Figure 1 The structure of divide and conquer method

整個(gè)求解過程可看作邏輯上的樹型結(jié)構(gòu)，如圖1所示。在葉子結(jié)點(diǎn)進(jìn)行特征求解，在非葉子結(jié)點(diǎn)上進(jìn)行特征合并。

1.3 求解注意事項(xiàng)

降階：合并過程中如果D的主對(duì)角線上出現(xiàn)了相同的元素或者X向量中出現(xiàn)0元素時(shí)，文獻(xiàn)[6]中提出要先進(jìn)行deflation操作，通過Givens正交旋轉(zhuǎn)變換進(jìn)行迭代逼近前的預(yù)處理，部分特征值和特征向量就能原地得到，從而避免了特征區(qū)間內(nèi)迭代不收斂的情況。

正交性：如果特征值前后間隔過小，使用原公式求解得到的特征向量會(huì)逐漸丟失正交性。 文獻(xiàn)[2]給出了特征向量的改進(jìn)求法，引入了新的計(jì)算方式來(lái)修正計(jì)算過程中結(jié)果的正交性。

2 基于Cilk的節(jié)點(diǎn)內(nèi)并行

2.1 節(jié)點(diǎn)內(nèi)并行算法

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

算法1分治算法基于Cilk的節(jié)點(diǎn)內(nèi)并行。

dc_solver(T):if (size T==threshold): exec_solve(T); return result;else: cilk_spawn dc_solver(T1); cilk_spawn dc_solver(T2); cilk_sync; exec_merge(); return result;end

先判斷輸入矩陣的規(guī)模，如果矩陣規(guī)模小于等于指定閾值，調(diào)用QR或其他算法進(jìn)行特征求解，返回結(jié)果；否則執(zhí)行Cilk任務(wù)劃分，劃分矩陣T得子矩陣T1、T2，并行對(duì)T1、T2遞歸調(diào)用算法1，求出子矩陣的特征分解后，進(jìn)行排序和迭代逼近等操作，合并子矩陣的特征值和特征向量，得到T的特征分解，返回結(jié)果。

Cilk在程序遞歸執(zhí)行過程中不斷劃分新的任務(wù)，產(chǎn)生Cilk線程去處理，并把其分配到空閑的核上，和父線程并行執(zhí)行。 Cilk通過任務(wù)竊取(Work-Stealing)的方式來(lái)執(zhí)行調(diào)度：當(dāng)一個(gè)核完成自己的任務(wù)而其他核還有很多任務(wù)未完成，此時(shí)會(huì)將忙的核上的任務(wù)重新分配給空閑的核。

算法1計(jì)算過程中涉及矩陣乘法、排序、迭代逼近等操作，可以利用Cilk提供的數(shù)據(jù)并行機(jī)制加速，Cilk通過數(shù)據(jù)劃分形成一個(gè)個(gè)可供調(diào)度的線程級(jí)任務(wù)并行執(zhí)行。 空閑線程通過工作密取從其他線程獲取一部分工作量，Cilk利用貪心策略來(lái)進(jìn)行任務(wù)調(diào)度分配，能夠避免單核上出現(xiàn)任務(wù)過載和饑餓等待的問題，同時(shí)能將密取的次數(shù)控制在最低水平，減少密取帶來(lái)的性能開銷。

圖2 節(jié)點(diǎn)內(nèi)任務(wù)并行流程Figure 2 Task flow of divide and conquer method inside nodes

2.2 算法分析

一個(gè)Cilk線程可看作在同步、生成新線程或返回結(jié)果之前所能執(zhí)行的最長(zhǎng)序列化指令集(http:∥www.cilkplus.org/cilk-documentation-full)。 根據(jù)此定義，求解過程中，遞歸調(diào)用T1之前的過程可以看作一個(gè)Cilk線程，記為任務(wù)A， 遞歸調(diào)用T2的過程可以看作任務(wù)B，同步及合并過程可以看作任務(wù)C。 現(xiàn)假設(shè)程序只進(jìn)行了4次嵌套調(diào)用，則全局任務(wù)劃分后生成的任務(wù)流程如圖2所示。

程序從有向無(wú)環(huán)圖的最上面開始執(zhí)行，對(duì)應(yīng)遞歸的最外層。每一個(gè)出度大于1的結(jié)點(diǎn)都會(huì)派生出新的Cilk線程，執(zhí)行新的任務(wù)。設(shè)程序執(zhí)行過程產(chǎn)生的所有的Cilk線程數(shù)為W(圖2中節(jié)點(diǎn)數(shù))，關(guān)鍵路徑長(zhǎng)度為S(圖2中虛線部分路徑長(zhǎng)度)。 則程序的并行性Parallelism可表示為

Parallelism=W/S。

(4)

程序的執(zhí)行時(shí)間Time滿足

Time>max(T(W/P),T(S)),

(5)

其中:P是處理器核數(shù)；T(W/P)表示所有任務(wù)平均到每個(gè)核上運(yùn)行的時(shí)間；T(S)表示關(guān)鍵路徑上的任務(wù)運(yùn)行花費(fèi)的總時(shí)間。

3 基于MPI/Cilk的多節(jié)點(diǎn)混合并行

3.1 混合并行算法

實(shí)現(xiàn)分治算法的MPI多節(jié)點(diǎn)并行，需要把原始矩陣按圖1結(jié)構(gòu)劃分成小的子矩陣分發(fā)到葉子結(jié)點(diǎn)上進(jìn)行初步計(jì)算，然后通過MPI通信，每棵子樹的父結(jié)點(diǎn)收集子結(jié)點(diǎn)的計(jì)算結(jié)果后進(jìn)行匯總，把整理后的結(jié)果重新發(fā)送到子結(jié)點(diǎn)進(jìn)行合并操作，重復(fù)此過程，直到返回到圖1中的根結(jié)點(diǎn)。 為了減少進(jìn)程間通信，充分利用Cilk任務(wù)并行機(jī)制，需要對(duì)原矩陣進(jìn)行較粗粒度的劃分，以榨取單節(jié)點(diǎn)處理器的計(jì)算性能。 并且在消息傳遞時(shí)，只使用一個(gè)線程作為主控線程，以減少通信開銷。

對(duì)于n階對(duì)稱三對(duì)角矩陣T在N個(gè)進(jìn)程上的特征求解，假設(shè)N=2q，n=2p，算法初始時(shí)把原矩陣劃分為N個(gè)子矩陣，每個(gè)子矩陣階數(shù)為2p-q，設(shè)為2k，下面給出分治算法的在多進(jì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)圖1中樹的深度)，每次循環(huán)先把進(jìn)程分組并確定主控進(jìn)程來(lái)進(jìn)行結(jié)果收發(fā)。對(duì)每個(gè)MPI進(jìn)程，需要根據(jù)自己的進(jìn)程類別，按序選擇執(zhí)行下列步驟：

1) 組內(nèi)進(jìn)程向主控進(jìn)程發(fā)送本進(jìn)程求得的所有局部特征值，如果本進(jìn)程是上一輪循環(huán)的主控進(jìn)程，則還要向本輪主控進(jìn)程發(fā)送上一輪計(jì)算得到的特征向量矩陣。

2) 主控進(jìn)程收集屬于同一組進(jìn)程中所有求得的局部特征值和特征向量，對(duì)特征值進(jìn)行排序，并保留排序后的位置索引數(shù)組。 按照1.2節(jié)的過程拼接Z1、Z2成一個(gè)新向量X，根據(jù)位置索引將X排序,將排好序的特征值均勻地劃發(fā)到組內(nèi)各進(jìn)程，同時(shí)也把向量X傳到組內(nèi)各進(jìn)程。

3) 組內(nèi)進(jìn)程接收主控進(jìn)程發(fā)送過來(lái)的排好序的部分特征值和向量X，求解secular方程和特征向量，并發(fā)送結(jié)果給主控進(jìn)程。

4) 主控進(jìn)程更新當(dāng)前循環(huán)的特征向量矩陣，接著執(zhí)行下一輪循環(huán)。

3.2 算法分析

在對(duì)特征值進(jìn)行排序后，同時(shí)記錄排序后各位置元素的原來(lái)位置的索引，避免對(duì)特征向量重復(fù)排序，方便下一輪特征向量的計(jì)算。 算法2中步驟3)和4)中涉及大量計(jì)算密集型操作，占據(jù)程序執(zhí)行的主要時(shí)間，同樣可以通過Cilk數(shù)據(jù)并行機(jī)制劃分多個(gè)線程級(jí)子任務(wù)并行執(zhí)行。

在第p輪循環(huán)中，從0號(hào)進(jìn)程開始,每2p+1個(gè)進(jìn)程分為一組，每隔2p+1個(gè)進(jìn)程被選為主控線程，每個(gè)小組內(nèi)共進(jìn)行3(2p+1-1)次通信，每輪循環(huán)的通信次數(shù)是3N(2p+1-1)/2p+1。

若單節(jié)點(diǎn)內(nèi)劃分的數(shù)據(jù)粒度適當(dāng)粗時(shí)，就能利用Cilk機(jī)制充分發(fā)揮多核處理器的計(jì)算和調(diào)度能力，得到較好的計(jì)算性能，同時(shí)減少進(jìn)程間通信次數(shù)，降低通信開銷。

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

我們?cè)谥锌圃骸霸奔?jí)超算上針對(duì)該混合模型進(jìn)行了數(shù)值實(shí)驗(yàn)，實(shí)驗(yàn)運(yùn)行在CPU II計(jì)算隊(duì)列上，節(jié)點(diǎn)上搭載的是Intel E5-2680 V3芯片，每塊芯片包含12個(gè)主頻為2.5 GHz的計(jì)算核心。 MPI程序使用Intel mpiicpc編譯器編譯，編譯時(shí)鏈接Cilk Plus運(yùn)行，計(jì)算過程中使用了部分MKL庫(kù)中的函數(shù)。 實(shí)驗(yàn)對(duì)象是30 000階實(shí)對(duì)稱三對(duì)角矩陣，矩陣劃分求解的規(guī)模閾值設(shè)為200階，最終求出全部特征值和特征向量。 實(shí)驗(yàn)通過環(huán)境變量設(shè)置使用4和8個(gè)線程進(jìn)行實(shí)驗(yàn)，為方便對(duì)比，已測(cè)出實(shí)驗(yàn)在單節(jié)點(diǎn)單線程環(huán)境下平均串行時(shí)間為882.53 s。

表1是純MPI方法和MPI/Cilk方法在4～64個(gè)核參與計(jì)算下所用時(shí)間比較。 從表1可以看出，對(duì)于同一矩陣，使用相同數(shù)目的核參與計(jì)算時(shí)，MPI/Cilk混合方法所用時(shí)間要少于純MPI方法，性能要優(yōu)于純MPI方法，而8線程的MPI/Cilk方法性能要優(yōu)于4線程的，這表明8線程的混合模型有更好的負(fù)載均衡。 在較少節(jié)點(diǎn)參與計(jì)算時(shí)，獲得加速效果較為明顯，這說明對(duì)數(shù)據(jù)進(jìn)行粗粒度劃分時(shí)能取得較好性能。

表2是MPI模型和MPI/Cilk模型程序的加速比對(duì)比。 表2數(shù)據(jù)表明，MPI/Cilk方法具有比純MPI方法更高的加速比，而8線程MPI/Cilk方法的加速比開始時(shí)和4線程接近，隨著可供調(diào)度的核數(shù)增多，加速效果逐漸超過了4線程方法。 MPI/Cilk方法的加速比在一開始上升較快，隨著進(jìn)程數(shù)的增加，加速效果隨之下降，但仍比純MPI方法變緩更慢，擁有更好的擴(kuò)展性。 對(duì)于給定規(guī)模的矩陣，當(dāng)參與計(jì)算的核數(shù)超過一定數(shù)量時(shí)，計(jì)算效果的提升越來(lái)越有限，出現(xiàn)這種現(xiàn)象主要是因?yàn)橛?jì)算任務(wù)粒度變細(xì)和進(jìn)程間通信開銷增加導(dǎo)致的。

表3和表4分別給出了Cilk和openMP這兩種模型的計(jì)算時(shí)間和加速比的比較，其中openMP 模型采用數(shù)據(jù)并行方式，通過指導(dǎo)語(yǔ)句對(duì)涉及計(jì)算的部分進(jìn)行了并行化。從結(jié)果可以看出，在參與計(jì)算的核數(shù)較少時(shí)，openMP模型的效果要優(yōu)于Cilk模型，原因是計(jì)算粒度較大時(shí)，數(shù)據(jù)并行占主導(dǎo)地位。 而隨著參與計(jì)算的核數(shù)增多，Cilk模型的效果逐漸超過openMP模型，此時(shí)可供調(diào)度的核數(shù)增多，Cilk動(dòng)態(tài)調(diào)度的優(yōu)勢(shì)逐漸體現(xiàn)。

表1 MPI模型和MPI/Cilk模型程序的時(shí)間Table 1 Time of pure MPI method and MPI/Cilk method

表2 MPI模型和MPI/Cilk模型程序的加速比Table 2 Speedup of pure MPI method and MPI/Cilk

表3 MPI/openMP模型和MPI/Cilk模型程序的時(shí)間Table 3 Time of MPI/openMP method and MPI/Cilk method

表4 MPI/openMP模型和MPI/Cilk模型程序的加速比Table 4 Speedup of MPI/openMP method and MPI/Cilk method

5 總結(jié)

本文針對(duì)實(shí)對(duì)稱三對(duì)角矩陣特征求解的分治算法，提出了一種基于MPI/Cilk的混合并行實(shí)現(xiàn)算法。在節(jié)點(diǎn)間，通過粗粒度計(jì)算任務(wù)的劃分，在相同的求解效率下減少了所需要的計(jì)算進(jìn)程和進(jìn)程間的通信次數(shù)。在節(jié)點(diǎn)內(nèi)，利用Cilk的任務(wù)并行機(jī)制提高CPU的利用率，改善了負(fù)載均衡，提高了并行度。 實(shí)驗(yàn)結(jié)果體現(xiàn)了該算法的性能。 該算法還可進(jìn)一步研究，如運(yùn)用數(shù)據(jù)并行機(jī)制加速特征合并過程，研究Cilk線程啟動(dòng)數(shù)和數(shù)據(jù)劃分粒度對(duì)性能的影響以及同openMP的任務(wù)并行機(jī)制對(duì)比等。