田 冪,胡 亮,車喜龍
(吉林大學 計算機科學與技術學院,長春 130012)
粒子群優化算法(PSO)具有穩健、 有效、 簡易和適應性廣的特點,應用廣泛[1-4].標準粒子群優化算法(SPSO)[5]是對原始算法的直接擴展.盡管SPSO算法擁有許多優點,但在大規模問題,如在大維度和大群數上,仍需較長時間尋找解決方案.SPSO算法的結構存在內在的并行性,通過直接刪除更新之間的依賴,編寫同步PSO,同時允許粒子不嚴格依賴最新信息[6],放寬同步限制,可在并行計算架構下實現SPSO算法.
圖形處理單元(GPU)因其并行計算機制和快速的浮點操作能力,在科學計算上展現出較強的優勢,因此應用廣泛.與其他加速平臺相比,如CellBE和FPGAs,GPU已在HPC研究領域占主導地位[7],同時使用CUDA開發工具[8]可輕松地進行編程.
利用可編程的GPU加速SPSO算法的可行性與有效性已被證實[9-10],然而這些研究都沒有利用最新的硬件版本[11].Fermi 同樣是雙架構設計,即CUDA計算架構和圖形架構,兩種架構可靈活切換,但Fermi幾乎是重新設計并更注重通用計算架構.本文利用Fermi GPU的同時,在最新架構下實現了并行版本的SPSO算法.

