許樂，安虹，陳俊仕，張鵬飛，武錚

（中國科學(xué)技術(shù)大學(xué)計算機科學(xué)與技術(shù)學(xué)院，合肥 230026）

0 概述

近幾十年來，計算流體力學(xué)（Computational Fluid Dynamics，CFD）呈現(xiàn)蓬勃發(fā)展的態(tài)勢，在學(xué)術(shù)界與各學(xué)科相互交叉，科研成果日新月異，在工業(yè)界與各領(lǐng)域深度結(jié)合，成為輔助工程設(shè)計的新興技術(shù)。非結(jié)構(gòu)網(wǎng)格是有限元和有限體積等數(shù)值方法中最常用的離散網(wǎng)格，在CFD 計算中具有重要意義，需要基于硬件架構(gòu)精細調(diào)優(yōu)以保證其良好的性能優(yōu)勢。非結(jié)構(gòu)網(wǎng)格計算一般表示為稀疏矩陣運算，其更新時的隨機訪存加劇了數(shù)據(jù)存儲的離散性。隨著網(wǎng)格規(guī)模的增加，離散訪存的規(guī)模也在成倍增加，在帶寬受限的計算系統(tǒng)中成為主要性能瓶頸。離散訪存的計算特點使得并行化時也會出現(xiàn)多個計算任務(wù)同時訪問相同元素的寫后讀沖突和寫寫沖突問題。

神威·太湖之光計算機系統(tǒng)［1］是世界上第一臺峰值性能超過十億億次量級的超算系統(tǒng)，使用完全自主研制的申威26010 異構(gòu)眾核處理器［2］。但非結(jié)構(gòu)網(wǎng)格在眾核處理器上的計算普遍存在離散訪存、數(shù)據(jù)依賴等問題，并行化難度較高。此外，有依賴關(guān)系的算子和對稱矩陣的計算進一步提高了在眾核處理器上實現(xiàn)高性能非結(jié)構(gòu)網(wǎng)格計算的難度。

為了解決非結(jié)構(gòu)網(wǎng)格計算中有依賴關(guān)系算子在眾核處理器上的優(yōu)化問題，本文對大量稀疏網(wǎng)格數(shù)據(jù)進行分析，從網(wǎng)格本身結(jié)構(gòu)和數(shù)據(jù)之間的關(guān)系出發(fā)，提出自適應(yīng)和無依賴的任務(wù)劃分策略，使任務(wù)劃分方法與具體算子不產(chǎn)生綁定關(guān)系，從而提高對不同類型算子的普適性。根據(jù)主從核架構(gòu)的特點，本文提出N 階對角染色算法平衡主從核計算，并在從核計算時摒棄傳統(tǒng)的寄存器通信操作，便于擴展到新一代神威平臺。此外，考慮到計算訪存重疊技術(shù)是申威處理器的常見優(yōu)化策略，本文利用該技術(shù)進一步提升計算效率。

1 研究背景

1.1 申威異構(gòu)眾核處理器

神威·太湖之光計算機系統(tǒng)是我國首臺完全自主研發(fā)的世界第一超算系統(tǒng)，也是我國目前使用最廣泛的高性能計算平臺之一，為經(jīng)濟和社會發(fā)展提供了有力支撐。神威·太湖之光超算系統(tǒng)由高速計算系統(tǒng)和輔助計算系統(tǒng)及配套的互連網(wǎng)絡(luò)和存儲系統(tǒng)組成，配備精準(zhǔn)的資源調(diào)度系統(tǒng)和豐富的并行編程環(huán)境。系統(tǒng)由40 960 塊申威26010 處理器構(gòu)成，內(nèi)存空間為1 024 TB，訪存總帶寬為4 473 TB/s，峰值運算速度達到125PFLOPs，比其他同量級超算系統(tǒng)節(jié)能60%以上。

申威26010 處理器是我國通過自主核心技術(shù)研制的全新異構(gòu)眾核處理器，支持64 位申威指令集。申威處理器由4 個同構(gòu)核組構(gòu)成，每個核組內(nèi)有1 個控制核心（主核）和64 個計算核心（從核），共享統(tǒng)一編址的8 GB 主存，如圖1 所示。

