一種基于優化存儲格式的DLB_GaBP算法

陳振武，黃 婧，蘭添才，鄭漢垣

(1.龍巖學院 數學與信息工程學院，福建 龍巖 364012；2.龍巖學院 大數據挖掘與應用福建省重點實驗室，福建 龍巖 364012;3.龍巖學院 傳播與設計學院，福建 龍巖 364012)

1 概 述

現代多核處理機或多處理機系統速度能力已實現了每秒千萬億次浮點運算的量級，因此，要提高大規模的數值計算的速度與效率，最終均歸結于操作系統的任務流調度算法、存儲空間存儲格式的合理定義與分隔。

對于多核處理機環境中的任務流調度算法可分為全局隊列調度算法和局部隊列調度算法兩種。局部隊列調度算法是為每個處理機核心維護一個局部任務等待隊列，處理機的利用率較低，在大規模并行計算中一般不予使用；全局隊列調度算法從處理機的任務均衡性、存儲格式的合理性方面綜合考慮充分發揮系統中各處理機核的綜合性能，當某一處理機空閑時，系統立即從全局任務等待隊列中選取就緒任務分配給該處理機執行，從而提高處理機的利用率，達到提高系統的運算速度與效率的目的。

吉林大學的耿曉中提出了一種動態負載均衡模型，并將影響多核處理器負載均衡的因素分為5類：多核系統的負載均衡環境、用戶提交的任務屬性、系統的負載評價、系統所采用的調度策略以及系統的調度評價指標[1]，為多核動態調度因素提供依據。

文獻[2]針對動態調度算法的基本原理和執行過程，結合負載均衡策略，提出了一種基于異構多核處理器的STDS動態任務調度算法，并通過選擇調度時間、負載均衡和任務等待時間三個方面證明了STDS算法在保證速度性能的同時有較理想的內核負載均衡效果。

賈燕成等提出了負載動態均衡規則[3]是在任務的調度過程中根據處理器的狀態，將較大的任務分解并動態映射到其他節點上，并與該節點的原始任務組進行組合形成后繼任務循環執行，使各個節點的任務達到動態平衡，可以避免節點的空閑等待，也縮短了各個節點通信開銷，最大化提高了并行效率。

近年來，隨著互聯網+、物聯網、大數據云計算應用領域的不斷開展，計算機系統的數據運算量也在不斷增加。如電力系統運行控制的潮流計算問題[4]、環境監控、實時評價問題、地震實時監測等領域，均涉及到大規模數值計算問題。而對這些大規模數值問題的求解，大多都是通過相應的數學建模，然后將問題的求解轉換成為相應建立起來的稀疏線性方程組Ax=b的求解[5-7]。

迭代法是在稀疏線性方程組求解中常用的方法，尤其是近年來，基于Krylov子空間方法[8-9]的迭代法得到了廣泛應用，Ori-Shental等人基于遞歸更新的概率推理算法提出了GaBP迭代算法[10]。

文獻[11-13]針對GaBP迭代算法具有計算復雜度低、并行度高的特性，在GaBP的迭代加速優化方法的基礎上給出了對應的多種GaBP迭代加速優化算法，并從動態松馳因子的GaBP算法和MannGaBP迭代加速優化算法的實驗，構造多核處理機動態因子。

文中針對多核處理機環境中的大規模稀疏矩陣線性方程組的存儲格式、處理機間的任務均衡性、并行算法的優化等進行研究，在優化稀疏矩陣存儲格式的基礎上，依據文獻[3]提出的負載動態均衡規則，對多核處理機環境的并行GaBP算法進行優化，設計實現具有動態負載均衡特性的多核并行GaBP算法(DLB_GaBP算法)。經過實驗驗證，該算法在對大規模稀疏矩陣線性方程的求解過程中，具有更高的加速效果、更快的計算速度、更好的環境適應能力。

2 高斯置信傳播(GaBP)算法

2.1 算法思想

GaBP算法主要包括：初始化和迭代(包括消息累加、消息傳播與更新、求解向量三個過程[13])兩個步驟。具體可描述如下：

步驟一：初始化。

(2)初始化。

步驟二：迭代計算。

(1)消息累加。

//N(i)為第i行節點編號的集合(不含i);

(2)消息傳播與更新。

(3)求解向量。

2.2 GaBP算法的改進

3 存儲格式定義

文獻[16]在充分分析了ELL、CSR、COO、ELL+COO混合的HYB(Hybrid)等多種存儲格式的優缺點的基礎上，提出了ELL+CSR混合的HEC存儲格式，并進行了相應的實驗對比驗證，也發現了HEC存儲格式比CSR格式多存儲了2組分別用以存儲格式本身所需設定的值和列數據，對讀取帶來一定的開銷。

