基于OpenCL的顯卡加速射線能譜數(shù)據(jù)處理

劉 端，葛良全，張慶賢，谷 懿

（成都理工大學(xué)，成都 610059）

0 前言

現(xiàn)代儀器譜解析，由原來的特征峰解析趨向于全譜解析，多采用矩陣運算實現(xiàn)。因此計算機的數(shù)據(jù)處理能力和計算時間，是儀器譜解析的重要指標。為了提高儀器譜解析的效率，一般需要對算法進行精簡。但是在數(shù)學(xué)基礎(chǔ)上，精簡量非常的有限，以矩陣運算的復(fù)雜度為例，其精簡程度有限［5］而且實際在計算機上的執(zhí)行效率也不高。另外一種方法，就是提高計算機的計算能力。目前計算機多采用雙核或者四核CPU，其數(shù)據(jù)計算能力大大加強。同時在系統(tǒng)中擁有用于圖像處理的GPU，其具有并行計算的特征［3］。但是目前進行核數(shù)據(jù)處理軟件中，多數(shù)仍然采用串行單線程為主，沒有完全發(fā)揮多數(shù)據(jù)處理器計算機的計算能力。

在計算能力不斷提高的計算機系統(tǒng)中，除了傳統(tǒng)的計算設(shè)備CPU外，其它的設(shè)備如GPU、DSP也逐漸進入通用計算行列。由于GPU同CPU相比較，具有其獨特的數(shù)據(jù)并行計算能力，GPU目前已經(jīng)在科學(xué)計算領(lǐng)域扮演著越來越重要的角色。但不同廠家不同型號的設(shè)備的差異性，給開發(fā)增加了不少難度，由此出現(xiàn)了一些GPGPU開發(fā)環(huán)境，如OpenCL、CDUA和DirectCompute等，相比于其它兩個OpenCL支持設(shè)備最為廣泛。OpenCL（Open Computing Language，開放運算語言）是一個面向異構(gòu)系統(tǒng)通用目的并行編程的開放式的免費標準，它規(guī)定了一套標準的應(yīng)用程序接口API，可適用于不同設(shè)備［4］。目前常見的如Intel公司的多核CPU、NVIDIA公司的顯卡、AMD公司的多核CPU和顯卡都可支持［4］。

極大似然估計的儀器譜解析方法是在文獻［1］中使用的。極大似然估計方法用來對NaI（Tl）探測的譜線進行解析，通過對合成譜線的解析，說明該方法可以有效地分解能量差為2／3FHWM的重疊峰。但該方法計算量較大，需要經(jīng)過5 000次以上的迭代計算才能得到結(jié)果，對于1 024道的譜線解析將耗時幾十分鐘，從而不能達到實際應(yīng)用的目的。作者在本文使用OpenCL在NVIDIA公司的Quadro FX 2700M顯卡，以開發(fā)并行化的極大似然估計算法為例，提高算法的運行效率，減少算法的執(zhí)行時間，介紹并行計算在核數(shù)據(jù)處理中的應(yīng)用。

1 OpenCL的計算模型構(gòu)架

OpenCL是面向異構(gòu)平臺的編程標準，將各種計算設(shè)備（如CPU、GPU、DSP）組織為統(tǒng)一的平臺，并提供完整的并行編程框架。OpenCL程序運行分為Host端與計算設(shè)備，Host端程序主要負責管理協(xié)調(diào)各個不同的計算設(shè)備；計算設(shè)備上運行的程序稱為kernel。通常多個相似的kernel組成一個工作空間，由同一個計算設(shè)備并行執(zhí)行，其中每個kernel稱為work－item。每個工作空間由NDRange進行索引（如圖1所示），work－item又被分給不同work－group進行管理。同一個work－group中的work－item都由同一個計算單元執(zhí)行。在同一個work－group中，可以使用其內(nèi)部高速的loacl memory進行數(shù)據(jù)共享和線程同步操作，但是不同work－group中的work－item，只能通過速度較慢的global memory共享數(shù)據(jù)并且不能進行線程同步。

作者在本文使用的是NVIDIA Quadro FX 2700M顯卡，其主要計算構(gòu)架如圖2所示，在圖2中的SM（Stream Multiprocessor流多處理器），即前文提到的計算單元。該顯卡有六個SM，每個SM擁有八個SP（Stream Processor流處理器。OpenCL的work－item運行在SP上，而同一work－group中的work－item運行于同一SM。FX 2700M顯卡共擁有四十八個SP（即通常所說的48核），其核心頻率為1 350MHz，理論浮點計算能力為190.8GFLOPS（十億次浮點運算每秒）。

2 計算實現(xiàn)

極大依然法［1］的計算主要是對一方程組采用Richardson－Lucky（R－L）迭代公式進行求解，表達式為：

該計算主要是求式（1）右邊的與x相乘的式子的值：

則t（n）（i）＝g（i）／u（n）（i）。

式（2）計算過程如下：

步驟①先求u（n）＝Hx（n）；步驟②再由t（n）（i）＝g（i）／u（n）（i）求t（n）；步驟③再求 HTt（n）得到結(jié)果。

此過程中有兩次矩陣與向量的乘法運算，此運算的實現(xiàn)是關(guān)鍵。