圖1 申威26010 處理器架構(gòu)Fig.1 The architecture of SW26010 processor

在申威26010 處理器中：主核負責(zé)任務(wù)分發(fā)和調(diào)度，工作頻率為1.45 GHz，L1 Cache 為32 KB，L2 Cache（數(shù)據(jù)和指令Cache 混合）為256 KB；從核負責(zé)稠密計算，工作頻率為1.45 GHz，指令Cache 為16 KB。從核采用64 KB 的局部存儲（Local Data Memory，LDM）代替硬件Cache，需要用戶手動完成數(shù)據(jù)的換入換出，有利于充分利用片上存儲空間，但也給編程帶來極大挑戰(zhàn)。

申威處理器98%的計算性能來源于從核陣列，因此，挖掘從核架構(gòu)特性、充分利用從核計算資源十分重要。稀疏矩陣運算的計算強度遠低于稠密矩陣，為達到較高的性能，需要更高訪存帶寬的支持，而申威處理器訪存性能不佳，因此，必須充分利用從核訪存機制來盡可能降低開銷。

從核可以通過gld/gst 指令直接對主存進行訪問，但基準(zhǔn)測試顯示gld/gst 的延遲很高，達到上百節(jié)拍數(shù)，帶寬低于1.5 GB/s，因此，在密集訪存時一般不作考慮。在通常情況下，從核使用DMA 操作連續(xù)訪問主存數(shù)據(jù)可獲得明顯的性能提升。DMA 操作帶寬如表1 所示。

表1 DMA 操作帶寬Table 1 The bandwidth of DMA operation

DMA 操作是指從核LDM 和主存之間的數(shù)據(jù)傳輸，它只能由從核線程發(fā)起，有單從核模式（PE_MODE）、行模式（ROW_MODE）、廣播模式（BCAST_MODE）、廣播行模式（BROW_MODE）和行集合模式（RANK_MODE）5 種，一般情況下常用單從核模式。此外，從核還支持跨步DMA，即按照一定跨步間隔連續(xù)訪問主存。

在每個核組中，64 個從核構(gòu)成8×8 的從核陣列，從核間的數(shù)據(jù)交換通過寄存器通信機制（Register Level Communication，RLC）進行，該機制以生產(chǎn)者/消費者模式運行，各從核以1 個向量長度為單位在其行或列上進行數(shù)據(jù)廣播和接收。如圖2 所示，源從核首先將256 位對齊的數(shù)據(jù)加載到寄存器中，然后通過發(fā)送緩沖區(qū)（Send Buffer）將它們發(fā)送到從核通信網(wǎng)格；目的從核通過接收緩沖區(qū)（Receive Buffer）從通信網(wǎng)格中獲取數(shù)據(jù)，并將其加載到本地寄存器。

圖2 寄存器通信機制原理Fig.2 RLC principle

寄存器的通信延遲通常低至幾個周期，如表2所示，這使得從核間的數(shù)據(jù)可以快速交換，但每個從核只能通過行廣播向同一行中的一個或多個從核發(fā)送數(shù)據(jù)，或通過列廣播在同一列中發(fā)送數(shù)據(jù)，這給數(shù)據(jù)的自由傳輸帶來極大限制，不利于開發(fā)有復(fù)雜依賴關(guān)系的從核程序。

表2 寄存器通信延遲Table 2 Register communication latency

1.2 相關(guān)研究

近年來，CFD 數(shù)值模擬軟件作為高性能計算領(lǐng)域中的重要應(yīng)用，已經(jīng)廣泛部署于眾多超算平臺上。商業(yè)軟件因其平臺適配性強和穩(wěn)定性高，曾是非結(jié)構(gòu)網(wǎng)格計算軟件的主流，如Fluent［3］、UMS3D［4-5］、FUN3D［6-7］、TAU［8］、CFD++［9］、NSU3D［10］等，但其內(nèi)部源代碼不對外公開，很難精準(zhǔn)解決用戶需求。相比而言，開源CFD軟件可以滿足用戶不同的開發(fā)需求，如OpenFOAM（Open source Field Operation And Manipulation）［11］、Featflow［12］、Gerris［13］、Code_Saturne［14］等，其中，OpenFOAM是目前應(yīng)用范圍最廣、可擴展性最強、解算器最全的開源軟件包。