為了充分利用GaBP算法的特殊性，文中定義動態負載均衡多核并行GaBP的數據結構，在數據存儲格式及數據讀取方法上采用列壓縮存儲格式，格式本身分為兩種形式：改進的列壓縮存儲格式(modified compressed spare column，MCSC)或優化的列壓縮存儲格式(optimized compressed spare column，OCSC)。

3.1 矩陣數據參數設定

設：n階稀疏矩陣矩陣Α，有非零元素nz個；

另設：存儲單元數組ADN、ALN、AUN、NL、NJL、NU、NJU。

其中，ADN用于存放按順序排列的n個對角矩陣元素值ai；

ALN用于按列存儲Α左下三角(nz-n)/2個非零元素值aij(i>j)，NL用于存儲行號，NJL用于存儲非零元素對應列起始行號(NJL=1～n)；

AUN用于按列存儲Α右上三角(nnz-n)/2個非零元素值aij(i

3.2 存儲格式定義

由于文中利用UFgett中的稀疏矩陣進行實驗驗證[17]，Α按行號、列號、元素值三個元素組進行存儲，可建立如下變換存儲格式的左端矩陣數據結構：

條件一：若Α不對稱

則：按ADN，ALN，NL，NJL，AUN，NU，NJU的順序的MCSC格式存儲。

條件二：若Α對稱

則：按AD，AL，NBL，NBJL順序的OCSC格式存儲。

條件三：若Α是對角稠密矩陣，則非零元素的存取與查找，采用OCSC格式存儲結合順序查找算法可以非常直接地進行。

條件四：若Α是非對角稠密矩陣，非零元素，一般可采用二分查找方法進行。依據GaBP算法稀疏性的特性，形狀對稱矩陣P、μ、Α(其中Α是左端數值對稱矩陣)，利用OCSC存儲它們，可以將NL、NJL與ALN、ADN的存儲格式進行一樣的定義。但由于矩陣μij是形狀對稱的，這需將其上三角元素形成矩陣存儲，可以借助MCSC格式進行存儲，因此，矩陣μij可對應為：μAD，μAL，μNL，μNJL，μAU，μNU，μNJU。

3.3 存儲格式的比較

圖1給出了稠密矩陣存儲格式與MCSC、OCSC存儲格式在存儲空間的占比情況。

圖1 存儲空間比值

從圖1可以明顯看出，MCSC和OCSC與稠密存儲空間的比值隨著矩陣階數、矩陣的稀疏性和非零元素數量規模的增大而增大，說明MCSC和OCSC格式所需的存儲空間遠少于稠密存儲空間，且OCSC格式所需的存儲空間比ACSC格式的存儲空間更少。

4 基于列壓縮存儲格式的動態負載均衡

4.1 負載均衡存儲格式劃分

對稀疏線性方程組求解的并行算法首先需對矩陣進行數據劃分，劃分的方式可以有列劃分、行劃分。文中采用按列數相同的列塊順序進行存儲格式劃分，確保各線程進行計算時的列塊中的列數相等。

由于Α的每一列中所包含非零元個數不一定相同，對應分配到各線程中的列塊中包含的非零元數往往也不一定相同，很明顯，各線程對列的計算量與該列中的包含非零元素的個數成正比，這很容易引發線程負載的不均衡，線程負載不均衡勢必進一步引發部分處理機的資源空閑浪費。因此，列劃分應盡量使分配給各個線程計算列塊中所包含的非零元素的個數差異性不大，以確保形成均衡的負載。

4.2 動態列矩陣塊的劃分算法

系統線程調度、運行過程是動態的，為確保建立動態的負載均衡，而現代操作系統針對線程的調度都是建立在具體算法的基礎上，下面是列矩陣塊的動態劃分算法步驟。

1:partition(seg[],NZ[],n,p) //定義參數，seg[]是存放矩陣列塊劃分方式的數組，NZ[]是存放非零元素的數組，p為線程數；

2:k=0,q=0,l=0;seg[0]=0 //初始化；

3:for allj∈0,1,…,ndo

4:l=l+NZ[j] //NZ[j]為第j列非零元素個數；

5:end for

6:for allj∈0,1,…,ndo //迭代劃分列塊；

7:k=k+NZ[j]

8:ifk>[q+1]*l/pthen

9:seg[q+1]=j;q=q+1

10:end if

11:seg[p]=n+1

12:end for

由于p為線程數，其結果是對應于q=0,1,…,p-1時，將稀疏矩陣的第segq到segq+1-1列分配到線程號q中，這就產生兩種運行狀態需要進行處理。

節點狀態處理：每一輪迭代結束均判定節點狀態，如果處于消息動態均衡，則不作計算處理，以減少列塊中待計算節點的數量。

動態均衡處理：為防止因減少計算節點引起新的負載不均衡，新一輪的迭代需重新對列塊進行劃分，并重新分配給相應線程，讓每一線程重新達到相同的計算量。以此類推，達到動態負載均衡。

4.3 DLB_GaBP算法

DLB_GaBP算法的實現實質上就是將GaBP算法中的消息累加、消息傳播與更新、求解向量分別加入到OpenMP的并行制導控制段中，構成三段多核并行控制段，即基于OpenMP的多核并行GaBP算法。該算法可表述如下：

Initialization: //初始化

1:init():

Iteration

2:repeat

5:if(iters mod D_LOOP==0)then//列塊劃分迭代步數條件判斷

6:partition(seg,NZ,n,Nthreads)

7:end if

8:#pragma omp parallel private(tid,p,j)threads(NTHREADS)

9:{tid=omp_get_thread_num()

10: for(i=seg[tid];j

&&Ti==0;i++)do

11:if(Ti==0)then

13:for allj∈N(i)jdo

15:end for

16:end if

17:end for

18:}

19:#pragma omp parallel private(tid,p,j)threads(NTHREADS)

20:{tid=omp_get_thread_num()

21: for(i=seg[tid];j

23:xi=μi/Pi//xi為收斂精度

24:if(Ti==0)and(xiconverged)then //xi滿足設定的收斂精度，該節點的Ti值設定為1

25:Ti=1;NZj=0

26:end if

27:end for

28:}

29:until convergence:allxiconverged

30:output:x*=[xi]

為能夠盡量好地提高計算效率，文中通過迭代步來控制上述的有效組合，具體是通過將迭代的次數限制一定數量范圍內，超過即重新進行列塊劃分。即在Partition()函數前面添加了迭代步控制，每隔D_LOOP迭代步重新進行列矩陣塊的劃分。D_LOOP的取值一般是依據劃分列塊中的列數，當然也可以是定量或變量，這只要在算法中進行相應的定義即可。

當收斂精度xi滿足設定的精度要求時，返回該節點Ti=1，NZ[j]=0，下一輪迭代計算時跳過該節點，同時，負載均衡處理中也不統計該節點的鄰居節點數，這可以大大減少節點的計算負載量，但是需要進行列塊的重新劃分以實現新的負載均衡。

新一輪迭代計算前進行新的列塊均衡劃分，總體上減少了計算量，也保持了各線程的計算量相對均衡，進而解決了因線程同步造成的大量線程的長時間等待問題，同時加速了迭代計算。

5 數值試驗

5.1 實驗環境

實驗在高性能集群系統中進行，系統環境配置如下：

CPU:Intel Xeon E5-2650；內存:96 GB；

眾核加速卡:Intel MIC卡(61核,8 G共享存儲);

操作系統：Redhat Linux 6.3 X64；

編譯器:支持Intel MIC應用開發的Intel編譯器。

5.2 測試結果及分析

測試數據采用稀疏矩陣集(UFget)[17]，文中所提算法是基于對角占優線性方程組的求解，選取的用于測試算法的稀疏矩陣集見表1。

表1 用于測試算法的稀疏矩陣集列表

圖2給出了DLB_GaBP算法與并行GaBP算法在相同的計算環境中的計算時間對比，可以看出DLB_GaBP算法具有更好的計算效率。

圖2 DLB_GaBP與GaBP算法計算時間對比

圖3給出了DLB_GaBP算法在相同運算規模而不同的運算線程或核心數環境中的加速比對比，從圖中可看出隨著算例計算規模的加大加速比也明顯增大，表明算法具有良好的線程數或核數變化的加速效果。

圖3 DLB_GaBP在不同線程數或核數中的加速比

圖4是DLB_GaBP與GaBP算法前后加速比對比圖，算例矩陣名稱是按其非零元素規模大小排列。

從圖中可清楚看出，兩種算法的加速比隨著非零元數量的增大而提高，尤其是非零元素規模不大的環境中，DLB_GaBP的加速優勢并未呈現出來，隨著非零元素個數的增加，即計算問題數據規模的增大，DLB_GaBP算法明顯比GaBP算法有更好的加速比。

結合用于測試算法的稀疏矩陣集列表1，觀察圖4的DLB_GaBP算法的加速比對比曲線，稀疏矩陣的行數小于104，非零元小于105的小規模計算問題，算法并沒有加速效果，這是由于多核并行算法在啟動多核并行占用的時間遠比小規模算例實際求解計算的時間長。這說明問題求解的規模越大，算法的優越性越明顯。

圖4 DLB_GaBP與GaBP算法加速比對比

6 結束語

通過建立列數相同的列塊順序進行存儲格式定義劃分，以確保各線程進行計算時的列塊中的列數相等的基礎上，實現在有限的內存空間對大規模稀疏線性方程組的左端矩陣存儲，可降低算法的空間復雜度。同時通過挖掘基于GaBP的迭代加速優化方法本身的特性，建立具有負載均衡特性的DLB_GaBP算法，具有適應性更好、迭代收斂速度更快的優勢，更適用于大規模數值問題的計算，具有一定的推廣應用價值。