圖1 SPSO算法流程Fig.1 Flowchart of SPSO
SPSO系統初始化是一群隨機的粒子,并通過更新粒子尋找多維解決空間的最優值.每個粒子依據其自己發現的本地最近解決方案進行移動.相鄰最佳解決方案由與其相鄰的左鄰居和右鄰居發現.每次迭代中粒子更新其速度和位置直到滿足終止條件.SPSO算法的基本流程如圖1所示.
程序可分為兩部分:主機端部分和設備端部分.主機端部分在CPU上執行,管理數據的傳輸,上下文的預處理和后處理,調度設備端部分的執行; 設備端部分在GPU上執行,由一個或多個稱為kernel的函數組成, 每個kernel是程序的一個并行單元,由許多并行的線程執行.Fermi架構支持不同kernel的同時執行,它的流水技術可減少應用上下文切換帶來的開銷,同時支持kernel間的快速通信.
線程是設備端部分程序的最小執行單位,以grid和block的形式組成.一個block由許多線程組成,同一block中的線程共享數據并同步它們的執行.同一kernel的一些block組成一個grid.線程以warp的形式進行調度和執行,同一warp內的所有線程可認為是同時運行的,存儲訪問也以warp形式工作.
Fermi GPU以圖形處理集群可擴展陣列(GPC)為基礎,每個GPC由4個流水多重處理器(SM)組成.每個SM有32個CUDA核,16個Load/Store(LD/ST)單元和4個特殊函數單元(SFU).CUDA核由全流水整數算術邏輯單元(ALU)和浮點單元(FPU)組成.一個GPU執行一個或多個grid; 一個SM執行一個或多個block; 一個ALU/FPU/SFU執行一個線程.
在GPU芯片外,有64位的GDDR5 DRAM,可視為是global memory,local memory,constant memory 和texture memory的組合,由主機端緩存和讀寫.Fermi對local memory,shared memory和global memory使用統一的地址空間完成加載和存儲操作.在GPU芯片內,有3種類型的快速存儲器:register files,shared memory和cache,它們由設備端進行讀寫操作.
SPSO算法在GPU上的執行過程如下.
算法1
1) 在global memory中申請所需的變量空間;
2) 將數據由主機內存拷貝到global memory內相應的變量中;
fori=1 to Iter do //并行執行for語句中步驟
3) 計算每個粒子的適應值;
4) 更新每個粒子的最佳適應值和最佳位置;
5) 更新粒子群的最佳適應值和最佳位置;
6) 更新每個粒子的速度和位置;
end for
7) 將結果數據由global memory拷貝回主機內存.
本文使用的主要參數如下:Iter表示SPSO算法的迭代次數;D表示每個粒子的維度;N表示粒子的數目.
在Fermi上加速SPSO算法,需要考慮計算量與block維度的選擇、 存儲的使用和控制指令流的使用.
2.1.1 計算量與block維度的選擇 如果計算規模很小,使用GPU進行計算將不能彌補傳輸數據帶來的延時,即使用GPU計算的時間甚至比使用CPU計算的時間更長.所以要避免較小的計算規模,充分利用每個SM的計算能力.
單純考慮計算規模還不夠,block的大小最終影響SM上可運行的線程數量, 所以要選擇合適的block維度,不僅要保證SM上可運行的線程數量,還要兼顧存儲的使用情況.
2.1.2 存儲的使用 對于global memory,如果可執行異步拷貝,在完成操作前,主機端就會獲得返回值,從而CPU線程可繼續執行余下操作, 所以可采用異步拷貝達到加速的目的, 同時應該合理布局global memory上的數據影響GPU的表現.
在每個SM中,register是有限的.當超過這個限制時,SM就會使用local memory存儲.因為local memory的訪問速度要慢于register的訪問速度,所以應合理聲明變量以減少local memory的使用.
在Fermi架構中,有一塊64 Kb的RAM,可分成16 Kb的cache和48 Kb的shared memory,或48 Kb的cache和16 Kb的shared memory.前者有更大的shared memory,相應的cache較小,cache命中率就要降低;后者中,cache有足夠的空間提高命中率.可針對不同的情況選擇不同的劃分方式以提高運行速度.
2.1.3 控制指令流的使用 控制流(if,switch,while和for)使線程跳到不同的分支.GPU在這種情況下,不同執行路徑必須串行執行,且在所有的分支執行結束后,所有的線程才能回到相同的執行路徑上.所以應減少分支數,可將一些控制流移到主機端執行,對于無法在主機端執行的控制流進行適當合并.
本文給出算法1的主要實現過程.計算能力為2.x的GPU技術參數如下:每個SM上可運行的最大線程數為1 536; 每個SM上可運行的最大block數為8; 32位的register數量為32 K.
2.2.1 并行部分的實現 本文共實現4個經典的enchmark test函數,包括Sphere,Rastrigin,Griewangk和Rosenbrock.基于本文的加速策略,對算法1的步驟3),編寫4個不同的kernel函數分別對應4個enchmark test函數,在主機端執行switch語句選擇不同的kernel函數.對于算法1的步驟4)~6),4個enchmark test函數的操作是一致的,所以分別編寫一個kernel函數.在kernel函數里,尤其是步驟5)中,不可避免使用if語句,要盡量合并這些if語句.
2.2.2 線程的分布和shared memory的選擇 在算法1的并行部分中,用一個線程處理一個粒子,即一個線程代表一個粒子.在計算能力為2.x的GPU中,每個SM上可運行的最大線程數是1 536,SM的數量是15.所以為了最大限度地利用SM的計算能力,一次至少需要23 040個粒子.
在執行并行部分前,使用標志cudaFunc-CachePreferL1將RAM分成16 Kb cache 和 48 Kb shared memory,或 48 Kb cache和16 Kb shared memory,分別用存儲1和存儲2表示.在算法1的步驟6)中,需要使用最多的shared memory,設計使用5個大小為D×N的float類型數組,用于存儲適應值和位置數據.為兼顧SM的計算能力和存儲能力,并方便實驗對比,將block的維度設計成3種:8×8,16×8 和 16×16,分別用維度1、 維度2和維度3表示.
2.2.3 Global memory和變量的使用 在算法1的步驟2)中,需要依次拷貝每個粒子的適應值、 最優適應值、 位置、 最優位置和速度及粒子群的最優適應值和最優位置,同時還有計算時所需的隨機數, 使用cudaMemcpyAsync( )依次完成上述數據的異步拷貝.
設計數據在global memory中的存儲方式,按維度依次存儲每個粒子,即可認為是邏輯上的二維數組.在邏輯二維數組中,第j個粒子的第i維數據下標為(i,j),表示在大小為D×N的數組中下標為(i×N+j)的數據.這樣的設計保證在同一個block中可同時計算每個粒子及其維度.
在CUDA中有一些內建的變量如threadid和blockid,在kernel函數中,避免不必要的變量聲明并最大限度地使用內建變量.
2.2.4 粒子的限制 使用
(1)
計算關于shared memory的使用情況[6].在計算能力為2.x的GPU中,shared memory最大為48 Kb,結合式(1),可得:
(2)
不等式(2)的前半部分表示粒子位置、 速度和最優位置及粒子群的最優位置所占存儲的大小,后半部分表示適應值、 粒子最優適應值和粒子群最優適應值所占存儲的大小.一個float類型數據的大小是32位,根據式(2),D的取值是o(103).為實現“無限個”粒子,在設計中通過使用循環完成.
實驗條件:Windows XP操作系統,Visual Studio 2005編程工具,cudaToolkit3.1.GTX480,Inter(R)Core(TM)2 Duo CPU E5700@2.93 GHz和Dell Optiplex360.
實驗選取4個函數:Sphere,Rastrigin,Griewangk和Rosenbrock, 分別表示如下:
參數設定:N=105,D=50,Iter=5 000;時間1表示文獻[10]的運行時間,時間2表示多核CPU運行時間,時間3表示本文優化后CPU的運行時間; 數值1表示CPU上的適應值,數值2表示多核CPU上的適應值, 數值3表示GPU上的適應值.
在2.2.2中所述的3種block維度選擇和2種RAM劃分情況下,在f1~f4函數上運行SPSO算法,運行時間分別列于表1~表4.

