有限元GPU加速計(jì)算的實(shí)現(xiàn)方法

張健飛　沈德飛

摘要：研究基于GPU的有限元求解中的總剛矩陣生成和線性方程組求解問(wèn)題.通過(guò)對(duì)單元著色和分組完成總剛矩陣的生成，并以行壓縮存儲(chǔ)（Compressed Sparse Row，CSR）格式存儲(chǔ)，用預(yù)處理共軛梯度法求解所生成的大規(guī)模線性稀疏方程組.在CUDA（Compute Unified Device Architecture）平臺(tái)上完成程序設(shè)計(jì)，并用GT430 GPU對(duì)彈性力學(xué)的平面問(wèn)題和空間問(wèn)題進(jìn)行試驗(yàn).結(jié)果表明，總剛矩陣生成和方程組求解分別得到最高11.7和8的計(jì)算加速比.

關(guān)鍵詞：GPU計(jì)算； 有限元法； 剛度矩陣； 預(yù)處理共軛梯度法

中圖分類號(hào)： TB115.7；TP311

文獻(xiàn)標(biāo)志碼：B

0 引 言

作為一種求解微分方程或積分方程的微分方法，有限元法[1]以其高度的適應(yīng)性，成為現(xiàn)代工程設(shè)計(jì)和結(jié)構(gòu)分析的重要方法之一，并在土木、水利、汽車、機(jī)械、航空航天、核工業(yè)和大地勘測(cè)等眾多領(lǐng)域應(yīng)用廣泛.隨著科學(xué)技術(shù)的不斷發(fā)展，工程問(wèn)題的規(guī)模和復(fù)雜程度相應(yīng)提高，也對(duì)有限元計(jì)算提出更大規(guī)模、更快速度的要求.有限元法的基本思想是“化整為零，積零為整”，與并行計(jì)算技術(shù)的“分而治之”的基本原則相協(xié)調(diào)，因此，對(duì)于大規(guī)模的有限元結(jié)構(gòu)分析，可將效率低下的串行有限元分析改進(jìn)為并行有限元分析.

自從2006年NVIDIA正式發(fā)布用于通用計(jì)算的統(tǒng)一計(jì)算架構(gòu)CUDA（Compute Unified Device Architecture）平臺(tái)[2]后，圖形處理器的體系架構(gòu)得到迅速發(fā)展和完善，性能得到大幅提高，使得GPU不僅能高效應(yīng)用于計(jì)算機(jī)圖形處理，而且其強(qiáng)大的計(jì)算能力也能很好地適用于高性能計(jì)算[3]，大大推動(dòng)基于GPU通用計(jì)算的研究，并廣泛應(yīng)用于醫(yī)學(xué)成像、生物信息學(xué)、計(jì)算結(jié)構(gòu)力學(xué)、計(jì)算流體力學(xué)、計(jì)算金融、地震勘探、地理信息系統(tǒng)以及電影和動(dòng)畫制作等領(lǐng)域.商業(yè)有限元分析軟件ANSYS和Abaqus中也應(yīng)用GPU技術(shù).[4]

