多自由度結(jié)構(gòu)靜風(fēng)響應(yīng)的GPU并行計(jì)算

楊智誠(chéng)+饒瑞

摘要： 根據(jù)某大型雙層柱面網(wǎng)殼風(fēng)致靜力響應(yīng)計(jì)算的有限元模型，建立基于GPU的MATLAB快速并行計(jì)算平臺(tái)，實(shí)現(xiàn)CUDA框架下多自由度結(jié)構(gòu)風(fēng)致靜力位移響應(yīng)的快速求解.數(shù)值計(jì)算表明，與傳統(tǒng)的CPU串行計(jì)算相比，通過(guò)GPU實(shí)現(xiàn)的大型矩陣的求逆、乘法、除法等運(yùn)算速度得到大幅提高，位移計(jì)算獲得23倍的最大加速比；結(jié)果誤差對(duì)比分析也表明基于GPU的計(jì)算結(jié)果能夠滿足工程精度要求.

關(guān)鍵詞： 雙層網(wǎng)殼； 風(fēng)致靜力響應(yīng)； 節(jié)點(diǎn)位移； GPU； MATLAB； CUDA

中圖分類(lèi)號(hào)： TU311.1文獻(xiàn)標(biāo)志碼： B

Abstract： According to the finite element model which is built for the calculation of static windinduced response of a largescale doublelayer cylindrical reticulated shell， a fast parallel computing platform is built in MATLAB by GPU， and the fast solution for the windinduced static displacement responses is implemented with the frame of CUDA. The numerical calculation indicates that， compared with the traditional serial computing of CPU， the calculation speed of the operation of inversion， multiplication， division of large matrix is greatly improved， and the maximum speedup ratio of displacement calculation is improved by 23 times. The comparison and analysis on the errors show that the results based on GPU calculation can meet the accuracy requirements of engineering projects.

Key words： doublelayer reticulated shell； static windinduced response； nodal displacement； GPU； MATLAB； CUDA

