多自由度結構靜風響應的GPU并行計算

楊智誠+饒瑞

摘要： 根據某大型雙層柱面網殼風致靜力響應計算的有限元模型，建立基于GPU的MATLAB快速并行計算平臺，實現CUDA框架下多自由度結構風致靜力位移響應的快速求解.數值計算表明，與傳統的CPU串行計算相比，通過GPU實現的大型矩陣的求逆、乘法、除法等運算速度得到大幅提高，位移計算獲得23倍的最大加速比；結果誤差對比分析也表明基于GPU的計算結果能夠滿足工程精度要求.

關鍵詞： 雙層網殼； 風致靜力響應； 節點位移； GPU； MATLAB； CUDA

中圖分類號： TU311.1文獻標志碼： 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

基金項目： 廣州市屬高校“羊城學者”項目（12A004S）； 高等學校博士點基金（博導類）（20124410110005）； 珠江科技新星專項（2012J2200093）

作者簡介： 楊智誠（1992—），男，廣東湛江人，碩士研究生，研究方向為結構抗風設計，（Email）328157170@qq.com0引言

隨著現代建筑結構的復雜程度越來越高，采用有限元法進行結構分析時，大量自由度導致大型高階剛度矩陣出現.在結構位移響應的求解計算過程中，大型矩陣運算成為影響計算效率的瓶頸.MATLAB是目前工程計算中自編程序的主流工具，可通過矩陣求逆或左除求解結構靜力響應，但大型矩陣常使該類計算在普通電腦上出現因內存溢出而無法計算或耗時較長等問題.

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

本文基于某大跨雙層柱面網殼結構（桿件總數為10 080個）的風致靜力位移求解的有限元模型，研究在MATLAB平臺上通過GPU加速大型結構靜力響應求解的方法，并分析不同算法的求解效率和誤差.

1工程背景和需加速函數的定位

1.1工程背景

某發電廠干煤棚為雙層柱面網殼結構，縱向長度為140 m，橫向跨度為103 m，高為40 m，沿縱向邊緣為剛性支承，見圖1.該結構為正放四角錐柱面網殼，桿件總數和節點總數分別為10 080和2 592個，桿件彈性模量為206 GPa，密度為7 850 kg/m3.基于風洞實驗獲得某一典型風向角下的實際風壓時程[56]，在此基礎上計算等效靜力風載荷[7].

根據有限元法，結構的節點位移求解公式[8]為U=K-1·F（1）式中：U為節點位移矩陣，階數為7 776×1，即結構總自由度數為7 776；K-1為總體剛度矩陣的逆矩陣，階數為7 776×7 776；F為節點的等效靜力風載荷.

MATLAB的相應實現代碼為U=inv（K）*F （2）式中：inv為MATLAB求逆函數.

1.2確定需加速的函數

在MATLAB平臺上編制相應有限元計算程序，包含以下主要步驟.

1）讀入相關數據，如節點數、節點坐標數組、節點力數組等.

2）計算節點力數組F.

3）計算單元剛度矩陣k，函數名為SpaceTrussElementStiffness.

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

5）引入邊界條件[9]，并計算節點位移矩陣U.

通過MATLAB的代碼效率調試工具Profiler獲得主要函數的運行時間、調用次數等信息，統計耗時前3位的函數，包括代碼行序、函數名稱、調用次數、總耗時及占整個耗時的比例，見表1.顯然，計算節點位移﹝式（2）﹞的運行時間最長，占總計算時間的94.7%，因此考慮將最耗時部分交由GPU實現，從而提高求解速度.

次數耗時/

s時間

基于GPU的求解將通過2種方式實現：

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

U=g_mul（k_inv，F） （3）式中：g_inv和g_mul均為自編的調用GPU實現矩陣運算的函數；g_inv為矩陣求逆運算；g_mul為矩陣乘法運算.

2）將矩陣求逆和乘法轉換為矩陣除法，然后將矩陣除法交由GPU實現，主要代碼為U=g_div_l（K，F） （4）式中：g_div_l為自定義的調用GPU實現矩陣左除的函數.

通過以上2種方法實現在GPU下的并行求解.

2GPU算法實現

2.1軟硬件環境

計算在普通個人電腦上完成，其基本配置如下.

在MATLAB平臺上，可以自編多線程CUDA程序或調用NVIDIA公司提供的函數庫（如CUBLAS和CULA）實現GPU加速.自編CUDA程序在VC平臺下編寫，因此需要實現CUDA程序與MATLAB的對接.NVIDIA公司提供NVMEX插件，使得在MATLAB平臺上能夠編譯CUDA程序，以生成可被MATLAB調用的動態鏈接庫.

2.2算法實現

2.2.1矩陣乘法

基于GPU的矩陣乘法主要通過調用CUBLAS函數庫實現，其中包含大量已經編寫好的關于線性代數運算的程序.

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

算法實現過程中先調用CUBLASCreate函數載入CUBLAS函數庫，調用CUDAMalloc函數為A，B和C開辟顯存空間，并通過CUBLASSetMatrix函數將A和B輸入顯存空間，然后調用CUBLASSgemm函數實現矩陣相乘，計算完成后通過CUBLASGetMatrix函數將C輸出到主機內存，最后通過CUDAFree函數釋放顯存空間，并調用CUBLASDestroy函數釋放所有與CUBLAS函數庫內容相關的資源.

