集群下Cholesky分解的核外預取算法

劉 鳳,劉青昆

(遼寧師范大學 計算機與信息技術學院，遼寧 大連 116081)

隨著科學技術的迅猛發展，人們需要處理的數據量迅速增長。大規模并行應用涉及的數據量是非常驚人的，內存容量常常不能滿足涉及到大數據量計算問題的存儲需求。因待處理數據無法全部讀入內存，所以只能保存數據到外部磁盤系統。當程序運行時，某個時刻只能將部分數據調入內存并參加計算，并且在適當的時候寫回外存，通過某種策略實現內外存數據的交換。由于運算過程中數據并沒有全部存放在內核存儲器中，所以稱之為核外計算。核外計算程序中包含大量的文件操作，因訪問磁盤數據的速度比較慢，所以在處理大數據量問題時，I/O性能顯然要超過CPU性能而成為重要的限制因素。因此，如何對核外數組進行合理調度與高效訪問是解決核外計算應用問題的關鍵。

利用三角分解式求解對稱正定方程組是一種有效方法。單行卷簾存儲Cholesky分解[1]使每個節點所分配的數據最多相差一行，但是這種劃分方式的通信開銷比較大。多行卷簾存儲Cholesky分解[2]可以減少通信開銷,充分利用節點的緩存，但是沒有充分考慮緩存的效率問題，因為緩存可以對程序性能帶來巨大的改變。遞歸算法[3]通過將矩陣分塊使各個子矩陣的運算能夠在高速緩存中進行，以提高運算效率，同時遞歸算法良好的數據局部性和矩陣遞歸的分塊，使其非常適合分層多級存儲的計算機結構。但是在遞歸調度算法中，非對角塊的運算不能充分利用計算機的分級存儲結構。分塊Cholesky分解[4-5]通過精心的挑選塊的大小，使每個塊都能充分地利用一級緩存，并且合理地劃分塊的大小還能使每個塊常駐內存。參考文獻文[4]充分利用了系統硬件的并行機制以使通信與計算重疊，減少了節點的等待時間。同時通過矩陣元素的重排使每個塊都被連續存儲，以充分利用計算機的多層次存儲結構，減少數據拷貝和內部的變化。但是這些算法只在數據量較小的情況下適用；當數據量較大時，數據將無法一次性讀入內存，因此引進了核外計算的概念。

本文對Cholesky分解的并行算法進行了深入的研究，給出的核外算法通過預取[6-7]的方法使得數據在實際使用之前從硬盤移到緩存中,從而進一步提高了緩存的命中率，減少了文件讀取的時間。

1 基本概念

1.1 Cholesky分解

Cholesky分解用于求解系數矩陣對稱正定的線性方程組。

式(1)中，A為實對稱正定矩陣，且A的所有順序主子式均不為零，X和b為矩形矩陣或向量。為了避免平方根法Cholesky分解 (A＝LLT)的開方運算且擴大使用范圍，將矩陣對稱正定矩陣A分解成A＝LDLT。當A是正定對稱矩陣時，A的分解是唯一的。其中D是對角線元素全為正的對角矩陣，L是單位下三角矩陣。其求解公式如下：

1.2 預取方法

在核外計算中數據存儲在磁盤上，只有CPU對數據進行處理時才將其讀入內存緩沖區，處理完成后寫回文件。這樣，核外計算過程就可以描述成：讀入核外數組文件的一塊數據到內存緩沖區，進行計算等處理，處理完成后將其寫回數組文件；讀取下一塊數據，進行數據處理，處理完成后寫回文件。如此反復執行，直到整個核外數組處理完成后結束。顯而易見，這是一個順序流水線，執行I/O操作與CPU進行計算是順序執行的；由于二者在執行前需彼此等待，所以程序的執行效率非常低。于是考慮到，當I/O操作的數據與CPU執行計算的數據無相關性時，可以把對下一個數據塊的I/O操作與對當前數據塊的計算重疊起來，這樣就可以隱藏I/O的延遲，達到縮短整個程序運行時間的目的。圖1所示為數據預取前后的對比。

本算法采用雙緩沖區的方式，即為一個核外數組分配兩個內存緩沖區buffer[2]，輪流存放本次和下一次讀入的數據塊。用預取的方法對第s個子文件更新，其偽碼如下：

圖1 數據預取前后的對比

(1)主節點將第一個子文件讀入內存數組buffer[0]中，并將其廣播出去；

(2)依次用已分解的子文件更新第s個子文件，即對ka=1,2,……,s-1循環：

若該節點是主節點,做以下工作：①讀入第ka個文件到buffer[ka%2]中；②發送第ka個文件;

若該節點是從節點,則做以下工作：①用第ka-1個文件更新第s個文件；②接收第ka個文件。

