田 冪,胡 亮,車喜龍

(吉林大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春 130012)

粒子群優(yōu)化算法(PSO)具有穩(wěn)健、 有效、 簡(jiǎn)易和適應(yīng)性廣的特點(diǎn),應(yīng)用廣泛[1-4].標(biāo)準(zhǔn)粒子群優(yōu)化算法(SPSO)[5]是對(duì)原始算法的直接擴(kuò)展.盡管SPSO算法擁有許多優(yōu)點(diǎn),但在大規(guī)模問題,如在大維度和大群數(shù)上,仍需較長(zhǎng)時(shí)間尋找解決方案.SPSO算法的結(jié)構(gòu)存在內(nèi)在的并行性,通過(guò)直接刪除更新之間的依賴,編寫同步PSO,同時(shí)允許粒子不嚴(yán)格依賴最新信息[6],放寬同步限制,可在并行計(jì)算架構(gòu)下實(shí)現(xiàn)SPSO算法.

圖形處理單元(GPU)因其并行計(jì)算機(jī)制和快速的浮點(diǎn)操作能力,在科學(xué)計(jì)算上展現(xiàn)出較強(qiáng)的優(yōu)勢(shì),因此應(yīng)用廣泛.與其他加速平臺(tái)相比,如CellBE和FPGAs,GPU已在HPC研究領(lǐng)域占主導(dǎo)地位[7],同時(shí)使用CUDA開發(fā)工具[8]可輕松地進(jìn)行編程.

利用可編程的GPU加速SPSO算法的可行性與有效性已被證實(shí)[9-10],然而這些研究都沒有利用最新的硬件版本[11].Fermi 同樣是雙架構(gòu)設(shè)計(jì),即CUDA計(jì)算架構(gòu)和圖形架構(gòu),兩種架構(gòu)可靈活切換,但Fermi幾乎是重新設(shè)計(jì)并更注重通用計(jì)算架構(gòu).本文利用Fermi GPU的同時(shí),在最新架構(gòu)下實(shí)現(xiàn)了并行版本的SPSO算法.

1 SPSO算法與Fermi架構(gòu)

1.1 SPSO算法

圖1 SPSO算法流程Fig.1 Flowchart of SPSO

SPSO系統(tǒng)初始化是一群隨機(jī)的粒子,并通過(guò)更新粒子尋找多維解決空間的最優(yōu)值.每個(gè)粒子依據(jù)其自己發(fā)現(xiàn)的本地最近解決方案進(jìn)行移動(dòng).相鄰最佳解決方案由與其相鄰的左鄰居和右鄰居發(fā)現(xiàn).每次迭代中粒子更新其速度和位置直到滿足終止條件.SPSO算法的基本流程如圖1所示.

1.2 Fermi架構(gòu)

程序可分為兩部分：主機(jī)端部分和設(shè)備端部分.主機(jī)端部分在CPU上執(zhí)行,管理數(shù)據(jù)的傳輸,上下文的預(yù)處理和后處理,調(diào)度設(shè)備端部分的執(zhí)行； 設(shè)備端部分在GPU上執(zhí)行,由一個(gè)或多個(gè)稱為kernel的函數(shù)組成， 每個(gè)kernel是程序的一個(gè)并行單元,由許多并行的線程執(zhí)行.Fermi架構(gòu)支持不同kernel的同時(shí)執(zhí)行,它的流水技術(shù)可減少應(yīng)用上下文切換帶來(lái)的開銷,同時(shí)支持kernel間的快速通信.

線程是設(shè)備端部分程序的最小執(zhí)行單位,以grid和block的形式組成.一個(gè)block由許多線程組成,同一block中的線程共享數(shù)據(jù)并同步它們的執(zhí)行.同一kernel的一些block組成一個(gè)grid.線程以warp的形式進(jìn)行調(diào)度和執(zhí)行,同一warp內(nèi)的所有線程可認(rèn)為是同時(shí)運(yùn)行的,存儲(chǔ)訪問也以warp形式工作.