CFD 軟件中的核心部分是高精度數(shù)值求解器，越來越多的研究人員開始使用異構(gòu)加速方法來改進線性方程組的求解。BOLZ等［15］首次在GPU 上實現(xiàn)了高效的SpMV 算子，此后，基于ELLPACK格式，BELL等［16］設(shè)計HYB存儲格式實現(xiàn)了SpMV，而VáZQUEZ等［17］、MONAKOV等［18］和CHOI等［19］分別設(shè)計了ELLPACK-R、sliced-ELLPACK 和blocked ELLPACK 存儲格式。

基于CSR格式，KOZA等［20］提出了CSMR格式以提高數(shù)據(jù)重用，GREATHOUSE等［21］和ASHARI等［22］分別提出CSR-adaptive和ACSR格式以解決負載均衡問題。此外，MERRILL等［23］和LIU等［24］還分別提出MCSR和CSR5格式，取得了良好的性能提升。對于分塊問題，BULU?等［25］、ASHARI等［26］、LIANG等［27］和YAN等［28］分別提出了CSB、BRC、HCC和BCCOO存儲格式。

在國產(chǎn)申威處理器上，文獻［29］基于CSR 格式提出動靜態(tài)優(yōu)化方法，其相比主核實現(xiàn)取得了6 倍的加速效果，文獻［30］進一步提出雙邊多級劃分方法，相比主核實現(xiàn)取得了12 倍以上的加速效果。文獻［31］基于非結(jié)構(gòu)網(wǎng)格實現(xiàn)了稀疏下三角方程求解器，文獻［32］提出基于排序思想的通用眾核優(yōu)化算法以減少非結(jié)構(gòu)網(wǎng)格計算中的隨機訪存，隨后，文獻［33］提出兩階段優(yōu)化方法以克服大規(guī)模不規(guī)則訪存和帶寬利用率低的問題。

2 非結(jié)構(gòu)網(wǎng)格計算

OpenFOAM 是一款對連續(xù)介質(zhì)力學(xué)問題進行數(shù)值計算的開源C++類庫，因其模塊化和可定制程度高，目前已成為超算平臺上主流的CFD 軟件。OpenFOAM基于C++語言開發(fā)，利用操作符重載、繼承和模板等面向?qū)ο筇匦?，支持?shù)據(jù)預(yù)處理、數(shù)據(jù)后處理和自定義偏微分方程求解，框架內(nèi)提供網(wǎng)格生成、有限體積法、線性方程組求解、輸入輸出處理等功能。

如圖3 所示，OpenFOAM 中非結(jié)構(gòu)網(wǎng)格一般通過鄰接矩陣表示，而非結(jié)構(gòu)網(wǎng)格的稀疏性又使得鄰接矩陣為稀疏矩陣。稀疏運算計算強度低，在眾核處理器上仍有很大的優(yōu)化空間，因此，本文基于申威處理器提出非結(jié)構(gòu)網(wǎng)格計算的通用加速框架。

圖3 非結(jié)構(gòu)網(wǎng)格及其矩陣表示Fig.3 The unstructured gird and its matrix representation

2.1 基礎(chǔ)數(shù)據(jù)結(jié)構(gòu)

OpenFOAM 中稀疏矩陣按照LDU 格式存儲，矩陣上除對角線元素外各非零元素用一個三元組（行號，列號，數(shù)值）表示，如圖4 所示。

圖4 LDU 存儲格式Fig.4 LDU storage format

整個稀疏矩陣分為對角部分（Diag）、上三角部分（Upper）和下三角部分（Lower），對角線元素即網(wǎng)格頂點數(shù)據(jù)，上三角和下三角元素為網(wǎng)格邊數(shù)據(jù)。上三角元素的行索引（Row）對應(yīng)下三角元素的列索引（Col），而上三角元素的列索引對應(yīng)下三角元素的行索引，數(shù)據(jù)可以使用相同的索引數(shù)組存儲。在從核并行計算時，雖然矩陣上三角每行元素按序排列，但是下三角每行元素并不連續(xù)，訪問下三角元素時很難滿足空間局部性。