目前，基于GPU的有限元分析在求解由有限元生成的稀疏線性方程組的研究中比較多，主要是由于在求解方程組時(shí)數(shù)據(jù)量大而且計(jì)算比較集中，便于并行運(yùn)算.NATHAN等[5]討論稀疏矩陣的數(shù)據(jù)結(jié)構(gòu)，并探討幾種有效的基于CUDA的高效的稀疏矩陣與向量相乘的方法；AIL等[6]探討針對(duì)多GPU的基于CUDA的快速共軛梯度法，并探討共軛梯度法中最耗時(shí)的稀疏矩陣與向量相乘的操作；JOLDES等[7]研究基于GPU的以六面體單元為主的混合網(wǎng)格及在有限元中尋求穩(wěn)定解的問(wèn)題，并通過(guò)CUDA實(shí)現(xiàn)基于非線性力學(xué)模型的自動(dòng)模擬神經(jīng)外科的過(guò)程；PAWEL等[8]探討基于GPU的三維有限元數(shù)值積分算法和計(jì)算方面的內(nèi)容；CECKA等[9]探討基于GPU的有限元法剛度矩陣組裝方法，評(píng)估每種方法的優(yōu)缺點(diǎn)；李熙銘[10]驗(yàn)證基于CUDA的復(fù)電阻率問(wèn)題，并詳細(xì)研究復(fù)共軛梯度法；胡耀國(guó)[11]運(yùn)用單元分組的方式計(jì)算得出有限元法中的總剛矩陣，并研究基于GPU的共軛梯度法，因當(dāng)時(shí)的資源限制，其在總剛矩陣生成部分的加速效果并不明顯.

本文利用GPU強(qiáng)大的并行計(jì)算能力和CPU的高效邏輯處理能力，將有限元法中計(jì)算量很大的單剛矩陣計(jì)算和總剛矩陣生成及線性方程組的求解交給GPU運(yùn)算，CPU負(fù)責(zé)相應(yīng)數(shù)據(jù)的前處理和后處理及整個(gè)分析過(guò)程中的邏輯關(guān)系.應(yīng)用CUDA平臺(tái)完成基于GPU的有限元分析程序，在本文的計(jì)算平臺(tái)上運(yùn)行，并用彈性力學(xué)平面問(wèn)題和空間問(wèn)題的有限元求解測(cè)試其運(yùn)行效果.在總剛矩陣生成部分，為避免計(jì)算時(shí)線程沖突，采用文獻(xiàn)[9]和[11-12]中都提到的對(duì)結(jié)構(gòu)單元進(jìn)行著色的方法.在求解大規(guī)模線性方程組時(shí)充分利用現(xiàn)有資源，在完全調(diào)用現(xiàn)有庫(kù)函數(shù)的情況下實(shí)現(xiàn)基于GPU的預(yù)條件共軛梯度法.

1 CUDA編程模型

NVIDIA的CUDA并行計(jì)算模型使用C語(yǔ)言作為開(kāi)發(fā)語(yǔ)言，同時(shí)也支持其他編程語(yǔ)言或應(yīng)用程序接口.一個(gè)完整的CUDA程序包含連續(xù)運(yùn)行在CPU端的主程序及用CUDA編寫的運(yùn)行在并行GPU設(shè)備上的內(nèi)核程序.由圖1可知，在CUDA架構(gòu)下，應(yīng)用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）組成.運(yùn)行流程：在CPU端準(zhǔn)備數(shù)據(jù)，然后傳到GPU端進(jìn)行并行計(jì)算，最后將計(jì)算好的數(shù)據(jù)再傳回CPU端.

CPU端的程序主要負(fù)責(zé)數(shù)據(jù)的準(zhǔn)備和一些邏輯運(yùn)算，其中Kernel函數(shù)是整個(gè)應(yīng)用程序的關(guān)鍵.Kernel函數(shù)以線程網(wǎng)格（Grid）的形式組織，每個(gè)線程網(wǎng)格由若干個(gè)線程塊（Block）組成，而每個(gè)線程塊又由若干個(gè)線程（Thread）組成.在執(zhí)行時(shí)，Kernel函數(shù)以線程塊為單位，各線程塊間并行執(zhí)行，不同線程塊間只能通過(guò)全局顯存進(jìn)行數(shù)據(jù)共享，同一線程塊內(nèi)的線程之間可以通過(guò)共享內(nèi)存通信，即在Kernel函數(shù)中存在2個(gè)層次的并行：在線程網(wǎng)格中的線程塊間并行和在線程塊中的線程間并行.