收稿日期： 2015[KG*9〗03[KG*9〗26修回日期： 2015[KG*9〗04[KG*9〗22

基金項(xiàng)目： 廣州市屬高?！把虺菍W(xué)者”項(xiàng)目（12A004S）； 高等學(xué)校博士點(diǎn)基金（博導(dǎo)類(lèi)）（20124410110005）； 珠江科技新星專(zhuān)項(xiàng)（2012J2200093）

作者簡(jiǎn)介： 楊智誠(chéng)（1992—），男，廣東湛江人，碩士研究生，研究方向?yàn)榻Y(jié)構(gòu)抗風(fēng)設(shè)計(jì)，（Email）328157170@qq.com0引言

隨著現(xiàn)代建筑結(jié)構(gòu)的復(fù)雜程度越來(lái)越高，采用有限元法進(jìn)行結(jié)構(gòu)分析時(shí)，大量自由度導(dǎo)致大型高階剛度矩陣出現(xiàn).在結(jié)構(gòu)位移響應(yīng)的求解計(jì)算過(guò)程中，大型矩陣運(yùn)算成為影響計(jì)算效率的瓶頸.MATLAB是目前工程計(jì)算中自編程序的主流工具，可通過(guò)矩陣求逆或左除求解結(jié)構(gòu)靜力響應(yīng)，但大型矩陣常使該類(lèi)計(jì)算在普通電腦上出現(xiàn)因內(nèi)存溢出而無(wú)法計(jì)算或耗時(shí)較長(zhǎng)等問(wèn)題.

基于GPU的并行計(jì)算是指在NVIDIA公司推出的統(tǒng)一計(jì)算設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）下，通過(guò)調(diào)用GPU多線程和函數(shù)庫(kù)完成單指令的并發(fā)執(zhí)行，以實(shí)現(xiàn)大數(shù)據(jù)快速處理的計(jì)算方法.[1]與傳統(tǒng)的串行計(jì)算相比，使用GPU并行計(jì)算能夠大幅度提高計(jì)算效率且在普通個(gè)人電腦上就能實(shí)現(xiàn).MATLAB的主要計(jì)算模式是基于CPU的串行計(jì)算，為充分利用GPU的計(jì)算優(yōu)勢(shì)，其2010b版本開(kāi)始支持GPU計(jì)算.[2]劉紹波等[3]實(shí)現(xiàn)MATLAB在CUDA框架下對(duì)GPU擴(kuò)展使用，并應(yīng)用于巖土工程計(jì)算.張健飛等[4]將GPU應(yīng)用于有限元計(jì)算中，在CUDA框架下實(shí)現(xiàn)預(yù)處理共軛梯度法求解大規(guī)模線性方程組.

本文基于某大跨雙層柱面網(wǎng)殼結(jié)構(gòu)（桿件總數(shù)為10 080個(gè)）的風(fēng)致靜力位移求解的有限元模型，研究在MATLAB平臺(tái)上通過(guò)GPU加速大型結(jié)構(gòu)靜力響應(yīng)求解的方法，并分析不同算法的求解效率和誤差.

1工程背景和需加速函數(shù)的定位

1.1工程背景

某發(fā)電廠干煤棚為雙層柱面網(wǎng)殼結(jié)構(gòu)，縱向長(zhǎng)度為140 m，橫向跨度為103 m，高為40 m，沿縱向邊緣為剛性支承，見(jiàn)圖1.該結(jié)構(gòu)為正放四角錐柱面網(wǎng)殼，桿件總數(shù)和節(jié)點(diǎn)總數(shù)分別為10 080和2 592個(gè)，桿件彈性模量為206 GPa，密度為7 850 kg/m3.基于風(fēng)洞實(shí)驗(yàn)獲得某一典型風(fēng)向角下的實(shí)際風(fēng)壓時(shí)程[56]，在此基礎(chǔ)上計(jì)算等效靜力風(fēng)載荷[7].

根據(jù)有限元法，結(jié)構(gòu)的節(jié)點(diǎn)位移求解公式[8]為U=K-1·F（1）式中：U為節(jié)點(diǎn)位移矩陣，階數(shù)為7 776×1，即結(jié)構(gòu)總自由度數(shù)為7 776；K-1為總體剛度矩陣的逆矩陣，階數(shù)為7 776×7 776；F為節(jié)點(diǎn)的等效靜力風(fēng)載荷.

MATLAB的相應(yīng)實(shí)現(xiàn)代碼為U=inv（K）*F （2）式中：inv為MATLAB求逆函數(shù).

1.2確定需加速的函數(shù)

在MATLAB平臺(tái)上編制相應(yīng)有限元計(jì)算程序，包含以下主要步驟.

1）讀入相關(guān)數(shù)據(jù)，如節(jié)點(diǎn)數(shù)、節(jié)點(diǎn)坐標(biāo)數(shù)組、節(jié)點(diǎn)力數(shù)組等.

2）計(jì)算節(jié)點(diǎn)力數(shù)組F.

3）計(jì)算單元?jiǎng)偠染仃噆，函數(shù)名為SpaceTrussElementStiffness.

4）組裝總體剛度矩陣K，函數(shù)名為SpaceTrussAssemble.

5）引入邊界條件[9]，并計(jì)算節(jié)點(diǎn)位移矩陣U.

通過(guò)MATLAB的代碼效率調(diào)試工具Profiler獲得主要函數(shù)的運(yùn)行時(shí)間、調(diào)用次數(shù)等信息，統(tǒng)計(jì)耗時(shí)前3位的函數(shù)，包括代碼行序、函數(shù)名稱(chēng)、調(diào)用次數(shù)、總耗時(shí)及占整個(gè)耗時(shí)的比例，見(jiàn)表1.顯然，計(jì)算節(jié)點(diǎn)位移﹝式（2）﹞的運(yùn)行時(shí)間最長(zhǎng)，占總計(jì)算時(shí)間的94.7%，因此考慮將最耗時(shí)部分交由GPU實(shí)現(xiàn)，從而提高求解速度.

次數(shù)耗時(shí)/

s時(shí)間

基于GPU的求解將通過(guò)2種方式實(shí)現(xiàn)：

1）將串行代碼中的矩陣求逆與乘法分別交由GPU實(shí)現(xiàn)，主要代碼為K-1=g_inv（K）

U=g_mul（k_inv，F(xiàn)） （3）式中：g_inv和g_mul均為自編的調(diào)用GPU實(shí)現(xiàn)矩陣運(yùn)算的函數(shù)；g_inv為矩陣求逆運(yùn)算；g_mul為矩陣乘法運(yùn)算.

2）將矩陣求逆和乘法轉(zhuǎn)換為矩陣除法，然后將矩陣除法交由GPU實(shí)現(xiàn)，主要代碼為U=g_div_l（K，F(xiàn)） （4）式中：g_div_l為自定義的調(diào)用GPU實(shí)現(xiàn)矩陣左除的函數(shù).

通過(guò)以上2種方法實(shí)現(xiàn)在GPU下的并行求解.

2GPU算法實(shí)現(xiàn)

2.1軟硬件環(huán)境

計(jì)算在普通個(gè)人電腦上完成，其基本配置如下.

在MATLAB平臺(tái)上，可以自編多線程CUDA程序或調(diào)用NVIDIA公司提供的函數(shù)庫(kù)（如CUBLAS和CULA）實(shí)現(xiàn)GPU加速.自編CUDA程序在VC平臺(tái)下編寫(xiě)，因此需要實(shí)現(xiàn)CUDA程序與MATLAB的對(duì)接.NVIDIA公司提供NVMEX插件，使得在MATLAB平臺(tái)上能夠編譯CUDA程序，以生成可被MATLAB調(diào)用的動(dòng)態(tài)鏈接庫(kù).

2.2算法實(shí)現(xiàn)

2.2.1矩陣乘法

基于GPU的矩陣乘法主要通過(guò)調(diào)用CUBLAS函數(shù)庫(kù)實(shí)現(xiàn)，其中包含大量已經(jīng)編寫(xiě)好的關(guān)于線性代數(shù)運(yùn)算的程序.

首先在CUDA程序中編寫(xiě)MEXFunction接口函數(shù)，將MATLAB中需要處理的數(shù)據(jù)傳入CUDA程序，再調(diào)用CUBLAS函數(shù)庫(kù)中的CUBLASSgemm函數(shù)實(shí)現(xiàn)矩陣相乘.根據(jù)NVIDIA CUBLAS Library文檔說(shuō)明，該函數(shù)主要進(jìn)行的運(yùn)算是C=αAB+βC（5）式中：A和B為待相乘的矩陣；C為結(jié)果矩陣；α和β為常因數(shù)，取α=1，β=0.

算法實(shí)現(xiàn)過(guò)程中先調(diào)用CUBLASCreate函數(shù)載入CUBLAS函數(shù)庫(kù)，調(diào)用CUDAMalloc函數(shù)為A，B和C開(kāi)辟顯存空間，并通過(guò)CUBLASSetMatrix函數(shù)將A和B輸入顯存空間，然后調(diào)用CUBLASSgemm函數(shù)實(shí)現(xiàn)矩陣相乘，計(jì)算完成后通過(guò)CUBLASGetMatrix函數(shù)將C輸出到主機(jī)內(nèi)存，最后通過(guò)CUDAFree函數(shù)釋放顯存空間，并調(diào)用CUBLASDestroy函數(shù)釋放所有與CUBLAS函數(shù)庫(kù)內(nèi)容相關(guān)的資源.

2.2.2矩陣求逆

對(duì)非奇異矩陣A和單位矩陣B，用高斯消元方法能夠?qū)/B轉(zhuǎn)化成B/A-1.[10]從CUDA架構(gòu)的編程思路出發(fā)，將算法中的指令盡量設(shè)計(jì)為少分支、可多線程并行執(zhí)行模式，以便于GPU進(jìn)行成塊數(shù)據(jù)運(yùn)算.矩陣求逆主要包括歸一和消元2種操作[11]，算法描述如下.

1）在MATLAB中對(duì)矩陣A和B初始化并通過(guò)接口函數(shù)將數(shù)據(jù)輸入CUDA程序中.

2）在GPU中為A和B分配顯存空間，并把A和B拷貝到GPU的顯存空間中.

3）調(diào)用內(nèi)核函數(shù)[12]計(jì)算.

4）計(jì)算完成后，將結(jié)果從GPU顯存拷貝到到主機(jī)內(nèi)存中.

5）釋放顯存空間.

其中，內(nèi)核函數(shù)的主要求解步驟如下.

1）用一組臨時(shí)標(biāo)量記錄A中第i行的對(duì)角線元素.

另外，CUDA函數(shù)庫(kù)提供矩陣求逆運(yùn)算函數(shù)，可通過(guò)編寫(xiě)CUDA程序直接調(diào)用該函數(shù).基于以上方法，能夠?qū)崿F(xiàn)大型矩陣在GPU下的求逆運(yùn)算.

