劉 侃，王欣亮，許 平，薛 巍,2+

1.清華大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)系，北京100086

2.國(guó)家超級(jí)計(jì)算無(wú)錫中心，江蘇 無(wú)錫214100

+通訊作者E-mail:xuewei@mail.tsinghua.edu.cn

1 引言

計(jì)算機(jī)數(shù)值模擬是科學(xué)和工程領(lǐng)域廣泛使用的研究方法。通過(guò)將真實(shí)的問(wèn)題經(jīng)過(guò)數(shù)學(xué)建模，可以用數(shù)值計(jì)算的方式對(duì)這些問(wèn)題進(jìn)行模擬。在計(jì)算機(jī)中，連續(xù)的真實(shí)世界往往用有限、離散的數(shù)據(jù)來(lái)表示。隨著人們對(duì)數(shù)值模擬精度的需求越來(lái)越高，高分辨率的數(shù)值模擬成為必需。而高分辨率帶來(lái)的大數(shù)據(jù)量和大計(jì)算量，則對(duì)數(shù)值計(jì)算的性能提出了非常高的要求。

偏微分方程是真實(shí)問(wèn)題中常見的一類數(shù)學(xué)模型。對(duì)于含時(shí)的偏微分方程，主要的兩種求解方法是顯式方法和隱式方法。其中，隱式方法可以將微分方程的求解轉(zhuǎn)化為一系列線性方程組的求解。而對(duì)于一維空間的問(wèn)題，隱式方法中最常見的線性方程組便是三對(duì)角方程。而在二維和三維空間的問(wèn)題的求解中，常有需要求解大量一維子問(wèn)題的場(chǎng)景，如交替方向隱式方法[1-2]、迭代求解器的預(yù)調(diào)節(jié)器[3]等。因此，三對(duì)角方程求解器廣泛地被應(yīng)用于如計(jì)算機(jī)圖形學(xué)[1,4-5]、流體力學(xué)[2,4,6]、三次樣條計(jì)算[7]等眾多科學(xué)和工程領(lǐng)域。為了解決實(shí)際應(yīng)用中在高維空間中求解大量一維子問(wèn)題的需要，本文所討論的范圍是求解大量的、規(guī)模不太大的三對(duì)角方程的情況，求解器最重要的性能指標(biāo)是吞吐率，即大量三對(duì)角方程的總未知數(shù)數(shù)量與總求解時(shí)間之比。

目前，對(duì)于三對(duì)角方程，中央處理器（central processing unit，CPU）、圖形處理器（graphics processing unit，GPU）等主流硬件平臺(tái)上都提出了高度優(yōu)化的并行算法，同時(shí)有基礎(chǔ)數(shù)學(xué)庫(kù)支持。但是對(duì)于申威26010眾核處理器，還沒(méi)有一種算法能有效地利用其硬件特性來(lái)達(dá)到最大化的性能。針對(duì)這一現(xiàn)狀，本文綜合考慮的申威處理器的硬件特性，從計(jì)算量、訪存量和并行度等方面分析追趕法（Thomas algorithm）、循環(huán)消去（cyclic reduction，CR）算法和并行循環(huán)消去（parallel cyclic reduction，PCR）算法在申威處理器上可能達(dá)到的性能，設(shè)計(jì)了一種面向申威處理器的分布式CR 算法swDCR（Sunway-oriented distributive cyclic reduction）。其吞吐率相比主核上的追趕法達(dá)到了單精度43.9 倍和雙精度36.7 倍的加速，相比從核上的追趕法達(dá)到了單精度和雙精度均2.07倍的加速。

2 背景介紹

2.1 申威處理器架構(gòu)介紹

申威26010 眾核處理器是中國(guó)第一款自主研發(fā)的眾核處理器?；谏晖?6010 處理器的神威太湖之光超級(jí)計(jì)算機(jī)，在2016 年到2017 年間曾是世界上計(jì)算能力最高的超級(jí)計(jì)算機(jī)，至2018 年11 月世界排名仍居第三。其理論計(jì)算能力達(dá)125 PFlops，Linpack[8]實(shí)測(cè)計(jì)算能力達(dá)到93 PFlops。

申威26010 眾核處理器架構(gòu)如圖1 所示。一個(gè)處理器包含4個(gè)核組，每個(gè)核組具有765 GFlops的峰值計(jì)算能力和34.1 GB/s的理論訪存帶寬。每個(gè)核組包含一個(gè)主核（management processing element，MPE）和64 個(gè)從核（computation processing element，CPE），從核通過(guò)8×8網(wǎng)格的形式組織成從核陣列。