在Kernel函數(shù)的程序中，每個(gè)線程擁有自己的私有寄存器（register）和局部存儲(chǔ)器（local memory）；每個(gè)線程塊擁有一個(gè)共享存儲(chǔ)器（shared memory）；線程網(wǎng)格中所有線程都可以訪問(wèn)全局存儲(chǔ)器（global memory）.還有2種程序中所有線程都可以訪問(wèn)的只讀存儲(chǔ)器：常數(shù)存儲(chǔ)器（constant memory）和紋理存儲(chǔ)器（texture memory），它們分別為不同的應(yīng)用進(jìn)行優(yōu)化.其中，全局存儲(chǔ)器、常數(shù)存儲(chǔ)器和紋理存儲(chǔ)器中的值在一個(gè)內(nèi)核函數(shù)執(zhí)行完成后被繼續(xù)保持，可以被同一程序中的其他內(nèi)核函數(shù)調(diào)用.

2 有限元總剛矩陣

2.1 單元著色

為能夠并行計(jì)算單元?jiǎng)偠染仃嚽也划a(chǎn)生線程沖突，通過(guò)單元著色方法生成剛度矩陣.單元著色原則是：對(duì)于共享同一個(gè)節(jié)點(diǎn)的所有單元都著色為不同的顏色，即對(duì)于任一自由度中的單元沒(méi)有2個(gè)單元的顏色是相同的.

完成對(duì)結(jié)構(gòu)中的所有單元著色后，再進(jìn)行同一顏色的分組.

2.2 總剛矩陣的壓縮存儲(chǔ)

有限元法生成的總剛矩陣為大型稀疏矩陣，如果使用與稠密矩陣一樣的滿陣存儲(chǔ)法存儲(chǔ)，不僅存儲(chǔ)量和計(jì)算量大，而且會(huì)浪費(fèi)很多不必要的內(nèi)存空間，需采用壓縮存儲(chǔ)的方式存儲(chǔ)該稀疏矩陣，即只存儲(chǔ)稀疏矩陣中的非零元素.在進(jìn)行單元?jiǎng)偠扔?jì)算之前，先計(jì)算出單元中相互有貢獻(xiàn)的節(jié)點(diǎn)，對(duì)于相互沒(méi)有貢獻(xiàn)的節(jié)點(diǎn)，不給予存儲(chǔ)空間.壓縮存儲(chǔ)法中的行壓縮存儲(chǔ)（Compressed Sparse Row，CSR）法對(duì)矩陣的結(jié)構(gòu)沒(méi)有要求，而且以CSR格式存儲(chǔ)的稀疏矩陣能夠更好地滿足GPU的并行計(jì)算.有限元法生成的總剛矩陣通常是無(wú)規(guī)則的稀疏矩陣，因此采用CSR格式存儲(chǔ).

摘要：研究基于GPU的有限元求解中的總剛矩陣生成和線性方程組求解問(wèn)題.通過(guò)對(duì)單元著色和分組完成總剛矩陣的生成，并以行壓縮存儲(chǔ)（Compressed Sparse Row，CSR）格式存儲(chǔ)，用預(yù)處理共軛梯度法求解所生成的大規(guī)模線性稀疏方程組.在CUDA（Compute Unified Device Architecture）平臺(tái)上完成程序設(shè)計(jì)，并用GT430 GPU對(duì)彈性力學(xué)的平面問(wèn)題和空間問(wèn)題進(jìn)行試驗(yàn).結(jié)果表明，總剛矩陣生成和方程組求解分別得到最高11.7和8的計(jì)算加速比.

關(guān)鍵詞：GPU計(jì)算； 有限元法； 剛度矩陣； 預(yù)處理共軛梯度法

中圖分類號(hào)： TB115.7；TP311

文獻(xiàn)標(biāo)志碼：B

0 引 言