Fermi GPU以圖形處理集群可擴(kuò)展陣列(GPC)為基礎(chǔ),每個(gè)GPC由4個(gè)流水多重處理器(SM)組成.每個(gè)SM有32個(gè)CUDA核,16個(gè)Load/Store(LD/ST)單元和4個(gè)特殊函數(shù)單元(SFU).CUDA核由全流水整數(shù)算術(shù)邏輯單元(ALU)和浮點(diǎn)單元(FPU)組成.一個(gè)GPU執(zhí)行一個(gè)或多個(gè)grid； 一個(gè)SM執(zhí)行一個(gè)或多個(gè)block； 一個(gè)ALU/FPU/SFU執(zhí)行一個(gè)線程.

在GPU芯片外,有64位的GDDR5 DRAM,可視為是global memory，local memory，constant memory 和texture memory的組合,由主機(jī)端緩存和讀寫.Fermi對(duì)local memory，shared memory和global memory使用統(tǒng)一的地址空間完成加載和存儲(chǔ)操作.在GPU芯片內(nèi),有3種類型的快速存儲(chǔ)器：register files，shared memory和cache,它們由設(shè)備端進(jìn)行讀寫操作.

2 SPSO算法在Fermi架構(gòu)下的加速

SPSO算法在GPU上的執(zhí)行過(guò)程如下.

算法1

1) 在global memory中申請(qǐng)所需的變量空間；

2) 將數(shù)據(jù)由主機(jī)內(nèi)存拷貝到global memory內(nèi)相應(yīng)的變量中；

fori=1 to Iter do //并行執(zhí)行for語(yǔ)句中步驟

3) 計(jì)算每個(gè)粒子的適應(yīng)值；

4) 更新每個(gè)粒子的最佳適應(yīng)值和最佳位置；

5) 更新粒子群的最佳適應(yīng)值和最佳位置；

6) 更新每個(gè)粒子的速度和位置；

end for

7) 將結(jié)果數(shù)據(jù)由global memory拷貝回主機(jī)內(nèi)存.

本文使用的主要參數(shù)如下：Iter表示SPSO算法的迭代次數(shù)；D表示每個(gè)粒子的維度；N表示粒子的數(shù)目.

2.1 加速策略

在Fermi上加速SPSO算法,需要考慮計(jì)算量與block維度的選擇、 存儲(chǔ)的使用和控制指令流的使用.

2.1.1 計(jì)算量與block維度的選擇 如果計(jì)算規(guī)模很小,使用GPU進(jìn)行計(jì)算將不能彌補(bǔ)傳輸數(shù)據(jù)帶來(lái)的延時(shí),即使用GPU計(jì)算的時(shí)間甚至比使用CPU計(jì)算的時(shí)間更長(zhǎng).所以要避免較小的計(jì)算規(guī)模,充分利用每個(gè)SM的計(jì)算能力.

單純考慮計(jì)算規(guī)模還不夠,block的大小最終影響SM上可運(yùn)行的線程數(shù)量， 所以要選擇合適的block維度,不僅要保證SM上可運(yùn)行的線程數(shù)量,還要兼顧存儲(chǔ)的使用情況.

2.1.2 存儲(chǔ)的使用 對(duì)于global memory,如果可執(zhí)行異步拷貝,在完成操作前,主機(jī)端就會(huì)獲得返回值,從而CPU線程可繼續(xù)執(zhí)行余下操作， 所以可采用異步拷貝達(dá)到加速的目的， 同時(shí)應(yīng)該合理布局global memory上的數(shù)據(jù)影響GPU的表現(xiàn).

在每個(gè)SM中,register是有限的.當(dāng)超過(guò)這個(gè)限制時(shí),SM就會(huì)使用local memory存儲(chǔ).因?yàn)閘ocal memory的訪問速度要慢于register的訪問速度,所以應(yīng)合理聲明變量以減少local memory的使用.

在Fermi架構(gòu)中,有一塊64 Kb的RAM,可分成16 Kb的cache和48 Kb的shared memory,或48 Kb的cache和16 Kb的shared memory.前者有更大的shared memory,相應(yīng)的cache較小,cache命中率就要降低;后者中,cache有足夠的空間提高命中率.可針對(duì)不同的情況選擇不同的劃分方式以提高運(yùn)行速度.

2.1.3 控制指令流的使用 控制流(if，switch，while和for)使線程跳到不同的分支.GPU在這種情況下,不同執(zhí)行路徑必須串行執(zhí)行,且在所有的分支執(zhí)行結(jié)束后,所有的線程才能回到相同的執(zhí)行路徑上.所以應(yīng)減少分支數(shù),可將一些控制流移到主機(jī)端執(zhí)行,對(duì)于無(wú)法在主機(jī)端執(zhí)行的控制流進(jìn)行適當(dāng)合并.