Fig.1 Architecture of Sunway 26010 many-core processor圖1 申威26010 眾核處理器架構(gòu)

每個(gè)主核具有一個(gè)32 KB 一級(jí)指令緩存、一個(gè)32 KB 一級(jí)數(shù)據(jù)緩存和一個(gè)256 KB 二級(jí)緩存，每個(gè)從核具有一個(gè)16 KB 指令緩存和64 KB 便箋式存儲(chǔ)器（scratch pad memory，SPM）。除了SPM 以外的緩存都是由硬件控制，并對(duì)軟件透明的。SPM 的訪問(wèn)延遲和帶寬與主核的一級(jí)緩存相當(dāng)，但區(qū)別在于軟件可控。

從核對(duì)主存的訪問(wèn)有兩種模式：全局讀入/寫出（global load/store，gload/gstore）和直接內(nèi)存訪問(wèn)（direct memory access，DMA）。gload/gstore完成主存和寄存器之間的數(shù)據(jù)傳輸，帶寬較低，適合細(xì)粒度的數(shù)據(jù)訪問(wèn)；DMA 完成主存和SPM 緩存之間的數(shù)據(jù)傳輸，在訪問(wèn)大段連續(xù)數(shù)據(jù)時(shí)帶寬較高，適合粗粒度的數(shù)據(jù)訪問(wèn)。根據(jù)Xu等的測(cè)試[9]，用一個(gè)核組64個(gè)從核運(yùn)行Stream Triad測(cè)試程序，測(cè)得DMA的帶寬為22.6 GB/s，而gload/gstore的帶寬只有不到1.5 GB/s。

在從核組成的8×8陣列上，同行或同列的兩個(gè)從核之間可以通過(guò)寄存器通信高速互傳數(shù)據(jù)。每個(gè)從核具有一個(gè)寄存器發(fā)送緩沖區(qū)、一個(gè)寄存器行接收緩沖區(qū)和一個(gè)寄存器列接收緩沖區(qū)。在進(jìn)行寄存器通信時(shí)，發(fā)送方從核會(huì)將數(shù)據(jù)放入發(fā)送緩沖區(qū)，并由核間互連將數(shù)據(jù)傳輸?shù)浇邮辗綇暮说男?列接收緩沖區(qū)中。在接收數(shù)據(jù)時(shí)，接收方從核會(huì)從行/列接收緩沖區(qū)取出數(shù)據(jù)。利用寄存器通信可以實(shí)現(xiàn)核間數(shù)據(jù)的高帶寬低延遲傳輸，根據(jù)Xu等的測(cè)試[9]，兩個(gè)從核之間的點(diǎn)到點(diǎn)通信延遲不超過(guò)11 個(gè)指令周期，單個(gè)核組的從核陣列上的寄存器通信集合帶寬可以高達(dá)637 GB/s。寄存器通信沒(méi)有低層緩存的訪問(wèn)開銷，無(wú)論是延遲還是帶寬，都遠(yuǎn)遠(yuǎn)優(yōu)于目前主流硬件平臺(tái)上通過(guò)共享緩存實(shí)現(xiàn)的核間通信。

申威處理器的寄存器通信機(jī)制存在來(lái)源混淆問(wèn)題。當(dāng)一個(gè)從核1可能要從同列的從核0和從核2接收數(shù)據(jù)時(shí)，來(lái)自從核0和從核2的數(shù)據(jù)都會(huì)被放入從核1的寄存器列接收緩沖區(qū)中，當(dāng)從核1從寄存器列接收緩沖區(qū)中取數(shù)據(jù)時(shí)便不能區(qū)分其來(lái)源。在必要時(shí)，需要加同步使得兩次寄存器通信先后進(jìn)行，或?qū)l(fā)送源的編號(hào)作為數(shù)據(jù)的一部分進(jìn)行通信。

綜上所述，在設(shè)計(jì)申威處理器上的算法時(shí)需要考慮三個(gè)硬件特性：

SPM 緩存：軟件可控的特性使得它比硬件控制的緩存更靈活，但是需要考慮將哪些數(shù)據(jù)放入緩存中。