作為一種求解微分方程或積分方程的微分方法，有限元法[1]以其高度的適應(yīng)性，成為現(xiàn)代工程設(shè)計(jì)和結(jié)構(gòu)分析的重要方法之一，并在土木、水利、汽車、機(jī)械、航空航天、核工業(yè)和大地勘測(cè)等眾多領(lǐng)域應(yīng)用廣泛.隨著科學(xué)技術(shù)的不斷發(fā)展，工程問(wèn)題的規(guī)模和復(fù)雜程度相應(yīng)提高，也對(duì)有限元計(jì)算提出更大規(guī)模、更快速度的要求.有限元法的基本思想是“化整為零，積零為整”，與并行計(jì)算技術(shù)的“分而治之”的基本原則相協(xié)調(diào)，因此，對(duì)于大規(guī)模的有限元結(jié)構(gòu)分析，可將效率低下的串行有限元分析改進(jìn)為并行有限元分析.

自從2006年NVIDIA正式發(fā)布用于通用計(jì)算的統(tǒng)一計(jì)算架構(gòu)CUDA（Compute Unified Device Architecture）平臺(tái)[2]后，圖形處理器的體系架構(gòu)得到迅速發(fā)展和完善，性能得到大幅提高，使得GPU不僅能高效應(yīng)用于計(jì)算機(jī)圖形處理，而且其強(qiáng)大的計(jì)算能力也能很好地適用于高性能計(jì)算[3]，大大推動(dòng)基于GPU通用計(jì)算的研究，并廣泛應(yīng)用于醫(yī)學(xué)成像、生物信息學(xué)、計(jì)算結(jié)構(gòu)力學(xué)、計(jì)算流體力學(xué)、計(jì)算金融、地震勘探、地理信息系統(tǒng)以及電影和動(dòng)畫制作等領(lǐng)域.商業(yè)有限元分析軟件ANSYS和Abaqus中也應(yīng)用GPU技術(shù).[4]

目前，基于GPU的有限元分析在求解由有限元生成的稀疏線性方程組的研究中比較多，主要是由于在求解方程組時(shí)數(shù)據(jù)量大而且計(jì)算比較集中，便于并行運(yùn)算.NATHAN等[5]討論稀疏矩陣的數(shù)據(jù)結(jié)構(gòu)，并探討幾種有效的基于CUDA的高效的稀疏矩陣與向量相乘的方法；AIL等[6]探討針對(duì)多GPU的基于CUDA的快速共軛梯度法，并探討共軛梯度法中最耗時(shí)的稀疏矩陣與向量相乘的操作；JOLDES等[7]研究基于GPU的以六面體單元為主的混合網(wǎng)格及在有限元中尋求穩(wěn)定解的問(wèn)題，并通過(guò)CUDA實(shí)現(xiàn)基于非線性力學(xué)模型的自動(dòng)模擬神經(jīng)外科的過(guò)程；PAWEL等[8]探討基于GPU的三維有限元數(shù)值積分算法和計(jì)算方面的內(nèi)容；CECKA等[9]探討基于GPU的有限元法剛度矩陣組裝方法，評(píng)估每種方法的優(yōu)缺點(diǎn)；李熙銘[10]驗(yàn)證基于CUDA的復(fù)電阻率問(wèn)題，并詳細(xì)研究復(fù)共軛梯度法；胡耀國(guó)[11]運(yùn)用單元分組的方式計(jì)算得出有限元法中的總剛矩陣，并研究基于GPU的共軛梯度法，因當(dāng)時(shí)的資源限制，其在總剛矩陣生成部分的加速效果并不明顯.