表1 Sphere函數上的運行結果Table 1 Results of Sphere function

表2 Rastrigin函數上的運行結果Table 2 Results of Rastrigin function

表3 Griewangk函數上的運行結果Table 3 Results of Griewangk function

表4 Rosenbrock函數上的運行結果Table 4 Results of Rosenbrock function
實驗結果表明,當block的維度為16×8,RAM劃分成48 Kb的cache和16 Kb的shared memory時,所用的時間最短.原因在于48 Kb的cache和16 Kb的shared memory劃分方式,可極大提高cache的命中率,同時在這種block維度下,單位時間內SM上所運行的實際線程為1 024,所使用的shared memory大小為15 Kb,這樣不僅可充分利用SM的計算能力,還可充分利用shared memory.
文獻[10]在GPU上實現了SPSO算法,在相同條件下,本文算法與文獻[10]算法的對比結果列于表5~表7.

表5 Global memory優化結果Table 5 Performance of the optimization of global memory

表6 Register files優化結果Table 6 Performance of the optimization of register

表7 控制指令流優化結果Table 7 Performance of the optimization of control instructions stream
由表5~表7可見:
1) global memory上的優化實現并沒有得到很好的結果.由算法1可知,完成數據拷貝后,SPSO算法由GPU完成.在數據依次拷貝中和拷貝后在CPU上進行的操作很少,異步傳輸的效果并不明顯.
2) register上的優化實現沒有得到很好的結果.因為計算中中間變量的使用是不可控制的.
3) 控制指令流的優化實現得到了很好的結果.因為算法1中步驟3)沒有采用一個kernel函數里使用switch語句選擇4個enchmark test函數的方式,而是編寫了4個kernel函數,且控制流主要由CPU執行,GPU運行算法的并行部分,這樣充分利用了GPU的并行性.
應用本文中的block維度選擇和RAM劃分情況,結合上述優化實現,最終的加速結果列于表8.在Rastrigin和Griewangk上運行SPSO算法與多核進行比較,加速結果列于表9.通過式(2),本文得到了1 000維的粒子, 時間結果列于表10.函數適應值(fitness)隨迭代次數的變化情況列于表11和表12.

表8 最終優化結果Table 8 Final optimized results

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

表10 Rastrigin和Griewangk所需時間比較結果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)
綜上所述, 本文實現了Fermi架構上加速SPSO算法.首先,通過block的3種劃分和shared memory的2種情況,衡量Fermi架構中SM的shared memory使用和SM的計算能力; 其次,在global memory 和register及指令流的使用上進行優化.實驗結果表明,當block為16×8,shared memory為16 Kb,同時使用異步傳輸,減少并合并分支,最大限度使用register時,在函數f2上的加速比為1.579,在函數f3上的加速比為1.919.實驗中實現了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.