2.2 計算特征分析

非結(jié)構(gòu)網(wǎng)格能夠適應(yīng)各種復(fù)雜結(jié)構(gòu)的面網(wǎng)格劃分，是流體力學(xué)仿真軟件中最主要的空間離散化解決方案。非結(jié)構(gòu)網(wǎng)格計算存在3 種模式，即頂點狀態(tài)更新（點更新）、邊狀態(tài)更新（邊更新）和通過鄰居頂點與鄰接邊更新頂點狀態(tài)（邊點更新）。

如圖5 所示，邊更新計算特征為S(e)=f(S(e)，S'(v1)，S'(v2))，典型計算函數(shù)為S(e)+=∑(S'(vi))。

圖5 邊更新模式Fig.5 The edge update mode

如圖6 所示，點更新計算特征為S(v)=f(S(v)，S'(e1)，S'(e2)，S'(e3)，S'(e4)，S'(e5))，典型計算函數(shù)為S(v)+=∑(S'(ei))。

圖6 點更新模式Fig.6 The point update mode

如圖7 所示，邊點更新計算特征為S(v)=f(S(v)，S'(e1)，S'(e2)，S'(e3)，S'(e4)，S'(e5)，S''(v1)，S''(v2)，S''(v3)，S''(v4)，S''(v5))，典型計算函數(shù)為S(v)+=∑(S'(ei) ·S''(vi))。

圖7 邊點更新模式Fig.7 The edge and point update mode

3 種模式的共同特征在于狀態(tài)信息在頂點與邊之間流動。本文將邊視為基本單元，頂點視為連接單元，計算過程則為遍歷非結(jié)構(gòu)網(wǎng)格的所有邊，獲取每條邊上左右頂點狀態(tài)，并與邊自身狀態(tài)進行運算，最終用運算結(jié)果更新左右頂點狀態(tài)或邊自身狀態(tài)，如表3 所示。

表3 非結(jié)構(gòu)網(wǎng)格計算模式Table 3 Unstructured grid computing mode

由于數(shù)據(jù)關(guān)系的依賴性和對數(shù)據(jù)訪問的離散性等固有特點，導(dǎo)致非結(jié)構(gòu)網(wǎng)格計算具有局部性差、數(shù)據(jù)相關(guān)、離散訪存、并發(fā)度低等問題，在眾核處理器上進行優(yōu)化時難度較大，性能偏低。

3 基于申威處理器的眾核加速方法

3.1 算法思想

由于非結(jié)構(gòu)網(wǎng)格的稀疏性，導(dǎo)致算法在計算時對元素訪問并不連續(xù)，無法充分利用訪存帶寬。此外，非結(jié)構(gòu)網(wǎng)格的離散寬度（稀疏矩陣中非零元素與該行對角線元素之間的距離）較大，造成訪存間隔過大，難以滿足空間局部性。如圖8 所示，1 號邊和2 號邊的距離較大，遍歷時雖然行索引相同，但是列索引相距較遠，不滿足空間局部性，訪存性能較差。因此，本文采用分塊劃分的思想，將一段時間內(nèi)的訪存數(shù)據(jù)盡可能集中存儲在較快的存儲器上，降低從下層存儲器讀取數(shù)據(jù)的時間開銷。

圖8 非結(jié)構(gòu)網(wǎng)格的鄰接矩陣Fig.8 The adjacent matrix of unstructured grid

此外，非結(jié)構(gòu)網(wǎng)格計算存在大量的對稱矩陣，為節(jié)省存儲空間，一般僅保留上三角矩陣，但是在計算時需要對稱更新，因此，在眾核處理器上的并行化存在數(shù)據(jù)相關(guān)和輸出相關(guān)（寫沖突）的問題。例如，在圖8 中，1 號邊和2 號邊位于同一行，表示其對相同目標(biāo)結(jié)果進行更新，存在依賴關(guān)系，如果計算任務(wù)被分配至2 個不同的線程執(zhí)行，可能會發(fā)生寫后讀沖突。此外，1 號邊和6 號邊位于同一列，且1′號邊和6 號邊位于同一行，由于對稱更新，則1 號邊和1′號邊同時更新時如果其他計算任務(wù)對6 號邊更新，就會發(fā)生寫寫沖突。本文在分析非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)特點和計算特征后，提出并行度更高的無依賴任務(wù)劃分方法，將數(shù)據(jù)相關(guān)和輸出相關(guān)的計算分配到相同任務(wù)隊列。