本文利用GPU強(qiáng)大的并行計(jì)算能力和CPU的高效邏輯處理能力，將有限元法中計(jì)算量很大的單剛矩陣計(jì)算和總剛矩陣生成及線性方程組的求解交給GPU運(yùn)算，CPU負(fù)責(zé)相應(yīng)數(shù)據(jù)的前處理和后處理及整個(gè)分析過(guò)程中的邏輯關(guān)系.應(yīng)用CUDA平臺(tái)完成基于GPU的有限元分析程序，在本文的計(jì)算平臺(tái)上運(yùn)行，并用彈性力學(xué)平面問(wèn)題和空間問(wèn)題的有限元求解測(cè)試其運(yùn)行效果.在總剛矩陣生成部分，為避免計(jì)算時(shí)線程沖突，采用文獻(xiàn)[9]和[11-12]中都提到的對(duì)結(jié)構(gòu)單元進(jìn)行著色的方法.在求解大規(guī)模線性方程組時(shí)充分利用現(xiàn)有資源，在完全調(diào)用現(xiàn)有庫(kù)函數(shù)的情況下實(shí)現(xiàn)基于GPU的預(yù)條件共軛梯度法.

1 CUDA編程模型

NVIDIA的CUDA并行計(jì)算模型使用C語(yǔ)言作為開(kāi)發(fā)語(yǔ)言，同時(shí)也支持其他編程語(yǔ)言或應(yīng)用程序接口.一個(gè)完整的CUDA程序包含連續(xù)運(yùn)行在CPU端的主程序及用CUDA編寫的運(yùn)行在并行GPU設(shè)備上的內(nèi)核程序.由圖1可知，在CUDA架構(gòu)下，應(yīng)用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）組成.運(yùn)行流程：在CPU端準(zhǔn)備數(shù)據(jù)，然后傳到GPU端進(jìn)行并行計(jì)算，最后將計(jì)算好的數(shù)據(jù)再傳回CPU端.

CPU端的程序主要負(fù)責(zé)數(shù)據(jù)的準(zhǔn)備和一些邏輯運(yùn)算，其中Kernel函數(shù)是整個(gè)應(yīng)用程序的關(guān)鍵.Kernel函數(shù)以線程網(wǎng)格（Grid）的形式組織，每個(gè)線程網(wǎng)格由若干個(gè)線程塊（Block）組成，而每個(gè)線程塊又由若干個(gè)線程（Thread）組成.在執(zhí)行時(shí)，Kernel函數(shù)以線程塊為單位，各線程塊間并行執(zhí)行，不同線程塊間只能通過(guò)全局顯存進(jìn)行數(shù)據(jù)共享，同一線程塊內(nèi)的線程之間可以通過(guò)共享內(nèi)存通信，即在Kernel函數(shù)中存在2個(gè)層次的并行：在線程網(wǎng)格中的線程塊間并行和在線程塊中的線程間并行.

在Kernel函數(shù)的程序中，每個(gè)線程擁有自己的私有寄存器（register）和局部存儲(chǔ)器（local memory）；每個(gè)線程塊擁有一個(gè)共享存儲(chǔ)器（shared memory）；線程網(wǎng)格中所有線程都可以訪問(wèn)全局存儲(chǔ)器（global memory）.還有2種程序中所有線程都可以訪問(wèn)的只讀存儲(chǔ)器：常數(shù)存儲(chǔ)器（constant memory）和紋理存儲(chǔ)器（texture memory），它們分別為不同的應(yīng)用進(jìn)行優(yōu)化.其中，全局存儲(chǔ)器、常數(shù)存儲(chǔ)器和紋理存儲(chǔ)器中的值在一個(gè)內(nèi)核函數(shù)執(zhí)行完成后被繼續(xù)保持，可以被同一程序中的其他內(nèi)核函數(shù)調(diào)用.

2 有限元總剛矩陣

2.1 單元著色

為能夠并行計(jì)算單元?jiǎng)偠染仃嚽也划a(chǎn)生線程沖突，通過(guò)單元著色方法生成剛度矩陣.單元著色原則是：對(duì)于共享同一個(gè)節(jié)點(diǎn)的所有單元都著色為不同的顏色，即對(duì)于任一自由度中的單元沒(méi)有2個(gè)單元的顏色是相同的.

完成對(duì)結(jié)構(gòu)中的所有單元著色后，再進(jìn)行同一顏色的分組.

