馮 佩， 鐘 誠(chéng)， 韋 偉

(廣西大學(xué)計(jì)算機(jī)與電子信息學(xué)院，廣西南寧 530004)

0 引 言

線性方程組的求解是數(shù)值代數(shù)的基本問題之一，許多工程與科學(xué)計(jì)算問題最終都可歸結(jié)為求解大規(guī)模線性代數(shù)方程組[1]。Gauss-Seidel算法是求解線性方程組的一種經(jīng)典方法，它由Jacobi算法演化而來[2]。已有的求解線性方程組的Gauss-Seidel并行算法大多是基于SIMD模型或者M(jìn) IMD-SM模型設(shè)計(jì)的[3-5]，并且針對(duì)的是方程組的系數(shù)矩陣，是稀疏矩陣的情形[6]，而對(duì)于大規(guī)模的矩陣求解可能會(huì)造成額外的同步開銷或通信開銷?；诙嗪颂幚砥骱虵PGA混合結(jié)構(gòu)的系統(tǒng)[7]，近年來已成為一種新型的計(jì)算平臺(tái)。與單核處理器相比，多核處理器能夠提供更強(qiáng)的并行處理能力，支持線程級(jí)并行和緩存高效[8，9]。

本文充分利用多核處理器的多級(jí)存儲(chǔ)結(jié)構(gòu)和多線程并行技術(shù)，研究在多核計(jì)算機(jī)上設(shè)計(jì)實(shí)現(xiàn)存儲(chǔ)高效、可擴(kuò)展的并行求解n階線性方程組的Gauss-Seidel算法，并與已有算法進(jìn)行實(shí)驗(yàn)性能對(duì)比。

1 增廣矩陣的存儲(chǔ)方式

設(shè)線性方程組A X=B，系數(shù)矩陣為A，增廣矩陣為(A|B)，方程可以表示為:

求解n階線性方程組Gauss-Seidel算法的計(jì)算方法[10]如下:

設(shè)矩陣規(guī)模為A N×N，共享L2Cache的大小為C2，私有L1Cache的大小為C1，處理核數(shù)目為p。根據(jù)多核計(jì)算機(jī)的多級(jí)存儲(chǔ)層次，將增廣矩陣(A|B)按行劃分的方法如下:

(1)將主存中的N行矩陣數(shù)據(jù)依據(jù)L2Cache的大小進(jìn)行分組，計(jì)算出每組中矩陣的行數(shù)。

對(duì)主存中數(shù)據(jù)進(jìn)行分組時(shí)，假設(shè)主存中每組的大小為αC2，則可將主存分成組，其中N(N+ 1)為增廣矩陣(A|B)中元素個(gè)數(shù);α為L(zhǎng)2Cache的利用率(α=0.6，0.7，0.8，0.9)。

記每組中的行數(shù)為L(zhǎng)，每組中有L(N+1)個(gè)元素，則

(2)對(duì)每組中的L行矩陣數(shù)據(jù)進(jìn)行分塊，并計(jì)算出每塊中的矩陣規(guī)模。

假設(shè)每個(gè)L1Cache的可用大小為βC1，則每塊的大小不超過 pβC1，β為 L1Cache的利用率(β=0.5，0.6，0.7，0.8，0.9，1.0)。記每塊中的行數(shù)為n，每塊中有n(N+1)個(gè)元素，則

(3)將每塊中的n行矩陣數(shù)據(jù)進(jìn)行分段，并計(jì)算出每段中矩陣規(guī)模。

記每段的行數(shù)為S，每段中有S(N+1)個(gè)元素，則

2 算法的設(shè)計(jì)與分析

設(shè)多核計(jì)算機(jī)有p個(gè)處理核，old[1，…，N]和new[1，…，N]分別用來存放舊的數(shù)據(jù)值和當(dāng)前計(jì)算出的新值，c[1，…，N]存放新舊數(shù)據(jù)值的比較結(jié)果。

多核多線程并行求解n階線性方程組算法的基本設(shè)計(jì)思想如下:

(1)并行地將初值賦給2個(gè)共享數(shù)old[1，…，N]和new[1，…，N]。

(2)按行劃分方式將增廣矩陣(A|B)進(jìn)行分塊，計(jì)算出分配到每個(gè)處理核的行數(shù)目以及這些行在每組、每塊、每段中的位置，確定出各個(gè)處理核中每段的行在整個(gè)增廣矩陣(A|B)中的開始和結(jié)束位置。

(3)p個(gè)處理核并行地計(jì)算每段中的矩陣行。

(4)根據(jù)下標(biāo)對(duì)當(dāng)前值和舊值求差的絕對(duì)值，并存放在共享數(shù)組c[1，…，N]中。

(5)重復(fù)執(zhí)行步驟(3)、(4)，直到將每組、每塊、每段中的行數(shù)全部計(jì)算出新值為止。

由于算法在執(zhí)行過程中可能有不止一個(gè)處理核讀取old[1，…，N]的值，而本文實(shí)驗(yàn)使用的多核計(jì)算機(jī)是支持并發(fā)讀的，處理核只需要通過自己所要計(jì)算的行的下標(biāo)讀取old[1，…，N]數(shù)組中的值。因此，該算法第(1)步中并行賦初值所需時(shí)間復(fù)雜度為O(1)，第(2)步中計(jì)算劃分后的矩陣塊在每組、每塊、每段中的位置所需的時(shí)間復(fù)雜度為:

第(1)～(5)步所需的時(shí)間復(fù)雜度為:

算法總的時(shí)間復(fù)雜度為:

從(2)式可以看出，第k+1次迭代必須依次求出 x0k+1，x1k+1，x2k+1，…，xn-1k+1，計(jì)算過程具有較強(qiáng)的順序性。當(dāng)j＜i時(shí)，必須等待的計(jì)算完成后才能進(jìn)行計(jì)算，這將會(huì)導(dǎo)致較大的同步開銷;而且各個(gè)處理核在并行執(zhí)行時(shí)相互之間存在數(shù)據(jù)移動(dòng)，即處理核之間要相互訪問，容易造成較大的通信開銷。

為了減少同步等待，本文采用異步并行執(zhí)行方式，將方程組的增廣矩陣按多級(jí)存儲(chǔ)方式存儲(chǔ)到p個(gè)處理核的L1Cache中，p個(gè)處理核并行地對(duì)各自分配到的矩陣行進(jìn)行計(jì)算，同時(shí)充分利用共享的L2Cache，使得處理核之間不存在數(shù)據(jù)移動(dòng)，可大大減少通信開銷，因此多核多線程并行求解線性方程組算法的運(yùn)行速度比原Gauss-Seidel并行算法更快。

3 實(shí) 驗(yàn)

3.1 實(shí)驗(yàn)環(huán)境

實(shí)驗(yàn)多核計(jì)算機(jī)硬件平臺(tái)為主存容量2 GB、二級(jí)緩存容量12 GB、每個(gè)處理器核心的一級(jí)緩存容量32 kB的Intel(R)Xeon(R)E5405 2.0 GH z4核計(jì)算機(jī)，在 Red Hat Enterprise Linux 5操作系統(tǒng)支持下，采用OpenMP和C語(yǔ)言編程實(shí)現(xiàn)多核多線程并行求解n階線性方程組算法。

3.2 方程組系數(shù)矩陣

采用的方程組系數(shù)矩陣A和向量B為:

其中，i，j=1，2，…，N。實(shí)驗(yàn)中采用 x-0[i]初值為零向量，精度ε=0.01。

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

在擁有4個(gè)處理核的計(jì)算機(jī)上，使用動(dòng)態(tài)關(guān)閉處理核的方法，本文采用不同矩陣規(guī)模、不同處理核數(shù)目、不同線程數(shù)目，分別測(cè)試并行求解線性方程組算法的運(yùn)行時(shí)間、加速比和可擴(kuò)展性;測(cè)試了二級(jí)緩存的利用率對(duì)算法性能的影響，并且將本文給出的多核多線程并行求解n階線性方程組算法和原Gauss-Seidel并行算法進(jìn)行了對(duì)比。