本文提出一種N 階對角染色算法，非結(jié)構(gòu)網(wǎng)格邊線沿對角方向劃分為大小相同的方塊后，將有依賴關(guān)系的方塊染上同種顏色，分配到同一任務(wù)隊列中進行并行計算。然后，染色器不斷向外擴展并對其他類對角方塊染色。算法根據(jù)方塊內(nèi)元素密度決定染色階數(shù)，即向外擴展對角線的次數(shù)。該算法支持非結(jié)構(gòu)網(wǎng)格下大多數(shù)的算子模型，特別是有依賴關(guān)系的算子。算法執(zhí)行步驟如下：

1）網(wǎng)格預(yù)處理及自適應(yīng)劃分。獲得當(dāng)前頂點數(shù)和邊數(shù)，記錄邊線所連頂點。根據(jù)LDM 存儲空間等，針對不同網(wǎng)格拓撲自適應(yīng)確定分塊大?。ㄟ吘€所連頂點范圍）和染色階數(shù)，保證計算單元負載均衡。

2）分塊染色及重排整理。根據(jù)分塊大小，從對角塊向外對網(wǎng)格逐階染色，按照邊隨頂點走、一階一類色的原則，同時建立索引表記錄頂點的塊內(nèi)位置與全局位置關(guān)系，方便后續(xù)計算結(jié)果更新。

3）任務(wù)調(diào)度。同色塊分配至同一任務(wù)隊列，采用動態(tài)調(diào)度方法管理任務(wù)隊列以維持從核負載平衡。

4）訪存及計算。從核通過DMA 操作完成網(wǎng)格邊線重排序，將當(dāng)前隊列內(nèi)染色塊加載至LDM 中并執(zhí)行計算。同時，為了隱藏DMA 操作時間，從核在進行當(dāng)前計算時開始下一輪DMA 操作，使得計算與訪存重疊。在從核計算過程中，主核同時負責(zé)未染色塊計算，因為未染色塊更為稀疏，局部性更差，更適合通過主核計算，而從核通過DMA 對數(shù)據(jù)的換入換出往往會帶來更高的時間開銷，不利于發(fā)揮其性能優(yōu)勢。

3.2 算例分析

考慮到計算的常見性和代表性，本文以典型算子稀疏矩陣向量乘（Sparse Matrix Vector Multiplication，SpMV）為例分析眾核加速方法。

作為最常見的稀疏運算，雙精度SpMV 的計算強度僅為0.080～0.125 FLOPs/Byte，在帶寬受限的眾核處理器上性能較差。SpMV 算子描述如算法1 所示。

算法1SpMV 算子

在算法1 中：V為輸入/輸出向量值，即網(wǎng)格頂點狀態(tài)；E.row/E.col 和E.val 分別為稀疏矩陣的行/列索引和值，即網(wǎng)格邊狀態(tài)。頂點狀態(tài)的更新由相連邊狀態(tài)及其連接頂點狀態(tài)的乘積累加得到。本文所提算法執(zhí)行過程如下：

1）在算法2 中，設(shè)最小并行單位為Δ，即分塊最小頂點范圍。根據(jù)LDM 空間大小、DMA 特性和讀取的網(wǎng)格拓撲信息，自適應(yīng)確定分塊因子大小α。

算法2分塊因子判決