DMA：帶寬遠(yuǎn)高于gload/gstore，但只適用于粗粒度的內(nèi)存訪問(wèn)。應(yīng)盡可能避免零散的數(shù)據(jù)訪問(wèn)。

寄存器通信：可以實(shí)現(xiàn)核間高速的數(shù)據(jù)傳遞，但只能在同行或同列間傳輸，而且應(yīng)該避免來(lái)源混淆問(wèn)題。

2.2 三對(duì)角方程求解算法的相關(guān)工作

三對(duì)角方程是一類形如Ax=d的方程組，其中A是三對(duì)角矩陣，d是已知向量，需要求解未知向量x，如圖2（a）所示。在內(nèi)存中，矩陣A的三條對(duì)角線和右端向量d通常是連續(xù)存儲(chǔ)的，并且可以認(rèn)為存在a0=cn-1=0，以保證數(shù)據(jù)對(duì)齊，如圖2（b）所示。

Fig.2 Definition of tridiagonal equations圖2 三對(duì)角方程的定義

三對(duì)角方程的經(jīng)典解法是基于高斯消去法的追趕法[10]。追趕法分為前向分解和后向回代兩個(gè)階段，在前向分解中，通過(guò)式（1），從第一行到最后一行，將方程變?yōu)樯先切问健?/p>

在后向回代過(guò)程中，通過(guò)式（2），從最后一行到第一行，求出解。

如果不考慮向量化和緩存，那么追趕法的計(jì)算量和訪存量都是最小的，因此是最快的串行算法。追趕法在求解單個(gè)三對(duì)角方程時(shí)是完全串行的，雖然在線程級(jí)粗粒度并行上，多個(gè)方程可以并行地求解，但是難以有效地利用向量化等細(xì)粒度并行機(jī)制。

CR算法由Hockney在1965年提出[11]。該算法包含兩個(gè)階段：前向消去階段和后向回代階段。在前向消去階段，將方程的偶數(shù)行提取出，按照式（3）得到一個(gè)規(guī)模折半的方程。

對(duì)新的方程可遞歸地進(jìn)行消去，直到得到一個(gè)兩個(gè)未知數(shù)的方程時(shí)，可直接求解。

在后向回代階段，當(dāng)偶數(shù)行的解已被求出時(shí)，按照式（4）可求得奇數(shù)行的解。

CR 算法可就地進(jìn)行，即用新的a、b、c和d覆蓋內(nèi)存中偶數(shù)行的a、b、c和d的位置。但這會(huì)使新的方程在內(nèi)存中不連續(xù)，影響后續(xù)訪問(wèn)的效率。針對(duì)這個(gè)問(wèn)題，G?ddeke 和Strzodka 討論了數(shù)據(jù)布局的優(yōu)化和利用GPU 共享內(nèi)存的優(yōu)化[12]；Davidson 和Owens討論了GPU上的寄存器打包的優(yōu)化[13]。

在CR算法的基礎(chǔ)上，Hockney 和Jesshope 于1981 年提出了PCR 算法[14]。與CR 算法不同之處在于，PCR算法不僅提取出原方程的偶數(shù)行做消去，也提取出原方程的奇數(shù)行做消去，從而將原方程轉(zhuǎn)化為兩個(gè)規(guī)模折半的方程。將兩個(gè)新方程分別遞歸求解后便可得到全部解，不再需要后向回代。

但是這個(gè)版本的PCR 算法是不實(shí)用的，因?yàn)樗挠?jì)算復(fù)雜度相比追趕法和CR 算法的O(n)提高到了O(nlgn)。為了解決這一問(wèn)題，Sakharnykh 將它與追趕法結(jié)合提出PCR-并行追趕法（PCR-pThomas algorithm）[2,6]，將原來(lái)的lgn層遞歸分解改為常數(shù)層分解，之后并行地用追趕法求解所得到的常數(shù)個(gè)小規(guī)模的方程。這樣使得訪存和計(jì)算次數(shù)大幅減少，而且處理器的向量化運(yùn)算仍然可以得到充分利用，但是多線程并行較困難。Kim等通過(guò)利用GPU上的共享內(nèi)存進(jìn)一步優(yōu)化了內(nèi)存訪問(wèn)，提出了tiled-PCRpThomas 算法，其訪存量與經(jīng)典的追趕法一致[15]。Wang 等通過(guò)將這一算法推廣到英特爾至強(qiáng)融合架構(gòu)，利用寄存器優(yōu)化內(nèi)存訪問(wèn)，并分析了PCR層數(shù)的最佳選擇，提出了register-PCR-(half)-pThomas算法[16]。