(3)若該節點是從節點，則用第s-1個文件更新第s個文件。

2 用預取算法實現Cholesky分解

2.1 Cholesky分解并行算法的思想

將大規模矩陣A連續拆分成q個子文件，并將其儲存到主節點中。主節點一次調入內存一個子文件，依次對各子文件進行Cholesky分解，直至最后一個子文件分解完畢。所以只要拆分的子文件大小不超過空閑內存的范圍，算法就可以運行。其算法流程見圖2，具體步驟如下：

(1)文件拆分過程

由主節點機進行矩陣拆分存儲的操作，打開存儲于硬盤中的A.txt文件，依內存容量的需要進行拆分，讀取相應規模的數據，生成子文件，并將這些子文件j.txt(j=1,2,3)存儲回硬盤。

(2)A陣的Cholesky分解過程

圖2 Cholesky分解流程圖

考慮到Cholesky分解過程中的依賴關系，對其由上至下依次進行分解。由于在三角分解中，既需要分解后的上三角，又需要分解后的下三角，所以可以將分解后的上三角復制到下三角。這樣既可以減少運算量又可以減少存儲量。在分解過程中利用其并行性，主節點打開第一個子文件 1.txt后，賦值給 a6×15，再將 a廣播出到各個從節點。每個從節點經卷簾計算出分解因子后，廣播出去，如此地計算kb=2次，該子文件的Cholesky分解完畢。當1.txt計算結束后,主節點打開下一個子文件2.txt，賦值給a6×15，再將 a廣播到各個從節點。同時主節點還要打開分解完畢的子文件 1.txt,賦值給 aa6×15(即 buffer[0]),再將aa廣播到各個從節點。各個從節點用aa6×15卷簾去更新 a6×15，更新完畢后 a6×15再像上面的操作，每個從節點卷簾的計算出分解因子后，廣播出去，如此的計算kb=2次后，主節點打開下一個子文件3.txt，進行如文件2.tx的操作。

2.2 Cholesky分解并行算法的描述

將A連續拆分存儲到q個子文件中，依次對各子文件進行Cholesky分解，直至對最后一個子文件分解后，原大規模系數矩陣已被拆分到幾個小文件中存儲。將節點機分別標記為 P0，P1…Pnp，np代表從節點機個數，P0代表主節點，kb表示在一個子文件中單個節點機計算的行數，s是子文件依次分解的循環控制變量，ka表示讀取已分解子文件的循環控制變量，在拆分后 a∈R[kb×np]×n，aa∈R[kb×np]×n。 具 體 步 驟 如 下 ：

(1)文件拆分

主節點機打開存于硬盤中的A.txt文件，按內存容量的需要進行拆分，讀取相應規模的數據，生成q個子文件。

(2)A矩陣的Cholesky分解過程

①1.txt(s=0)的 Cholesky分解：(a)主節點打開 1.txt，讀取文件內容，賦值到a[np*kb][n]數組中，并將數組a廣播出去。(b)從節點 myid(myid=1,2,……,np)將數組 a中的第 i行(i%np=myid-1)進行 Cholesky分解。(c)從節點將已分解的對角元素賦給數組d[n]。(d)從節點用數組a已分解的元素重寫其下三角。(e)從節點將分解后的數組a寫回主節點。