2.2 總剛矩陣的壓縮存儲(chǔ)

有限元法生成的總剛矩陣為大型稀疏矩陣，如果使用與稠密矩陣一樣的滿陣存儲(chǔ)法存儲(chǔ)，不僅存儲(chǔ)量和計(jì)算量大，而且會(huì)浪費(fèi)很多不必要的內(nèi)存空間，需采用壓縮存儲(chǔ)的方式存儲(chǔ)該稀疏矩陣，即只存儲(chǔ)稀疏矩陣中的非零元素.在進(jìn)行單元?jiǎng)偠扔?jì)算之前，先計(jì)算出單元中相互有貢獻(xiàn)的節(jié)點(diǎn)，對(duì)于相互沒(méi)有貢獻(xiàn)的節(jié)點(diǎn)，不給予存儲(chǔ)空間.壓縮存儲(chǔ)法中的行壓縮存儲(chǔ)（Compressed Sparse Row，CSR）法對(duì)矩陣的結(jié)構(gòu)沒(méi)有要求，而且以CSR格式存儲(chǔ)的稀疏矩陣能夠更好地滿足GPU的并行計(jì)算.有限元法生成的總剛矩陣通常是無(wú)規(guī)則的稀疏矩陣，因此采用CSR格式存儲(chǔ).

摘要：研究基于GPU的有限元求解中的總剛矩陣生成和線性方程組求解問(wèn)題.通過(guò)對(duì)單元著色和分組完成總剛矩陣的生成，并以行壓縮存儲(chǔ)（Compressed Sparse Row，CSR）格式存儲(chǔ)，用預(yù)處理共軛梯度法求解所生成的大規(guī)模線性稀疏方程組.在CUDA（Compute Unified Device Architecture）平臺(tái)上完成程序設(shè)計(jì)，并用GT430 GPU對(duì)彈性力學(xué)的平面問(wèn)題和空間問(wèn)題進(jìn)行試驗(yàn).結(jié)果表明，總剛矩陣生成和方程組求解分別得到最高11.7和8的計(jì)算加速比.

關(guān)鍵詞：GPU計(jì)算； 有限元法； 剛度矩陣； 預(yù)處理共軛梯度法

中圖分類號(hào)： TB115.7；TP311

文獻(xiàn)標(biāo)志碼：B

0 引 言

作為一種求解微分方程或積分方程的微分方法，有限元法[1]以其高度的適應(yīng)性，成為現(xiàn)代工程設(shè)計(jì)和結(jié)構(gòu)分析的重要方法之一，并在土木、水利、汽車、機(jī)械、航空航天、核工業(yè)和大地勘測(cè)等眾多領(lǐng)域應(yīng)用廣泛.隨著科學(xué)技術(shù)的不斷發(fā)展，工程問(wèn)題的規(guī)模和復(fù)雜程度相應(yīng)提高，也對(duì)有限元計(jì)算提出更大規(guī)模、更快速度的要求.有限元法的基本思想是“化整為零，積零為整”，與并行計(jì)算技術(shù)的“分而治之”的基本原則相協(xié)調(diào)，因此，對(duì)于大規(guī)模的有限元結(jié)構(gòu)分析，可將效率低下的串行有限元分析改進(jìn)為并行有限元分析.

自從2006年NVIDIA正式發(fā)布用于通用計(jì)算的統(tǒng)一計(jì)算架構(gòu)CUDA（Compute Unified Device Architecture）平臺(tái)[2]后，圖形處理器的體系架構(gòu)得到迅速發(fā)展和完善，性能得到大幅提高，使得GPU不僅能高效應(yīng)用于計(jì)算機(jī)圖形處理，而且其強(qiáng)大的計(jì)算能力也能很好地適用于高性能計(jì)算[3]，大大推動(dòng)基于GPU通用計(jì)算的研究，并廣泛應(yīng)用于醫(yī)學(xué)成像、生物信息學(xué)、計(jì)算結(jié)構(gòu)力學(xué)、計(jì)算流體力學(xué)、計(jì)算金融、地震勘探、地理信息系統(tǒng)以及電影和動(dòng)畫制作等領(lǐng)域.商業(yè)有限元分析軟件ANSYS和Abaqus中也應(yīng)用GPU技術(shù).[4]