2.2.3矩陣除法

在MATLAB對(duì)GPU計(jì)算的原生支持中已經(jīng)支持矩陣除法運(yùn)算，因此只需通過(guò)GPUArray函數(shù)將待處理數(shù)據(jù)傳入GPU內(nèi)存中，然后重新加載MATLAB中的除法運(yùn)算，通過(guò)Gather函數(shù)將計(jì)算結(jié)果傳出到主機(jī)內(nèi)存，由MATLAB執(zhí)行顯示結(jié)果，實(shí)現(xiàn)過(guò)程代碼為G_A=GPUArray（A）

A和B為待除矩陣；C為結(jié)果矩陣；G_A，G_B和G_C為與A，B和C相對(duì)應(yīng)的在顯存空間中的數(shù)據(jù)；GPUArray為MATLAB中將數(shù)據(jù)傳入顯存空間的函數(shù)；Gather為與GPUArray相對(duì)應(yīng)的將數(shù)據(jù)從顯存空間傳出到主機(jī)內(nèi)存的函數(shù).

3計(jì)算效率與誤差分析

3.1GPU單精度和雙精度浮點(diǎn)運(yùn)算能力

根據(jù)浮點(diǎn)運(yùn)算能力的計(jì)算公式，核數(shù)×主頻×每個(gè)時(shí)鐘周期可執(zhí)行的浮點(diǎn)數(shù)運(yùn)算次數(shù)，計(jì)算得到GTX760的單精度浮點(diǎn)運(yùn)算能力為2.21 TFLOPS，雙精度浮點(diǎn)運(yùn)算能力為0.11 TFLOPS.[13]