②文件 j.txt(s=1,2,……,q-1,j=s+1)的 Cholesky分解：(a)主節點打開并讀取 j.txt內容，賦值到a且將其廣播出去。(b)用預取的方法對第j個子文件更新。即主節點依次打開已分解文件 ka.txt(ka=0,1,…,s-1)，賦值給buffer[ka%2]，并將 buffer[ka%2]廣播出去。在各個從節點根據buffer[ka%2]的值更新數組a值的同時，主節點讀取下一個已分解的文件。(c)用數組buffer[ka%2]的元素重寫數組a的下三角。(d)將數組a進行Cholesky分解。(e）將數組a本次分解的對角線元素賦值給數組d[n]。(f)用數組a本次分解的元素重寫其下三角。(g)將分解后的數組a寫回主節點。

2.3 算法的復雜度分析

對于點對點的通信，測量開銷使用乒乓法：節點0發送m個字節給節點1；節點1從節點0接收m個字節后，立即將消息發回節點0。總的時間除以2，即可得到點到點通信時間，也就是執行單一發送或接收操作的時間。通信開銷的解析表達式是消息長度m(字節)的線性函數：tcomm(m)=Ts+Tb×m;其中 Ts表示通信的啟動時間,Tb表示發送每個字節所需時間 (它是帶寬的倒數)。由于MPI_Bcast函數采用樹算法，所以一次MPI_Bcast的通信開銷為 tcomm(s)×logP。Cholesky分解的執行時間可分為子文件的更新時間和自身分解時間。因此，Cholesky分解時間：Tc=Tg+Tf，其中 Tg是子文件更新時間，Tf是自身的分解時間。在文件的更新時間中,又細化為讀文件時間Td、計算時間Tj和通信時間Tt。為了減少分解的執行時間，本文通過預取的方法使更新過程中的讀文件與計算重疊，使得 Cholesky分解的時間開銷：Tc

子文件的大小是通過參數kb(1≤kb≤n/np)而設定的，當kb=n/np時相當于無文件劃分并行求解；當kb=1時，在每個子文件中每個節點機只計算一行，此時的通信量最大(按照這種分配方法劃分數據后，每個處理器上須存儲的內容為kb×np×n規模的矩陣 A)。因此算法的并行執行時間：Tn=T1+Tc+T2，其中，T1為劃分子文件消耗的時間，Tc是Cholesky分解時間，T2為冗余的控制和管理開銷。

3 實驗測試與結果分析

實驗環境采用由5臺主頻為2.8 GHz的Intel Xeon CPU，內存為ECC DDR-2 SDRAM 2 GB的Dell PowerEdge 2850構建的集群。該集群運行Linux RH9操作系統，并且建立了MPI并行編程環境。本文測試的兩個程序分別為：帶狀循環劃分的核外程序out和帶狀循環劃分的核外預取程序pre。

表1表示當A的階數n=2 478,核外Cholesky分解各個部分的執行時間對比。表2表示了在節點數np分別為2,3,4,5時，Cholesky分解中的子文件更新各個部分的時間對比(kb=10,n=2 479)。圖3為核外預取程序相對于核外程序更新時間的效率提高百分比。

表1 矩陣規模n=2 478,Cholesky分解各部分執行時間

由表1可見，在求解過程中Cholesky分解中的更新部分占了大部分時間，所以提高更新部分的執行率能明顯縮短總的執行時間。

表2 out和pre的Cholesky分解更新時間對比

圖3 pre相對于out的更新時間效率提高百分比

由表2可見，預取方法使得子文件的更新時間平均縮短了15%左右。

圖3表明，核外預取分解并行算法的效率明顯好于核外分解并行算法，并且隨著節點個數的增加，提高的效率先增大后減小。這是因為隨著節點數目的增加，子文件容量增加，計算量增加，計算與文件讀取的重疊時間越多。在節點數np=3時，效率提高幅度最大。當節點數np>3時，隨著節點數的增加通信時間急劇增加，而計算時間逐漸減少，計算與文件讀取的重疊時間減少，所以執行效率提高的幅度降低。

本文通過對Cholesky分解并行算法的研究，表明預取方法是計算核外稠密線性方程組的有效方法，非常有利于縮短I/O與CPU速度間的差距。數據預取的方法可以做到計算與I/O并行，即在計算一塊數據的同時讀入下一塊要處理的數據。為了存放預取的數據，本文采用雙緩沖區的方式，即為一個核外數組分配兩個內存緩沖區，輪流存放本次和下一次讀入的數據塊。此外，預取方法同樣適用于QR分解、LU分解、高斯削去等其他線性代數問題。當然采用數據預取技術時，I/O操作時間與CPU執行時間的比例是保證預取效果的關鍵。

事實上，Cholesky分解也可以采用數據重用的方法。通過減少對子文件的讀取次數，可以進一步提高Cholesky分解算法的效率。由于目前分布存儲計算機的處理速度都很高，而其網絡通信速度較慢,所以用增加計算和通信粒度的辦法來降低通信成本。

[1]遲學斌.Transputer上Cholesky分解的并行實現[J].計算數學，1993(3):289-294.

[2]王順緒，周樹荃.卷簾行存儲下的一種并行Cholesky分解及其在PAR95上的實現[J].南京航空航天大學學報，1999，31(4)：428-432.

[3]陳建平.Jerzy Wasniewski，Cholesky分解遞歸算法與改進[J].計算機研究與發展，2001，38(8):923-926.

[4]GUSTAVSON F G,KARLSSON L,KAGSTRO B M.Distributed SBP cholesky factorization algorithms with nearoptimal scheduling[J].ACM Transactions on Mathematical Software,2009,36(2):1-25.

[5]ANDERSEN B S,GUNNELS J A,GUSTAVSON F G,et al.A fully portable high performance minimal storage hybrid format cholesky algorithm[J].ACM Transactions on Mathematical Software,2005,31(2):201-227.

[6]丁文魁，汪劍平，向華，等.p-HPF并行編譯系統核外計算的實現及優化策略[J].計算機學報，1999，22(10):1042-1049.

[7]姚維.Linux下一種磁盤節能的預取算法[J].計算機系統應用，2010,19(7):91-94.