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

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

0 前言

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

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

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

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

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

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

2 計(jì)算實(shí)現(xiàn)

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

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

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

式（2）計(jì)算過(guò)程如下：

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

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

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

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

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

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

t（n）的計(jì)算代碼：

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

HTt（n）計(jì)算代碼：

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

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

4 計(jì)算效果比較

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

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

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

表1 計(jì)算耗時(shí)比較（單位：s）Tab.1 Computation time comparison（uint：s）

5 結(jié)論

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

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

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

［3］張舒，褚艷利.GPU高性能運(yùn)算之CUDA［M］.北京：中國(guó)水利水電出版社，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］丁洪林.核輻射探測(cè)器［M］.哈爾濱：哈爾濱工程大學(xué)出版社，2010.

［7］裴少英，王南萍.航空γ能譜全譜數(shù)據(jù)的高斯分解方法［J］.科技導(dǎo)報(bà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é)報(bào)，1993（3）：215.

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

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

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

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