2.2 設(shè)計(jì)實(shí)現(xiàn)

本文給出算法1的主要實(shí)現(xiàn)過(guò)程.計(jì)算能力為2.x的GPU技術(shù)參數(shù)如下：每個(gè)SM上可運(yùn)行的最大線程數(shù)為1 536； 每個(gè)SM上可運(yùn)行的最大block數(shù)為8； 32位的register數(shù)量為32 K.

2.2.1 并行部分的實(shí)現(xiàn) 本文共實(shí)現(xiàn)4個(gè)經(jīng)典的enchmark test函數(shù),包括Sphere,Rastrigin,Griewangk和Rosenbrock.基于本文的加速策略,對(duì)算法1的步驟3),編寫4個(gè)不同的kernel函數(shù)分別對(duì)應(yīng)4個(gè)enchmark test函數(shù),在主機(jī)端執(zhí)行switch語(yǔ)句選擇不同的kernel函數(shù).對(duì)于算法1的步驟4)～6),4個(gè)enchmark test函數(shù)的操作是一致的,所以分別編寫一個(gè)kernel函數(shù).在kernel函數(shù)里,尤其是步驟5)中,不可避免使用if語(yǔ)句,要盡量合并這些if語(yǔ)句.

2.2.2 線程的分布和shared memory的選擇 在算法1的并行部分中,用一個(gè)線程處理一個(gè)粒子,即一個(gè)線程代表一個(gè)粒子.在計(jì)算能力為2.x的GPU中,每個(gè)SM上可運(yùn)行的最大線程數(shù)是1 536,SM的數(shù)量是15.所以為了最大限度地利用SM的計(jì)算能力,一次至少需要23 040個(gè)粒子.

在執(zhí)行并行部分前,使用標(biāo)志cudaFunc-CachePreferL1將RAM分成16 Kb cache 和 48 Kb shared memory,或 48 Kb cache和16 Kb shared memory,分別用存儲(chǔ)1和存儲(chǔ)2表示.在算法1的步驟6)中,需要使用最多的shared memory,設(shè)計(jì)使用5個(gè)大小為D×N的float類型數(shù)組,用于存儲(chǔ)適應(yīng)值和位置數(shù)據(jù).為兼顧SM的計(jì)算能力和存儲(chǔ)能力,并方便實(shí)驗(yàn)對(duì)比,將block的維度設(shè)計(jì)成3種：8×8,16×8 和 16×16,分別用維度1、 維度2和維度3表示.

2.2.3 Global memory和變量的使用 在算法1的步驟2)中,需要依次拷貝每個(gè)粒子的適應(yīng)值、 最優(yōu)適應(yīng)值、 位置、 最優(yōu)位置和速度及粒子群的最優(yōu)適應(yīng)值和最優(yōu)位置,同時(shí)還有計(jì)算時(shí)所需的隨機(jī)數(shù)， 使用cudaMemcpyAsync( )依次完成上述數(shù)據(jù)的異步拷貝.

設(shè)計(jì)數(shù)據(jù)在global memory中的存儲(chǔ)方式,按維度依次存儲(chǔ)每個(gè)粒子,即可認(rèn)為是邏輯上的二維數(shù)組.在邏輯二維數(shù)組中,第j個(gè)粒子的第i維數(shù)據(jù)下標(biāo)為(i,j),表示在大小為D×N的數(shù)組中下標(biāo)為(i×N+j)的數(shù)據(jù).這樣的設(shè)計(jì)保證在同一個(gè)block中可同時(shí)計(jì)算每個(gè)粒子及其維度.

在CUDA中有一些內(nèi)建的變量如threadid和blockid,在kernel函數(shù)中,避免不必要的變量聲明并最大限度地使用內(nèi)建變量.

2.2.4 粒子的限制 使用

(1)

計(jì)算關(guān)于shared memory的使用情況[6].在計(jì)算能力為2.x的GPU中,shared memory最大為48 Kb,結(jié)合式(1),可得:

(2)

