基于分治法求解對稱三對角矩陣特征問題的MPI/Cilk混合并行算法

朱京喬， 趙永華

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

0 引言

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

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

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

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

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

1.1 三對角矩陣的劃分

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

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

T=T′+βzzT，

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

1.2 分治法求解過程

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

(1)

此時T可以表示為

(2)

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

(3)

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

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

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

1.3 求解注意事項

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

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

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

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

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

算法1分治算法基于Cilk的節(jié)點內(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，并行對T1、T2遞歸調(diào)用算法1，求出子矩陣的特征分解后，進(jìn)行排序和迭代逼近等操作，合并子矩陣的特征值和特征向量，得到T的特征分解，返回結(jié)果。

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

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

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

2.2 算法分析

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

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

Parallelism=W/S。

(4)

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

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

(5)

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

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

3.1 混合并行算法

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

對于n階對稱三對角矩陣T在N個進(jìn)程上的特征求解，假設(shè)N=2q，n=2p，算法初始時把原矩陣劃分為N個子矩陣，每個子矩陣階數(shù)為2p-q，設(shè)為2k，下面給出分治算法的在多進(jìn)程間的混合并行實現(xiàn)。

算法2分治算法基于MPI/Cilk的節(jié)點間并行。

第一步，對于劃分到N個進(jìn)程上的每個初始子矩陣，調(diào)用算法1求出各部分特征值和特征向量。

第二步，進(jìn)行p次循環(huán)(p對應(yīng)圖1中樹的深度)，每次循環(huán)先把進(jìn)程分組并確定主控進(jìn)程來進(jìn)行結(jié)果收發(fā)。對每個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ā)送上一輪計算得到的特征向量矩陣。

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

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

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

3.2 算法分析

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

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

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

4 實驗結(jié)果與分析

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

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

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

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

表1 MPI模型和MPI/Cilk模型程序的時間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模型程序的時間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é)

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