田 冪,胡 亮,車喜龍

(吉林大學 計算機科學與技術學院,長春 130012)

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

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

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

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

1.1 SPSO算法

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

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

1.2 Fermi架構(gòu)

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

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

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

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

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

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

算法1

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

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

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

3) 計算每個粒子的適應值；

4) 更新每個粒子的最佳適應值和最佳位置；

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

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

end for

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

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

2.1 加速策略

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

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

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

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

在每個SM中,register是有限的.當超過這個限制時,SM就會使用local memory存儲.因為local memory的訪問速度要慢于register的訪問速度,所以應合理聲明變量以減少local memory的使用.

在Fermi架構(gòu)中,有一塊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在這種情況下,不同執(zhí)行路徑必須串行執(zhí)行,且在所有的分支執(zhí)行結(jié)束后,所有的線程才能回到相同的執(zhí)行路徑上.所以應減少分支數(shù),可將一些控制流移到主機端執(zhí)行,對于無法在主機端執(zhí)行的控制流進行適當合并.

2.2 設計實現(xiàn)

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

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

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

在執(zhí)行并行部分前,使用標志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類型數(shù)組,用于存儲適應值和位置數(shù)據(jù).為兼顧SM的計算能力和存儲能力,并方便實驗對比,將block的維度設計成3種：8×8,16×8 和 16×16,分別用維度1、 維度2和維度3表示.

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

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

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

2.2.4 粒子的限制 使用

(1)

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

(2)

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

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

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

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

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

3.1 Block維度和shared memory的影響

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

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

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

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

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

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

3.2 算法的運行時間

文獻[10]在GPU上實現(xiàn)了SPSO算法,在相同條件下,本文算法與文獻[10]算法的對比結(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)化實現(xiàn)并沒有得到很好的結(jié)果.由算法1可知,完成數(shù)據(jù)拷貝后,SPSO算法由GPU完成.在數(shù)據(jù)依次拷貝中和拷貝后在CPU上進行的操作很少,異步傳輸?shù)男Ч⒉幻黠@.

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

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

應用本文中的block維度選擇和RAM劃分情況,結(jié)合上述優(yōu)化實現(xiàn),最終的加速結(jié)果列于表8.在Rastrigin和Griewangk上運行SPSO算法與多核進行比較,加速結(jié)果列于表9.通過式(2),本文得到了1 000維的粒子， 時間結(jié)果列于表10.函數(shù)適應值(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所需時間比較結(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)

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