不等式(2)的前半部分表示粒子位置、 速度和最優(yōu)位置及粒子群的最優(yōu)位置所占存儲(chǔ)的大小,后半部分表示適應(yīng)值、 粒子最優(yōu)適應(yīng)值和粒子群最優(yōu)適應(yīng)值所占存儲(chǔ)的大小.一個(gè)float類型數(shù)據(jù)的大小是32位,根據(jù)式(2),D的取值是o(103).為實(shí)現(xiàn)“無(wú)限個(gè)”粒子,在設(shè)計(jì)中通過(guò)使用循環(huán)完成.

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

實(shí)驗(yàn)條件：Windows XP操作系統(tǒng),Visual Studio 2005編程工具,cudaToolkit3.1.GTX480,Inter(R)Core(TM)2 Duo CPU E5700@2.93 GHz和Dell Optiplex360.

實(shí)驗(yàn)選取4個(gè)函數(shù)：Sphere,Rastrigin,Griewangk和Rosenbrock， 分別表示如下：

參數(shù)設(shè)定：N=105,D=50,Iter=5 000;時(shí)間1表示文獻(xiàn)[10]的運(yùn)行時(shí)間,時(shí)間2表示多核CPU運(yùn)行時(shí)間,時(shí)間3表示本文優(yōu)化后CPU的運(yùn)行時(shí)間； 數(shù)值1表示CPU上的適應(yīng)值,數(shù)值2表示多核CPU上的適應(yīng)值， 數(shù)值3表示GPU上的適應(yīng)值.

3.1 Block維度和shared memory的影響

在2.2.2中所述的3種block維度選擇和2種RAM劃分情況下,在f1～f4函數(shù)上運(yùn)行SPSO算法,運(yùn)行時(shí)間分別列于表1～表4.

表1 Sphere函數(shù)上的運(yùn)行結(jié)果Table 1 Results of Sphere function

表2 Rastrigin函數(shù)上的運(yùn)行結(jié)果Table 2 Results of Rastrigin function

表3 Griewangk函數(shù)上的運(yùn)行結(jié)果Table 3 Results of Griewangk function

表4 Rosenbrock函數(shù)上的運(yùn)行結(jié)果Table 4 Results of Rosenbrock function

實(shí)驗(yàn)結(jié)果表明,當(dāng)block的維度為16×8,RAM劃分成48 Kb的cache和16 Kb的shared memory時(shí),所用的時(shí)間最短.原因在于48 Kb的cache和16 Kb的shared memory劃分方式,可極大提高cache的命中率,同時(shí)在這種block維度下,單位時(shí)間內(nèi)SM上所運(yùn)行的實(shí)際線程為1 024,所使用的shared memory大小為15 Kb,這樣不僅可充分利用SM的計(jì)算能力,還可充分利用shared memory.

3.2 算法的運(yùn)行時(shí)間

文獻(xiàn)[10]在GPU上實(shí)現(xiàn)了SPSO算法,在相同條件下,本文算法與文獻(xiàn)[10]算法的對(duì)比結(jié)果列于表5～表7.

表5 Global memory優(yōu)化結(jié)果Table 5 Performance of the optimization of global memory

表6 Register files優(yōu)化結(jié)果Table 6 Performance of the optimization of register

表7 控制指令流優(yōu)化結(jié)果Table 7 Performance of the optimization of control instructions stream

由表5～表7可見：

1) global memory上的優(yōu)化實(shí)現(xiàn)并沒有得到很好的結(jié)果.由算法1可知,完成數(shù)據(jù)拷貝后,SPSO算法由GPU完成.在數(shù)據(jù)依次拷貝中和拷貝后在CPU上進(jìn)行的操作很少,異步傳輸?shù)男Ч⒉幻黠@.

2) register上的優(yōu)化實(shí)現(xiàn)沒有得到很好的結(jié)果.因?yàn)橛?jì)算中中間變量的使用是不可控制的.

3) 控制指令流的優(yōu)化實(shí)現(xiàn)得到了很好的結(jié)果.因?yàn)樗惴?中步驟3)沒有采用一個(gè)kernel函數(shù)里使用switch語(yǔ)句選擇4個(gè)enchmark test函數(shù)的方式,而是編寫了4個(gè)kernel函數(shù),且控制流主要由CPU執(zhí)行,GPU運(yùn)行算法的并行部分,這樣充分利用了GPU的并行性.