雙精度浮點(diǎn)運(yùn)算能力只有單精度的1/20，因此為獲得最大的加速效率，基于GPU的數(shù)據(jù)處理采用單精度數(shù)據(jù).然而，MATLAB默認(rèn)的計(jì)算數(shù)據(jù)類(lèi)型為雙精度類(lèi)型.因此，為保證結(jié)果的準(zhǔn)確性，在獲得效率提高的同時(shí)需要分析基于GPU（單精度類(lèi)型數(shù)據(jù)）的結(jié)果與基于CPU（雙精度類(lèi)型數(shù)據(jù)）的結(jié)果之間的誤差.若誤差在允許范圍內(nèi)，則表明基于GPU的計(jì)算結(jié)果能夠用于工程計(jì)算.

3.2效率分析

通過(guò)有限元法計(jì)算雙層網(wǎng)殼的總體剛度矩陣、等效靜力風(fēng)載荷，并根據(jù)式（1）求解節(jié)點(diǎn)位移，計(jì)算過(guò)程分別由式（2）～（4）完成.同時(shí)定義計(jì)算效率測(cè)量指標(biāo)為并行加速比=串行計(jì)算時(shí)間/并行計(jì)算時(shí)間.位移計(jì)算結(jié)果見(jiàn)圖2～4.