圖1所示給出了階數(shù)為600、1 000、2 000、3 000和3 500，分別在單機(jī)多核計(jì)算機(jī)上處理核數(shù)為4核、3核、2核運(yùn)行時(shí)，隨著線程數(shù)目逐漸增加，本文的并行求解線性方程組算法所需的運(yùn)行時(shí)間。

圖1 本文算法在不同核數(shù)上的運(yùn)行時(shí)間

圖1的實(shí)驗(yàn)結(jié)果表明，運(yùn)行不同處理核數(shù)目，隨著線程數(shù)目的改變，本文給出的并行求解線性方程組算法的運(yùn)行時(shí)間各不相同;同時(shí)也可以看出，對(duì)于不同處理核數(shù)，其最優(yōu)線程數(shù)目也不同，一般均為處理核數(shù)目的2～3倍左右。當(dāng)線程數(shù)目超過最優(yōu)線程數(shù)時(shí)，算法的運(yùn)行時(shí)間會(huì)上升，這是因?yàn)楫?dāng)線程數(shù)目增加時(shí)，會(huì)造成線程對(duì)處理器核的競(jìng)爭(zhēng)，同時(shí)線程的啟/停以及上下文切換也會(huì)增加一些額外的時(shí)間開銷。

通過圖1的對(duì)比也可以看出，隨著線程數(shù)目的變化，線性方程組矩陣規(guī)模越大，算法運(yùn)行時(shí)間變化幅度越明顯，這是因?yàn)榫€性方程組矩陣規(guī)模較大時(shí)，算法對(duì)緩存的利用率提高，進(jìn)而提升本文給出的并行求解線性方程組算法的效率。

不同矩陣規(guī)模，隨著處理核數(shù)增加，采用最優(yōu)的線程數(shù)運(yùn)行時(shí)，本文算法的加速比和等效率曲線，如圖2所示。

圖2 本文算法的加速比和等效率曲線

圖2a給出了隨著處理核數(shù)的增加，不同矩陣規(guī)模采用最佳線程數(shù)運(yùn)行時(shí)，本文算法并行求解方程組所獲得的加速比。圖2a的實(shí)驗(yàn)結(jié)果表明，采用最優(yōu)的線程數(shù)，隨著處理核數(shù)的增多，運(yùn)行不同規(guī)模的線性方程組矩陣，本文給出的并行求解線性方程組Gauss-Seidel算法的加速比會(huì)略有增加，線性方程組矩陣規(guī)模越大，加速比的增加幅度越明顯。

圖2b的實(shí)驗(yàn)結(jié)果表明，當(dāng)處理核數(shù)固定時(shí)，并行算法效率隨著線性方程組矩陣規(guī)模增大而增大，處理核數(shù)越多，效率變化趨勢(shì)越明顯;同時(shí)可以看出當(dāng)并行執(zhí)行效率不變時(shí)，矩陣規(guī)模隨著處理核數(shù)的增加而亞線性增加，這說明本文提出的多核多線程并行求解線性方程組算法是可擴(kuò)展的。

運(yùn)行4個(gè)處理核，不同矩陣規(guī)模下，L2Cache利用率增加時(shí)，本文算法所需的時(shí)間如圖3所示。

圖3 L2Cacheα增加時(shí)本文算法所需時(shí)間

圖3的實(shí)驗(yàn)結(jié)果表明，L2Cache的不同利用率，對(duì)算法性能有一定影響。當(dāng)運(yùn)行4個(gè)處理核時(shí)，對(duì)于不同規(guī)模的線性方程組矩陣，從整體上來看，當(dāng)L2Cache的利用率α=0.5時(shí)，也即使用二級(jí)緩存容量的1/2時(shí)，本文給出的并行求解線性方程組Gauss-Seidel算法獲得的效果最佳。由此可以看出，選擇合適的緩存大小，合理地利用緩存，對(duì)并行求解線性方程組算法的性能會(huì)有所改善。