2）在算法3 中，根據(jù)分塊大小αΔ掃描邊線，先將主對角塊染色并建立邊索引表，例如第一塊頂點范圍在0～（αΔ-1）內(nèi)，第二塊頂點范圍在αΔ～（2αΔ-1）內(nèi)，同時將頂點全局位置轉(zhuǎn)換為塊內(nèi)位置，建立關(guān)系映射表，原因是塊內(nèi)計算時不能使用全局地址。按照同樣的方法向二階及以上階擴展，皆為雙側(cè)次對角塊，即包括上三角和對稱的下三角兩部分。隨后，算法4 從對角塊中挑選較稠密對角塊進行染色，挑選標(biāo)準(zhǔn)由塊內(nèi)節(jié)點密度和網(wǎng)格整體密度決定，根據(jù)大量測試后得出。未被挑選的非稠密塊則分配給主核任務(wù)隊列，依據(jù)主從核的不同特點實現(xiàn)任務(wù)分配。三階對角染色示意圖如圖9 所示。

圖9 三階對角染色Fig.9 Third-order diagonal dyeing

算法3分塊染色及重排整理

算法4挑選染色對角

3）建立從核任務(wù)隊列queue，將全部一階色塊分入隊列，從核從隊列中獲取任務(wù)并完成計算。類似地，在一階色塊完成計算后，其他同階色塊依次被分配到任務(wù)隊列，隊列內(nèi)從核運行狀態(tài)一致，從而避免寫后讀沖突和寫寫沖突。

4）在算法5 中，從核通過DMA 獲取塊內(nèi)頂點狀態(tài)和邊狀態(tài)以及索引表并完成更新，在計算時可以預(yù)取下一輪數(shù)據(jù)，從而使得DMA 時間被有效隱藏，如圖10 所示。在從核計算的同時，主核負責(zé)其他未染色塊的計算，從而實現(xiàn)主從核異步并行，進一步提升計算效率。

圖10 計算訪存異步重疊Fig.10 The asynchronous overlap of computation and memory access

算法5訪存及計算

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

本次實驗基于申威26010 眾核處理器，硬件參數(shù)如表4 所示，采用swg++編譯器編譯全部C/C++程序。

表4 申威26010 處理器硬件參數(shù)Table 4 The hardware parameters of SW26010 processor

4.1 不同網(wǎng)格的性能分析

為保證算法性能的可靠性，本文選取典型稀疏算子SpMV 作為標(biāo)準(zhǔn)測試算子，隨機輸入非結(jié)構(gòu)網(wǎng)格實例進行測試分析。圖11 和圖12 分別為SpMV算子在不同網(wǎng)格規(guī)模下加速算法與主核樸素算法的運行時間及加速比，以驗證加速算法的通用性和優(yōu)化效果。

圖11 不同網(wǎng)格規(guī)模下的優(yōu)化性能Fig.11 Optimization performance under different grid scales

圖12 不同網(wǎng)格規(guī)模下的加速效果Fig.12 Acceleration effect under different grid scales

從圖11 和圖12 中可以看出，隨著網(wǎng)格規(guī)模的增加，加速算法的加速效果基本保持穩(wěn)定，因為網(wǎng)格劃分根據(jù)輸入網(wǎng)格密度和拓撲自適應(yīng)調(diào)整，染色階數(shù)也根據(jù)當(dāng)前對角線密度自動判決，因此加速算法能在多數(shù)網(wǎng)格規(guī)模下保持穩(wěn)定的性能優(yōu)勢，通用性較強。相比于主核上運行的SpMV 算子，組合加速算法獲得了平均10 倍左右的加速比，最高加速比可達24 倍。

4.2 不同算子的性能分析

本文設(shè)計非結(jié)構(gòu)網(wǎng)格計算在申威眾核處理器上的通用加速方法，因此，需要選取多種算子進行綜合分析。SpMV 算子的加速比如圖12 所示，Integration算子和calcLudsFcc 算子的加速比分別如圖13 和圖14 所示。

圖13 Integration 算子在不同網(wǎng)格規(guī)模下的加速效果Fig.13 Acceleration effect of Integration operator under different grid sizes

圖14 calcLudsFcc 算子在不同網(wǎng)格規(guī)模下的加速效果Fig.14 Acceleration effect of calcLudsFcc operator under different grid sizes