應(yīng)用本文中的block維度選擇和RAM劃分情況,結(jié)合上述優(yōu)化實(shí)現(xiàn),最終的加速結(jié)果列于表8.在Rastrigin和Griewangk上運(yùn)行SPSO算法與多核進(jìn)行比較,加速結(jié)果列于表9.通過(guò)式(2),本文得到了1 000維的粒子， 時(shí)間結(jié)果列于表10.函數(shù)適應(yīng)值(fitness)隨迭代次數(shù)的變化情況列于表11和表12.

表8 最終優(yōu)化結(jié)果Table 8 Final optimized results

表9 與多核的比較結(jié)果Table 9 Results of optimization compared with those by multi-core

表10 Rastrigin和Griewangk所需時(shí)間比較結(jié)果Table 10 Time comparison results by Rastrigin and Griewangk

表11 Rastrigin的fitness隨迭代的變化情況(N=2 000,D=50)Table 11 Fitness changes of Rastrigin function with iteration (N=2 000,D=50)

綜上所述， 本文實(shí)現(xiàn)了Fermi架構(gòu)上加速SPSO算法.首先,通過(guò)block的3種劃分和shared memory的2種情況,衡量Fermi架構(gòu)中SM的shared memory使用和SM的計(jì)算能力； 其次,在global memory 和register及指令流的使用上進(jìn)行優(yōu)化.實(shí)驗(yàn)結(jié)果表明,當(dāng)block為16×8,shared memory為16 Kb,同時(shí)使用異步傳輸,減少并合并分支,最大限度使用register時(shí),在函數(shù)f2上的加速比為1.579,在函數(shù)f3上的加速比為1.919.實(shí)驗(yàn)中實(shí)現(xiàn)了1 000維的粒子.

表12 Griewangk的fitness隨迭代的變化情況(N=2 000,D=50)Table 12 Fitness changes of Griewangk function with iteration (N=2 000,D=50)

[1] Kennedy J,Eberhart R C.Particle Swarm Optimization [C]//Proc IEEE International Conference on Neural Networks.Piscataway: IEEE Press,1995: 1942-1948.

[2] Shayeghi A,Shayeghi H,Eimani K H.Application of PSO Technique for Seismic Control of Tall Building [J].International Journal of Electrical and Computer Engineering,2009,4(5): 293-300.

[3] Das M T,Dulger L C.Signature Verification (SV) Toolbox: Application of PSO-NN [J].Engineering Applications of Artificial Intelligence,2009,22(4/5): 688-694.

[4] FANG Hong-qing,CHEN Long,SHEN Zu-yi.Application of an Improved PSO Algorithm to Optimal Tuning of PID Gains for Water Turbine Governor [J].Energy Conversion and Management,2011,52(4): 1763-1770.

[5] Bratton D,Kennedy J.Defining a Standard for Particle Swarm Optimization [C]//IEEE Swarm Intelligence Symposium.Piscataway: IEEE Press,2007: 120-127.

[6] Mussi L,Daolio F,Cagnoni S.Evaluation of Parallel Particle Swarm Optimization Algorithms within the CUDA Architecture [J].Information Sciences,2011,181(20): 4642-4657.

[7] Weber R,Gothandaraman A,Hinde R J,et al.Comparing Hardware Accelerators in Scientific Applications: A Case Study [J].IEEE Transactions on Parallel and Distributed Systems,2011,22(1): 58-68.

[8] NVIDIA Corporation.CUDA Programming Guide for CUDA Toolkit 3.2 [EB/OL].[2012-03-08].http://developer.download.nvidia.com/compute/cuda/3\_2\_prod/toolkit/docs/CUDA\_C\_Programming\_Guide.pdf.

[9] Veronese L P,De,Krohling R A.Swarm’s Flight: Accelerating the Particles Using C-CUDA [C]//Proc IEEE Congress on Evolutionary Computation (CEC 2009).Piscataway: IEEE Press,2009: 3264-3270.

[10] ZHOU You,TAN Ying.GPU-Based Parallel Particle Swarm Optimization [C]//Proc IEEE Congress on Evolutionary Computation (CEC 2009).Piscataway: IEEE Press,2009: 1493-1500.

[11] NVIDIA Corporation.Whitepaper: NVIDIA’s Next Generation CUDA Compute Architecture: Fermi [EB/OL].2009-09-01.http://www.nvidia.com/object/fermi\_architecture.html.