由表2可以看出，基于CPU的MATLAB求逆和乘法求解節(jié)點(diǎn)位移的耗時(shí)最長(zhǎng)，而式（3）將求逆和乘法通過(guò)GPU實(shí)現(xiàn)，獲得一定程度上的加速，為原來(lái)串行計(jì)算的6倍左右，最后通過(guò)改進(jìn)算式，將求逆和乘法轉(zhuǎn)換成矩陣除法并且通過(guò)GPU實(shí)現(xiàn)矩陣除法運(yùn)算，獲得23倍的加速比，可大大提高計(jì)算效率.CPU和GPU的硬件結(jié)構(gòu)簡(jiǎn)圖見(jiàn)圖5.由此可以看出：CPU中存在大量的控制單元Control和計(jì)算單元ALU，同時(shí)還存在大塊的高級(jí)緩存Cache，其中計(jì)算單元所占比例不大.這是由于CPU是串行執(zhí)行模式，需要通過(guò)控制單元與快速緩存決定指令流與數(shù)據(jù)的連續(xù)處理，而GPU是單指令多線程模式，對(duì)單指令能夠同時(shí)并發(fā)處理，且其中存在大量的計(jì)算單元，盡管其主頻遠(yuǎn)低于CPU，但是當(dāng)處理的數(shù)據(jù)量較大時(shí)，GPU的并行處理效率遠(yuǎn)遠(yuǎn)超過(guò)CPU的串行處理效率.

GPU在高性能計(jì)算中，大數(shù)據(jù)對(duì)存儲(chǔ)容量的要求較高.本算例中使用顯卡GTX760的顯存容量為2 048 MB，其使用分為顯示和計(jì)算2部分.在正常顯示下，GPU的內(nèi)存使用量約為300 MB，而通過(guò)GPU實(shí)現(xiàn)快速計(jì)算時(shí)，大量數(shù)據(jù)將拷貝到顯存中，在計(jì)算式（3）和（4）時(shí)，顯存使用量分別達(dá)到1 680 MB和1 210 MB，占用量大.因此，在進(jìn)行更大規(guī)模計(jì)算時(shí)，必須考慮單塊顯卡的存儲(chǔ)容量限制，當(dāng)數(shù)據(jù)量過(guò)大導(dǎo)致無(wú)法放置于單塊顯卡時(shí)，需要考慮使用高性能計(jì)算顯卡或多GPU節(jié)點(diǎn)，從而實(shí)現(xiàn)大數(shù)據(jù)在GPU上的處理.

3.3誤差分析

以式（2）的求解結(jié)果為參照，分析式（3）和（4）得到的節(jié)點(diǎn)位移的相對(duì)誤差E，以對(duì)數(shù)坐標(biāo)表示，見(jiàn)圖6和7.

誤差分析發(fā)現(xiàn)，基于GPU的式（3）與式（4）求解節(jié)點(diǎn)位移的誤差相近，式（3）的誤差為4.41×10-4，式（4）的誤差為4.30×10-4，均在10-4左右，因此可以認(rèn)為基于GPU的節(jié)點(diǎn)位移求解結(jié)果與CPU計(jì)算結(jié)果基本一致，誤差在允許范圍內(nèi)，計(jì)算結(jié)果能夠滿足工程精度.

4結(jié)論

以等效靜力風(fēng)載荷下某大跨雙層柱面網(wǎng)殼的靜力位移求解為工程背景，闡述涉及大型矩陣結(jié)構(gòu)計(jì)算的GPU并行實(shí)現(xiàn)過(guò)程，通過(guò)分析得出以下結(jié)論.

1）基于GPU的快速并行計(jì)算能夠大幅提高計(jì)算效率，尤其是對(duì)大數(shù)據(jù)量的單指令運(yùn)算能夠表現(xiàn)出明顯的優(yōu)勢(shì).

2）由GPU計(jì)算所得到的單精度數(shù)據(jù)結(jié)果與MATLAB計(jì)算的雙精度數(shù)據(jù)結(jié)果存在一定誤差，但誤差較小，能夠滿足實(shí)際工程計(jì)算要求.