運(yùn)行4個(gè)處理核，不同矩陣規(guī)模下，本文算法和原Gauss-Seidel并行算法的時(shí)間對(duì)比，如圖4所示。

圖4 本文算法和原Gauss-Seidel并行算法時(shí)間對(duì)比

圖4的實(shí)驗(yàn)結(jié)果表明，當(dāng)線性方程組矩陣規(guī)模較小時(shí)，原Gauss-Seidel并行算法的執(zhí)行時(shí)間和本文給出的并行求解線性方程組算法的執(zhí)行時(shí)間相當(dāng)，但當(dāng)矩陣規(guī)模增大時(shí)，本文算法的執(zhí)行時(shí)間明顯少于原Gauss-Seidel并行算法。這說明，當(dāng)矩陣規(guī)模較大時(shí)，在多核計(jì)算機(jī)上運(yùn)行本文的并行求解線性方程組算法更有效。

4 結(jié)束語(yǔ)

本文將線性方程組的增廣矩陣按行進(jìn)行分塊并存儲(chǔ)在多核計(jì)算機(jī)各級(jí)緩存中，多個(gè)處理核并行執(zhí)行多線程對(duì)各段的矩陣行進(jìn)行處理，設(shè)計(jì)實(shí)現(xiàn)了高效、可擴(kuò)展的線程級(jí)并行求解線性方程組Gauss-Seidel算法，并通過實(shí)驗(yàn)測(cè)試獲得不同的增廣矩陣規(guī)模、運(yùn)行處理核數(shù)和并行線程數(shù)目的最佳匹配關(guān)系以及各級(jí)緩存利用的最佳方案。

下一步的研究工作將在同構(gòu)/異構(gòu)多核機(jī)群系統(tǒng)上，研究設(shè)計(jì)存儲(chǔ)高效、通信高效、加速比高、可擴(kuò)展性好、進(jìn)程級(jí)和線程級(jí)并行的線性方程組求解算法。

本文初稿首次刊登于《計(jì)算機(jī)技術(shù)與應(yīng)用進(jìn)展?2010》

[1] 尚月強(qiáng)，楊一都.基于PVM的稠密線性方程組網(wǎng)上并行求解[J].計(jì)算機(jī)工程與設(shè)計(jì)，2006，7(9):1591-1594.

[2] Salkuyeh D K.Generalized Jacobiand Gauss-Seidelmethods for solving linear system[J].Jou rnal ofAp plied Mathematicsand Computing，1998，5(2):341-349.

[3] Malek M.Undirected graphsmodels for system-level fau lt diagnosis[C]//Pro ceedingsof 7th Symposium on Com pu ter A rchitectures，1980:1-35.

[4] BianchiniR P，BuskensW.Implemen tation of online distribu ted sy stem-level diagnosis theory[J].IEEE T ransactions on Compu ters，1992，41(3):616-625.

[5] K rzy sztof D，And rzej P.G lobally optimal diagnosis in systems with random fau lts[J].IEEE T ransactions on Compu ter，1997，46(2):200-204.

[6] 李曉梅，蔣增榮.并行算法[M].長(zhǎng)沙:湖南科學(xué)技術(shù)出版社，2001:7-54.

[7] 胡文彬，吳劍旗，洪 一.多FPGA驗(yàn)證平臺(tái)引腳限制的解決方案[J].合肥工業(yè)大學(xué)學(xué)報(bào):自然科學(xué)版，2010，33(10): 1519-1522.

[8] Akhter S，Roberts J.Mu lti-core programm ing:increasing performance through softw are mu lti-threading[M].北京:電子工業(yè)出版社，2007:21-71.

[9] 周偉明.多核計(jì)算與程序設(shè)計(jì)[M].武漢:華中科技大學(xué)出版社，2009:4-124.

[10] 陳國(guó)良.并行算法的設(shè)計(jì)與分析[M].第3版.北京:高等教育出版社，2009:411-475.