目前，基于GPU的有限元分析在求解由有限元生成的稀疏線性方程組的研究中比較多，主要是由于在求解方程組時(shí)數(shù)據(jù)量大而且計(jì)算比較集中，便于并行運(yùn)算.NATHAN等[5]討論稀疏矩陣的數(shù)據(jù)結(jié)構(gòu)，并探討幾種有效的基于CUDA的高效的稀疏矩陣與向量相乘的方法；AIL等[6]探討針對(duì)多GPU的基于CUDA的快速共軛梯度法，并探討共軛梯度法中最耗時(shí)的稀疏矩陣與向量相乘的操作；JOLDES等[7]研究基于GPU的以六面體單元為主的混合網(wǎng)格及在有限元中尋求穩(wěn)定解的問(wèn)題，并通過(guò)CUDA實(shí)現(xiàn)基于非線性力學(xué)模型的自動(dòng)模擬神經(jīng)外科的過(guò)程；PAWEL等[8]探討基于GPU的三維有限元數(shù)值積分算法和計(jì)算方面的內(nèi)容；CECKA等[9]探討基于GPU的有限元法剛度矩陣組裝方法，評(píng)估每種方法的優(yōu)缺點(diǎn)；李熙銘[10]驗(yàn)證基于CUDA的復(fù)電阻率問(wèn)題，并詳細(xì)研究復(fù)共軛梯度法；胡耀國(guó)[11]運(yùn)用單元分組的方式計(jì)算得出有限元法中的總剛矩陣，并研究基于GPU的共軛梯度法，因當(dāng)時(shí)的資源限制，其在總剛矩陣生成部分的加速效果并不明顯.

本文利用GPU強(qiáng)大的并行計(jì)算能力和CPU的高效邏輯處理能力，將有限元法中計(jì)算量很大的單剛矩陣計(jì)算和總剛矩陣生成及線性方程組的求解交給GPU運(yùn)算，CPU負(fù)責(zé)相應(yīng)數(shù)據(jù)的前處理和后處理及整個(gè)分析過(guò)程中的邏輯關(guān)系.應(yīng)用CUDA平臺(tái)完成基于GPU的有限元分析程序，在本文的計(jì)算平臺(tái)上運(yùn)行，并用彈性力學(xué)平面問(wèn)題和空間問(wèn)題的有限元求解測(cè)試其運(yùn)行效果.在總剛矩陣生成部分，為避免計(jì)算時(shí)線程沖突，采用文獻(xiàn)[9]和[11-12]中都提到的對(duì)結(jié)構(gòu)單元進(jìn)行著色的方法.在求解大規(guī)模線性方程組時(shí)充分利用現(xiàn)有資源，在完全調(diào)用現(xiàn)有庫(kù)函數(shù)的情況下實(shí)現(xiàn)基于GPU的預(yù)條件共軛梯度法.

1 CUDA編程模型

NVIDIA的CUDA并行計(jì)算模型使用C語(yǔ)言作為開(kāi)發(fā)語(yǔ)言，同時(shí)也支持其他編程語(yǔ)言或應(yīng)用程序接口.一個(gè)完整的CUDA程序包含連續(xù)運(yùn)行在CPU端的主程序及用CUDA編寫的運(yùn)行在并行GPU設(shè)備上的內(nèi)核程序.由圖1可知，在CUDA架構(gòu)下，應(yīng)用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）組成.運(yùn)行流程：在CPU端準(zhǔn)備數(shù)據(jù)，然后傳到GPU端進(jìn)行并行計(jì)算，最后將計(jì)算好的數(shù)據(jù)再傳回CPU端.