3）節(jié)點(diǎn)位移求解是各類(lèi)結(jié)構(gòu)計(jì)算中的重要部分，可能成為結(jié)構(gòu)設(shè)計(jì)效率或優(yōu)化分析效率的瓶頸，基于GPU的并行計(jì)算方法能夠快速、準(zhǔn)確地為后續(xù)結(jié)構(gòu)分析提供有效的計(jì)算數(shù)據(jù).參考文獻(xiàn)：

[1]SHANE C. CUDA并行程序設(shè)計(jì)GPU編程指南[M]. 蘇統(tǒng)華， 李東， 李松澤， 等， 譯. 北京： 機(jī)械工業(yè)出版社， 2012： 1069.

[2]劉維. 實(shí)戰(zhàn)MATLAB之并行程序設(shè)計(jì)[M]. 北京： 北京航空航天大學(xué)出版社， 2012： 10120.

[3]劉紹波， 劉明貴， 張國(guó)華. 基于CUDA的加速M(fèi)ATLAB計(jì)算研究[J]. 計(jì)算機(jī)應(yīng)用研究， 2010， 27（6）： 21402143.

LIU Shaobo， LIU Minggui， ZHANG Guohua. Model of accelerating MATLAB computation based on CUDA[J]. Appl Res Comput， 2010， 27（6）： 21402143.

[4]張健飛， 沈德飛. 有限元GPU加速計(jì)算的實(shí)現(xiàn)方法[J]. 計(jì)算機(jī)輔助工程， 2014， 23（2）： 4145.

ZHANG Jianfei， SHEN Defei. Implementation method of GPUaccelerated finite element calculation[J]. Comput Aided Eng， 2014， 23（2）： 4145.

[5]葉孟洋. 大跨度干煤棚結(jié)構(gòu)風(fēng)載風(fēng)洞試驗(yàn)研究[D]. 上海： 同濟(jì)大學(xué)， 2006.

[6]黃友欽. 風(fēng)雪共同作用下大跨度屋蓋結(jié)構(gòu)的動(dòng)力穩(wěn)定[D]. 上海： 同濟(jì)大學(xué)， 2010.

[7]黃友欽， 岳啟哲， 傅繼陽(yáng)， 等. 基于應(yīng)力比法的雙層網(wǎng)殼結(jié)構(gòu)抗風(fēng)優(yōu)化[J]. 結(jié)構(gòu)工程師， 2013， 29（4）： 112117.HUANG Youqin， YUE Qizhe， FU Jiyang， et al. Wind resistant optimization of doublelayer reticulated shells by the stress ratio method[J]. Struct Engineers， 2013， 29（4）： 112117.

[8]楊茀康， 李家寶. 結(jié)構(gòu)力學(xué)[M]. 4版. 北京： 高等教育出版社， 2006： 18.

[9]王煥定， 焦兆平. 有限單元法基礎(chǔ)[M]. 2版. 北京： 高等教育出版社， 2010： 19.

[10]李慶揚(yáng). 數(shù)值分析[M]. 5版. 北京： 清華大學(xué)出版社， 2008： 11.

[11]劉麗， 沈杰， 李洪林. 基于GPU的矩陣求逆性能測(cè)試和分析[J]. 華東理工大學(xué)學(xué)報(bào)（自然科學(xué)版）， 2010， 36（6）： 812817.

LIU Li， SHEN Jie， LI Honglin. Performance testing and analysis for matrix inversion base on GPU[J]. J East China Univ Sci & Technol（Nat Sci）， 2010， 36（6）： 812817.

[12]RYOO S， RODRIGUES C I， BAGHSORKHI S S， et al. Optimization principles and application performance evaluation of a multithreaded GPU using CUDA[C] //Proc ACM SIGPLAN Symrposium Principles & Practice Parallel Programming. New York： ACM Press， 2008： 7382.

[13]AHMED E Z， ERIC M C， ALISTAIR R， et al. Performance evaluation of the NVIDIA GeForce8800 GTX GPU for machine learning[C] //Proc 8th Int Conf Comput Sci. New York： Springer Verlag， 2009： 466475.