另外，Bondeli[17]、Polizzi等[18]和Chang等[19]也分別從不同的角度討論了大規(guī)模三對(duì)角方程的求解，但都是采用分治策略，將它的求解轉(zhuǎn)化成很多小規(guī)模方程的求解。本文所討論的范圍是求解大量的、規(guī)模不太大的三對(duì)角方程的情況。

表1 總結(jié)了追趕法、CR 算法和PCR-pThomas 算法的優(yōu)缺點(diǎn)。

Table 1 Comparison of Thomas,CR algorithm and PCR-pThomas algorithm表1 追趕法、CR算法和PCR-pThomas算法的比較

3 申威眾核架構(gòu)的三對(duì)角方程并行求解算法分析

首先考慮追趕法。在申威處理器的從核上，單精度除法的平均每條指令周期數(shù)（cycles per instruction，CPI）約為加乘的15倍；若將一次除法折合成15個(gè)加乘，則單精度下標(biāo)量計(jì)算訪存比為0.917。同理，雙精度下將除法折合成30 次加乘，計(jì)算訪存比則為0.875。因此在申威處理器的從核上，追趕法是訪存受限的。

然后考慮PCR-pThomas 算法。一方面，申威處理器對(duì)單精度和雙精度浮點(diǎn)運(yùn)算的向量化寬度都只有4 個(gè)浮點(diǎn)數(shù)。如表2 所示，PCR-pThomas 算法中，最多只能相比追趕法減少0.5n次向量除法，同時(shí)還要引入額外的向量混洗（shuffle）操作，向量化給PCR-pThomas 算法帶來(lái)的提升空間小。另一方面，PCR-pThomas 算法和追趕法的訪存量相同，受到同樣的訪存帶寬限制。綜合以上兩點(diǎn)，可以發(fā)現(xiàn)，PCRpThomas算法很難達(dá)到比追趕法更高的性能，因此在這兩者之中選擇更簡(jiǎn)單的追趕法即可。

Table 2 Number of floating point operations in Thomas and PCR-pThomas algorithm表2 追趕法與PCR-pThomas算法的浮點(diǎn)運(yùn)算次數(shù)

追趕法在申威處理器上是訪存受限的，要想達(dá)到比追趕法更高的性能，必須減少訪存量。在追趕法中，4/9的訪存是讀寫中間變量，即前向消去之后的超對(duì)角線和右端向量。注意到申威處理器從核陣列上的緩存是軟件可控的，如果能將這些中間變量放在高速緩存中，那么三對(duì)角方程求解器的性能相比追趕法仍有80%的提升空間。

然而每個(gè)從核上的緩存大小是有限的，要將所有中間變量都保存在緩存中，單個(gè)從核所能求解的方程規(guī)模是有上限的。單個(gè)從核的緩存是64 KB，考慮到其他棧變量需要占用一定的空間，可求解的方程規(guī)模上限在單精度為2 048維，在雙精度下為1 024維；這樣的規(guī)模上限對(duì)于實(shí)際問(wèn)題來(lái)說(shuō)限制太大。對(duì)于更大的規(guī)模就需要更大的緩存，為此可以將多個(gè)從核的緩存聯(lián)合起來(lái)，也就是使用多個(gè)從核來(lái)求解一個(gè)方程。

由于涉及到多核，求解算法必須支持線程級(jí)并行，否則不同核之間只能串行執(zhí)行，每個(gè)核都會(huì)把很多時(shí)間浪費(fèi)在等待其他核上，使得效率較低。追趕法本身不能并行，而PCR-pThomas 算法為了減少計(jì)算量和訪存量，舍棄了線程級(jí)并行的機(jī)會(huì)。因此，應(yīng)該選擇支持線程級(jí)并行的CR算法。

本文以CR算法為原型，設(shè)計(jì)了適合申威處理器的swDCR算法。通過(guò)聯(lián)合最多64個(gè)從核，使得可求解的方程的規(guī)模上限提高到了單精度下131 072 維和雙精度下65 536 維。同時(shí)，利用申威處理器的寄存器通信機(jī)制，使得CR算法的前后依賴數(shù)據(jù)能在不同從核間有效地傳輸。

4 swDCR的實(shí)現(xiàn)與優(yōu)化