CPU端的程序主要負(fù)責(zé)數(shù)據(jù)的準(zhǔn)備和一些邏輯運(yùn)算，其中Kernel函數(shù)是整個(gè)應(yīng)用程序的關(guān)鍵.Kernel函數(shù)以線程網(wǎng)格（Grid）的形式組織，每個(gè)線程網(wǎng)格由若干個(gè)線程塊（Block）組成，而每個(gè)線程塊又由若干個(gè)線程（Thread）組成.在執(zhí)行時(shí)，Kernel函數(shù)以線程塊為單位，各線程塊間并行執(zhí)行，不同線程塊間只能通過(guò)全局顯存進(jìn)行數(shù)據(jù)共享，同一線程塊內(nèi)的線程之間可以通過(guò)共享內(nèi)存通信，即在Kernel函數(shù)中存在2個(gè)層次的并行：在線程網(wǎng)格中的線程塊間并行和在線程塊中的線程間并行.

在Kernel函數(shù)的程序中，每個(gè)線程擁有自己的私有寄存器（register）和局部存儲(chǔ)器（local memory）；每個(gè)線程塊擁有一個(gè)共享存儲(chǔ)器（shared memory）；線程網(wǎng)格中所有線程都可以訪問(wèn)全局存儲(chǔ)器（global memory）.還有2種程序中所有線程都可以訪問(wèn)的只讀存儲(chǔ)器：常數(shù)存儲(chǔ)器（constant memory）和紋理存儲(chǔ)器（texture memory），它們分別為不同的應(yīng)用進(jìn)行優(yōu)化.其中，全局存儲(chǔ)器、常數(shù)存儲(chǔ)器和紋理存儲(chǔ)器中的值在一個(gè)內(nèi)核函數(shù)執(zhí)行完成后被繼續(xù)保持，可以被同一程序中的其他內(nèi)核函數(shù)調(diào)用.

2 有限元總剛矩陣

2.1 單元著色

為能夠并行計(jì)算單元?jiǎng)偠染仃嚽也划a(chǎn)生線程沖突，通過(guò)單元著色方法生成剛度矩陣.單元著色原則是：對(duì)于共享同一個(gè)節(jié)點(diǎn)的所有單元都著色為不同的顏色，即對(duì)于任一自由度中的單元沒(méi)有2個(gè)單元的顏色是相同的.

完成對(duì)結(jié)構(gòu)中的所有單元著色后，再進(jìn)行同一顏色的分組.

2.2 總剛矩陣的壓縮存儲(chǔ)

有限元法生成的總剛矩陣為大型稀疏矩陣，如果使用與稠密矩陣一樣的滿陣存儲(chǔ)法存儲(chǔ)，不僅存儲(chǔ)量和計(jì)算量大，而且會(huì)浪費(fèi)很多不必要的內(nèi)存空間，需采用壓縮存儲(chǔ)的方式存儲(chǔ)該稀疏矩陣，即只存儲(chǔ)稀疏矩陣中的非零元素.在進(jìn)行單元?jiǎng)偠扔?jì)算之前，先計(jì)算出單元中相互有貢獻(xiàn)的節(jié)點(diǎn)，對(duì)于相互沒(méi)有貢獻(xiàn)的節(jié)點(diǎn)，不給予存儲(chǔ)空間.壓縮存儲(chǔ)法中的行壓縮存儲(chǔ)（Compressed Sparse Row，CSR）法對(duì)矩陣的結(jié)構(gòu)沒(méi)有要求，而且以CSR格式存儲(chǔ)的稀疏矩陣能夠更好地滿足GPU的并行計(jì)算.有限元法生成的總剛矩陣通常是無(wú)規(guī)則的稀疏矩陣，因此采用CSR格式存儲(chǔ).