3 OpenCL內(nèi)核程序的實現(xiàn)

矩陣與向量乘法式（3）對于每行值的計算，具有計算過程一致、結(jié)果相互不影響的特點，因此可以將其并行執(zhí)行。1 024道的數(shù)據(jù)可以分給1 024個線程分別進行計算。

FX 2700M顯卡中單個SM雖然擁有八個SP，但對數(shù)據(jù)的并行訪問是十六個work－item一起同時進行。基于對內(nèi)存訪問速度的考慮，最佳的訪問方式為每個work－group所含的work－item為十六的整數(shù)倍。由于每道的計算都要用到t（n），對于同一work－group內(nèi)的計算，可將t（n）存入快速的local memory，讓同一 work－group內(nèi)的work－item共享訪問以減少數(shù)據(jù)讀入次數(shù)。因此應(yīng)提高每個work－group的work－item擁有量，但是work－group數(shù)量減少會影響計算單元的使用效率。作者在本文將1 024個work－item分為十六個work－group，每個work－group含有六十四個work－item。

OpenCL內(nèi)核kernel程序可以用OpenCL C語言編寫。

t（n）的計算代碼：

由于要合并訪問global memory中數(shù)據(jù)需要把矩陣H按列存儲，可共享的向量x讀入快速local memory。

HTt（n）計算代碼：

同理，此次需要計算的矩陣為HT需要將其按列存儲，即H按行存儲。由于H擁有按列存儲和按行存儲兩種模式，所以在計算之前將其分別存儲為兩個矩陣。

進一步觀察H的值發(fā)現(xiàn)，其每列數(shù)據(jù)都是隨著離矩陣對角線距離的增加而逐漸減小，在距離大于一定值過后可以忽略不計。矩陣數(shù)據(jù)集中在矩陣對角線兩側(cè)，可將H當作帶狀矩陣處理，以降低計算量。

4 計算效果比較

本測試機器配置為CPU Intel Core 2T9600 2.8GHz，其理論浮點計算能力為22.4GFLOPS，內(nèi)存4G，顯卡為Quadro FX 2700M，顯存512M，理論浮點計算能力為190.8GFLOPS。處理1024道譜線耗時見表1。

運用該程序?qū)σ粌x器譜就行處理。譜線總道數(shù)為1 024，其中314至634道有三處重峰。見下頁，圖3至圖6是對重峰的分離效果。

由圖3至圖6可知，1號重峰500次便可分離，2號重峰1 000次就可分離，但是3號要100 000次才可分離，若用CPU程序要耗時近20min，而GPU只用了不到7s。

表1 計算耗時比較（單位：s）Tab.1 Computation time comparison（uint：s）

5 結(jié)論

作者在本文通過OpenCL利用GPU的強大的并行計算能力和優(yōu)化算法，降低了儀器譜全譜解析所耗時間，同單純的采用CPU計算相比，能夠有效地縮短計算時間。對同一譜線測試表明，采用GPU并行計算，其計算時間減少為采用單純CPU計算的1%。GPU并行計算，對核能譜進行全譜解析具有其自身的優(yōu)勢，增加了該方法的實際可行性。

［1］張慶賢，葛良全，谷懿，等.基于極大似然估計的NaI（Tl）晶體儀器譜解析方法研究［J］.核技術(shù)，2011（8）：11.

［2］張慶賢.航空γ能譜特征和儀器譜解析方法研究［D］.成都：成都理工大學(xué)，2010.

［3］張舒，褚艷利.GPU高性能運算之CUDA［M］.北京：中國水利水電出版社，2009.

［4］AAFTAB MUNSHI.The OpenCL Specification［M］.Khronos OpenCL Working Group，2011.

［5］COPPERSMITH D，WINOGRAD S.Matrix multiplication via arithmetic progressions［M］.Symbolic Compute，1990.

［6］丁洪林.核輻射探測器［M］.哈爾濱：哈爾濱工程大學(xué)出版社，2010.

［7］裴少英，王南萍.航空γ能譜全譜數(shù)據(jù)的高斯分解方法［J］.科技導(dǎo)報，2006（11）：11.

［8］IAEA.Airborne gamma ray spectrometer surveying［R］.IAEA Technical Reports，2003.

［9］ICRU.International Commission on Radiological U－nits and Measurements［R］.ICRU report 53，Bethesda，Maryland：U.S.A，1994.

［10］DUBOIS G，DE CORT M.Mapping 137Cs deposition：data validation methods and data interpretation［J］.Journal of Environmental Radioactivity，2001，53（3）：271.

［11］李惕碚，吳枚.高能天文中成像和解譜的直接方法［J］.天體物理學(xué)報，1993（3）：215.

［12］李惕碚.能不能用低分辨儀器實現(xiàn)高分辨觀測［J］.大自然探索，1999（1）：45.

［13］龐巨豐.γ能譜數(shù)據(jù)分析［M］.西安：陜西科學(xué)技術(shù)出版社，1993.

［14］王振國.基于極大似然原理的天然圖像盲恢復(fù)與重建［D］.鄭州：解放軍信息工程大學(xué)，2006.

［15］張賢達.矩陣分析與應(yīng)用［M］.北京：清華大學(xué)出版社，2004.