4.1 基礎(chǔ)實(shí)現(xiàn)

使用從核緩存中的48 KB來(lái)存儲(chǔ)中間數(shù)據(jù)，其中最初讀入的原始方程占32 KB。在CR算法的每次前向消去時(shí)，將得到的新方程寫到當(dāng)前空閑的空間，同時(shí)將原方程的奇數(shù)行打包到原方程所占用空間的前半部分。在消去完成后，原方程所占空間的后半部分就被騰出用于下一層消去，整個(gè)過(guò)程的中間變量便可以存放在這48 KB中；同時(shí)打包也可使回代過(guò)程的數(shù)據(jù)訪問(wèn)變得連續(xù)，便于利用向量化運(yùn)算。

經(jīng)過(guò)一定次數(shù)的消去后，每個(gè)線程所處理的方程規(guī)模都不超過(guò)4，此時(shí)使用追趕法求解。由于此時(shí)數(shù)據(jù)仍然分布在各個(gè)從核的緩存中，稱這一過(guò)程的算法為分布式追趕法。雖然它用到了多核，但是實(shí)際的執(zhí)行仍然是串行的，不過(guò)由于方程規(guī)模足夠小，這部分消耗的時(shí)間對(duì)整體性能影響很小。在追趕法的前向分解過(guò)程中，每個(gè)線程會(huì)等待前面的線程完成前向分解，并通過(guò)寄存器通信將前向消去之后的超對(duì)角線和右端向量的最后一個(gè)元素傳過(guò)來(lái)。每個(gè)線程完成前向分解后會(huì)將該線程處理部分的超對(duì)角線和右端向量的最后一個(gè)元素傳輸給后方線程。而后向回代過(guò)程也是同樣，先等待后方線程，再求解，然后將解的第一個(gè)元素傳輸給前方線程。

為了更高的吞吐率，解單個(gè)方程時(shí)所用到的從核數(shù)應(yīng)該盡可能得少，同時(shí)每個(gè)從核所處理的問(wèn)題規(guī)模不能超過(guò)其上限nmax。可以根據(jù)輸入規(guī)模來(lái)選擇單個(gè)方程所用的從核數(shù)，如式（5）所示。

若所得到的p大于64，即不可能把CR算法所有中間變量放入緩存中，而必須寫入主存。當(dāng)涉及到主存訪問(wèn)時(shí)，CR 算法訪存量大的問(wèn)題將極大地降低算法的效率。另外，大規(guī)模的三對(duì)角方程也不在本文討論范圍內(nèi)。在這種情況下，由于追趕法的訪存量更小，選擇回退到追趕法。

4.2 向量化

申威架構(gòu)下單精度和雙精度浮點(diǎn)運(yùn)算的向量化寬度都是4 個(gè)浮點(diǎn)數(shù)。為了有效地利用向量化來(lái)加速CR算法的消去與回代過(guò)程，應(yīng)該使得每個(gè)線程處理的原方程行數(shù)是8的倍數(shù)，經(jīng)消去后得到的新方程的行數(shù)是4 的倍數(shù)，以保證所有數(shù)據(jù)的對(duì)齊訪問(wèn)，并避免出現(xiàn)不能向量化的余數(shù)部分。雖然可以通過(guò)在方程前后補(bǔ)上單位矩陣來(lái)使得原始方程規(guī)模是8 的倍數(shù)，但是經(jīng)過(guò)一次消去后得到的新方程的規(guī)模只有原始方程的一半，不一定仍是8的倍數(shù)。

針對(duì)這一問(wèn)題，在swDCR算法每次消去之前，引入“借”的機(jī)制。在每次消去前，假設(shè)線程所處理的原方程的行數(shù)已經(jīng)是8的倍數(shù)，那么在需要時(shí)該線程會(huì)從后方線程“借”方程的4行來(lái)自己處理，從而保證在每次消去時(shí)處理的方程行數(shù)都是8 的倍數(shù)。如果這個(gè)“借”的過(guò)程發(fā)生了，那么在回代之后，也需要將所得到的解“還”給后方線程。

在每一層消去時(shí)，每個(gè)線程將計(jì)算自己所處理的方程行數(shù)和起始行號(hào)。若起始行號(hào)不是8的倍數(shù)，則當(dāng)前線程的前4 行歸前方線程處理；若終止行號(hào)不是8 的倍數(shù)，則后方線程的前4 行歸當(dāng)前線程處理。另外，若從后方線程“借”了4 行，這4 行也會(huì)被打包到原方程所占空間的前半部分，在回代時(shí)還會(huì)訪問(wèn)到。