2.2.2矩陣求逆

對非奇異矩陣A和單位矩陣B，用高斯消元方法能夠將A/B轉化成B/A-1.[10]從CUDA架構的編程思路出發，將算法中的指令盡量設計為少分支、可多線程并行執行模式，以便于GPU進行成塊數據運算.矩陣求逆主要包括歸一和消元2種操作[11]，算法描述如下.

1）在MATLAB中對矩陣A和B初始化并通過接口函數將數據輸入CUDA程序中.

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

3）調用內核函數[12]計算.

4）計算完成后，將結果從GPU顯存拷貝到到主機內存中.

5）釋放顯存空間.

其中，內核函數的主要求解步驟如下.

1）用一組臨時標量記錄A中第i行的對角線元素.

另外，CUDA函數庫提供矩陣求逆運算函數，可通過編寫CUDA程序直接調用該函數.基于以上方法，能夠實現大型矩陣在GPU下的求逆運算.

2.2.3矩陣除法

在MATLAB對GPU計算的原生支持中已經支持矩陣除法運算，因此只需通過GPUArray函數將待處理數據傳入GPU內存中，然后重新加載MATLAB中的除法運算，通過Gather函數將計算結果傳出到主機內存，由MATLAB執行顯示結果，實現過程代碼為G_A=GPUArray（A）

A和B為待除矩陣；C為結果矩陣；G_A，G_B和G_C為與A，B和C相對應的在顯存空間中的數據；GPUArray為MATLAB中將數據傳入顯存空間的函數；Gather為與GPUArray相對應的將數據從顯存空間傳出到主機內存的函數.

3計算效率與誤差分析

3.1GPU單精度和雙精度浮點運算能力

根據浮點運算能力的計算公式，核數×主頻×每個時鐘周期可執行的浮點數運算次數，計算得到GTX760的單精度浮點運算能力為2.21 TFLOPS，雙精度浮點運算能力為0.11 TFLOPS.[13]

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

3.2效率分析

通過有限元法計算雙層網殼的總體剛度矩陣、等效靜力風載荷，并根據式（1）求解節點位移，計算過程分別由式（2）～（4）完成.同時定義計算效率測量指標為并行加速比=串行計算時間/并行計算時間.位移計算結果見圖2～4.

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

GPU在高性能計算中，大數據對存儲容量的要求較高.本算例中使用顯卡GTX760的顯存容量為2 048 MB，其使用分為顯示和計算2部分.在正常顯示下，GPU的內存使用量約為300 MB，而通過GPU實現快速計算時，大量數據將拷貝到顯存中，在計算式（3）和（4）時，顯存使用量分別達到1 680 MB和1 210 MB，占用量大.因此，在進行更大規模計算時，必須考慮單塊顯卡的存儲容量限制，當數據量過大導致無法放置于單塊顯卡時，需要考慮使用高性能計算顯卡或多GPU節點，從而實現大數據在GPU上的處理.

3.3誤差分析

以式（2）的求解結果為參照，分析式（3）和（4）得到的節點位移的相對誤差E，以對數坐標表示，見圖6和7.

誤差分析發現，基于GPU的式（3）與式（4）求解節點位移的誤差相近，式（3）的誤差為4.41×10-4，式（4）的誤差為4.30×10-4，均在10-4左右，因此可以認為基于GPU的節點位移求解結果與CPU計算結果基本一致，誤差在允許范圍內，計算結果能夠滿足工程精度.

4結論

以等效靜力風載荷下某大跨雙層柱面網殼的靜力位移求解為工程背景，闡述涉及大型矩陣結構計算的GPU并行實現過程，通過分析得出以下結論.

1）基于GPU的快速并行計算能夠大幅提高計算效率，尤其是對大數據量的單指令運算能夠表現出明顯的優勢.

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

3）節點位移求解是各類結構計算中的重要部分，可能成為結構設計效率或優化分析效率的瓶頸，基于GPU的并行計算方法能夠快速、準確地為后續結構分析提供有效的計算數據.參考文獻：

[1]SHANE C. CUDA并行程序設計GPU編程指南[M]. 蘇統華， 李東， 李松澤， 等， 譯. 北京： 機械工業出版社， 2012： 1069.

[2]劉維. 實戰MATLAB之并行程序設計[M]. 北京： 北京航空航天大學出版社， 2012： 10120.

[3]劉紹波， 劉明貴， 張國華. 基于CUDA的加速MATLAB計算研究[J]. 計算機應用研究， 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]. 計算機輔助工程， 2014， 23（2）： 4145.

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

[5]葉孟洋. 大跨度干煤棚結構風載風洞試驗研究[D]. 上海： 同濟大學， 2006.

[6]黃友欽. 風雪共同作用下大跨度屋蓋結構的動力穩定[D]. 上海： 同濟大學， 2010.

[7]黃友欽， 岳啟哲， 傅繼陽， 等. 基于應力比法的雙層網殼結構抗風優化[J]. 結構工程師， 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]楊茀康， 李家寶. 結構力學[M]. 4版. 北京： 高等教育出版社， 2006： 18.

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

[10]李慶揚. 數值分析[M]. 5版. 北京： 清華大學出版社， 2008： 11.

[11]劉麗， 沈杰， 李洪林. 基于GPU的矩陣求逆性能測試和分析[J]. 華東理工大學學報（自然科學版）， 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.