經(jīng)過對上述2 種算子在不同網(wǎng)格規(guī)模下的測試發(fā)現(xiàn)，組合加速算法分別獲得了10.22 倍和6.82 倍的平均加速比，而且本文算法對不同算子模型的性能表現(xiàn)差異并不明顯，在有依賴和無依賴情況下都能取得穩(wěn)定的性能優(yōu)勢，說明算法在任務(wù)劃分和數(shù)據(jù)映射上并沒有以犧牲性能為代價，自適應(yīng)和無依賴任務(wù)劃分方法取得了良好效果。由于算子在從核的計算和訪存是基于染色后的數(shù)據(jù)塊，其加速效果與數(shù)據(jù)塊中的數(shù)據(jù)分布有關(guān)，在數(shù)據(jù)集中度較高的網(wǎng)格實例中，算子能獲得顯著的性能提升，可達20 多倍，但在數(shù)據(jù)非常離散的情況下效果一般。

4.3 不同優(yōu)化策略的性能分析

為了說明N 階對角染色算法和自適應(yīng)任務(wù)劃分方法的有效性，以SpMV 算子為例，分別采用非N 階對角染色的分塊算法和固定分塊大小的任務(wù)劃分方法進行對比實驗，并以主核樸素算法為基準(zhǔn)，實驗結(jié)果如圖15 所示。

圖15 不同優(yōu)化方法的加速效果Fig.15 Acceleration effect of different optimization methods

N 階對角染色算法通過分析對角塊密度來選擇是否染色當(dāng)前對角塊，而普通分塊算法缺少對角塊信息，易將過于稀疏的對角塊交由從核陣列計算。將全部矩陣塊映射到從核陣列會造成極大的性能損失，本文通過固定前100 階對角塊由從核計算來模擬非染色的普通分塊算法性能。自適應(yīng)劃分方法根據(jù)LDM 空間大小、DMA 特性和網(wǎng)格拓撲信息確定分塊大小，可以充分利用LDM 空間，而固定分塊大小則容易造成對LDM 空間的利用不足。實驗結(jié)果表明，非N 階對角染色的分塊算法平均加速比為2.64 倍，固定分塊大小的任務(wù)劃分方法平均加速比僅為1.8 倍，難以發(fā)揮眾核架構(gòu)的計算能力，甚至有負優(yōu)化效果出現(xiàn)。采用自適應(yīng)任務(wù)劃分的N 階對角染色算法能有效利用LDM 空間并根據(jù)塊內(nèi)密度選擇從核或主核來執(zhí)行計算，平均加速比可達10 倍。

5 結(jié)束語

為提升非結(jié)構(gòu)網(wǎng)格計算中有依賴關(guān)系算子在眾核處理器上的性能，本文針對非結(jié)構(gòu)網(wǎng)格的計算特點，提出一種眾核處理器上的通用加速方法，并基于申威26010 處理器架構(gòu)對其進行精細調(diào)優(yōu)。通過自適應(yīng)任務(wù)劃分方法將從核離散訪存組織為批量訪存，以降低訪存開銷。采用無依賴劃分策略避免計算時的數(shù)據(jù)沖突，通過N 階對角染色算法將計算任務(wù)分類調(diào)度執(zhí)行，從而有效利用主從核的架構(gòu)特點。此外，采用計算訪存重疊技術(shù)進一步提升計算性能。實驗結(jié)果表明，該方法在不同網(wǎng)格規(guī)模和不同算子模型下都取得了良好的加速效果，在有依賴和無依賴情況下均能保持穩(wěn)定的性能優(yōu)勢，證明了任務(wù)劃分方法的有效性。但是，本文所提方法對于數(shù)據(jù)分布極為分散的非結(jié)構(gòu)網(wǎng)格仍存在一定局限性，下一步將結(jié)合排序算法對網(wǎng)格數(shù)據(jù)進行重排整理，提升數(shù)據(jù)的局部性，增加在從核陣列計算的數(shù)據(jù)塊，從而滿足更多稀疏網(wǎng)格數(shù)據(jù)的眾核計算需求。