4.3 線程布局

無(wú)論是swDCR 算法的消去和回代過(guò)程，還是分布式追趕法的求解，都需要每個(gè)線程與其前后的線程進(jìn)行通信。考慮到從核間寄存器通信只能同行或同列傳輸，線程在從核上的布局需要滿足要求：每個(gè)線程的前一個(gè)線程和后一個(gè)線程都在其同行或同列的從核上。一種簡(jiǎn)單的線程布局如圖3（a）所示，其中紫色的線連接需要相互通信的兩個(gè)從核。

Fig.3 Thread layout on CPEs圖3 從核上的線程布局

這樣的線程布局會(huì)引起寄存器通信時(shí)的來(lái)源混淆問(wèn)題。針對(duì)這個(gè)問(wèn)題，在兩次寄存器通信之間加同步當(dāng)然也是一種解決方法，但是這樣的同步開銷其實(shí)是可以避免的?？紤]到從同行接收和從同列接收分別經(jīng)由的是行接收緩沖區(qū)和列接收緩沖區(qū)兩個(gè)不同的緩沖區(qū)，可以設(shè)計(jì)一種線程布局，使得對(duì)于每一個(gè)線程，其前一個(gè)線程和后一個(gè)線程的其中一個(gè)在同行的從核上，另一個(gè)在同列的從核上。新的線程布局如圖3（b）所示。

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

5.1 性能測(cè)試

使用申威處理器的一個(gè)核組進(jìn)行實(shí)驗(yàn)。一個(gè)核組具有765 GFlop/s 的峰值計(jì)算能力和34.1 GB/s 的理論訪存帶寬。實(shí)驗(yàn)會(huì)測(cè)試主核上的追趕法、從核上的追趕法以及不同優(yōu)化前后的swDCR算法。每個(gè)算法會(huì)測(cè)試規(guī)模從128 維到32 768 維的足夠量的三對(duì)角方程，且分別測(cè)試單精度和雙精度的問(wèn)題。每個(gè)數(shù)據(jù)都是在重復(fù)運(yùn)行200 次后取運(yùn)行時(shí)間的平均值，并以總未知數(shù)數(shù)量除以總時(shí)間計(jì)算吞吐率，以總訪存量除以總時(shí)間計(jì)算訪存帶寬。

圖4 展示了不同算法在申威26010 處理器單個(gè)核組上求解單精度和雙精度三對(duì)角方程的吞吐率和帶寬。根據(jù)實(shí)驗(yàn)結(jié)果，主核上的追趕法在單精度和雙精度下都只能獲得1 GB/s 的帶寬，從而吞吐率較低。而從核上的追趕法基本可以獲得21 GB/s 以上的帶寬，但是和本文所提出的swDCR算法相比，由于訪存量較多，吞吐率并不高，約為單精度下6.0×108未知數(shù)/s，雙精度下3.0×108未知數(shù)/s。

Fig.4 Throughput and bandwidth of different tridiagonal solvers圖4 不同的三對(duì)角求解算法吞吐率和帶寬

而本文所提出的swDCR 算法，可以獲得比追趕法更高的24 GB/s 帶寬，是理論帶寬的70%。同時(shí)更少訪存量使得吞吐率進(jìn)一步提高，在不同規(guī)模下平均達(dá)到單精度下1.23×109未知數(shù)/s，雙精度下6.2×108未知數(shù)/s。

圖5展示了不同的優(yōu)化對(duì)swDCR算法產(chǎn)生的性能提升。可以發(fā)現(xiàn)，向量化對(duì)于swDCR 的性能非常關(guān)鍵。在沒(méi)有向量化的版本中，和追趕法相比并無(wú)太大的吞吐率優(yōu)勢(shì)；而在加入向量化優(yōu)化后，性能得到了顯著的提升。在方程規(guī)模較小時(shí)，同步的影響較?。欢诜匠桃?guī)模較大，也就是所用的從核數(shù)較多時(shí)，同步對(duì)性能的影響逐漸顯現(xiàn)出來(lái)。

同時(shí)可以發(fā)現(xiàn)，性能曲線隨著方程規(guī)模的上升，性能會(huì)有階躍式的下降，斷點(diǎn)在4 096、8 192、16 384等位置。這主要是因?yàn)?，在選擇單個(gè)方程所用從核數(shù)時(shí)，在這些位置從核數(shù)翻倍，每個(gè)從核所處理的問(wèn)題規(guī)模有階躍式的減小，從而導(dǎo)致性能階躍式地降低。這一現(xiàn)象在未經(jīng)各項(xiàng)優(yōu)化的版本上最為明顯，而經(jīng)過(guò)各項(xiàng)優(yōu)化后，算法變得越來(lái)越受訪存帶寬限制，這樣的階躍式下降也變得越來(lái)越不明顯。同時(shí)也可以發(fā)現(xiàn)，當(dāng)每個(gè)從核所處理的問(wèn)題規(guī)模逐漸增大時(shí)，性能也逐漸提升。

與主核上的追趕法相比，swDCR 算法在單精度和雙精度下分別可以獲得43.9 倍和36.7 倍的加速。與從核上的追趕法相比，swDCR 算法在單精度和雙精度下均可獲得2.07倍的加速。這樣的加速效果超出了單從訪存角度考慮得到的1.8倍理論值，相比追趕法，CR 算法適合向量化這一優(yōu)勢(shì)能帶來(lái)額外的加速效果。

5.2 精度測(cè)試

swDCR在實(shí)現(xiàn)中是完全保精度的。為了驗(yàn)證其求解精度，選取5 個(gè)16 384 維的三對(duì)角矩陣A，令精確解x的所有元素全為1，計(jì)算右端向量d=Ax，以此作為測(cè)試算例。用不同算法求解得到x′后，計(jì)算前向方均根誤差（root mean square error，RMSE），如式（6）所示。

在申威處理器上測(cè)試了swDCR算法和從核上的追趕法，作為參考，在英特爾處理器Core i7上也測(cè)試了Math Kernel Library 中帶主元選取的LU 分解法（lower-upper triangular decomposition）。

表3展示了不同算法的求解誤差?？梢园l(fā)現(xiàn)，不同算法的誤差受對(duì)角占優(yōu)性、正定性和條件數(shù)的影響基本一致。而swDCR 在對(duì)角占優(yōu)的情況下，誤差和帶主元選取的LU 分解法在同一個(gè)數(shù)量級(jí)。但在非對(duì)角占優(yōu)的情況下有大于10 倍的誤差，這主要是因?yàn)樵贑R算法中，對(duì)角線元素作為除數(shù)出現(xiàn)在前向消去和后向回代的公式中，如果其絕對(duì)值較小則會(huì)產(chǎn)生較大誤差。

Fig.5 Performance improvement of different optimizations of swDCR圖5 swDCR的不同優(yōu)化所產(chǎn)生的性能提升

Table 3 RMSE of different algorithms表3 不同算法的方均根誤差

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

本文通過(guò)分析追趕法、CR算法和PCR算法在申威26010處理器上的性能，設(shè)計(jì)了基于申威處理器的swDCR 算法。通過(guò)聯(lián)合多個(gè)從核上的軟件可控緩存，將過(guò)程中的所有中間變量保存在緩存中，從而減少了主存的訪問(wèn)，提高了算法的效率。通過(guò)利用向量化運(yùn)算等方式，進(jìn)一步優(yōu)化算法的性能，使得該算法在訪存量已經(jīng)達(dá)到理論上最少的前提下能利用硬件理論帶寬的70%。其吞吐率相比主核上的追趕法達(dá)到了單精度43.9 倍和雙精度36.7 倍的加速，相比從核上的追趕法達(dá)到了單精度和雙精度均2.07倍的加速。

swDCR算法可聯(lián)合最多64個(gè)從核使得可求解的問(wèn)題規(guī)模達(dá)到單精度下131 072維和雙精度下65 536維。雖然所能處理的規(guī)模已經(jīng)能滿足絕大多數(shù)現(xiàn)有的實(shí)際應(yīng)用的需求，但是仍然是有限的。在進(jìn)一步工作中，算法的可擴(kuò)展性將要被著重考慮。一方面，對(duì)于未來(lái)可能出現(xiàn)的更大規(guī)模應(yīng)用，還需要對(duì)算法進(jìn)行更多改進(jìn)；另一方面，該算法能否在未來(lái)可能出現(xiàn)的架構(gòu)相似的處理器上獲得足夠的性能，還有待研究。