有限元GPU加速計算的實現方法

張健飛　沈德飛

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

關鍵詞：GPU計算； 有限元法； 剛度矩陣； 預處理共軛梯度法

中圖分類號： TB115.7；TP311

文獻標志碼：B

0 引 言

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

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

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

本文利用GPU強大的并行計算能力和CPU的高效邏輯處理能力，將有限元法中計算量很大的單剛矩陣計算和總剛矩陣生成及線性方程組的求解交給GPU運算，CPU負責相應數據的前處理和后處理及整個分析過程中的邏輯關系.應用CUDA平臺完成基于GPU的有限元分析程序，在本文的計算平臺上運行，并用彈性力學平面問題和空間問題的有限元求解測試其運行效果.在總剛矩陣生成部分，為避免計算時線程沖突，采用文獻[9]和[11-12]中都提到的對結構單元進行著色的方法.在求解大規模線性方程組時充分利用現有資源，在完全調用現有庫函數的情況下實現基于GPU的預條件共軛梯度法.

1 CUDA編程模型

NVIDIA的CUDA并行計算模型使用C語言作為開發語言，同時也支持其他編程語言或應用程序接口.一個完整的CUDA程序包含連續運行在CPU端的主程序及用CUDA編寫的運行在并行GPU設備上的內核程序.由圖1可知，在CUDA架構下，應用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）組成.運行流程：在CPU端準備數據，然后傳到GPU端進行并行計算，最后將計算好的數據再傳回CPU端.

CPU端的程序主要負責數據的準備和一些邏輯運算，其中Kernel函數是整個應用程序的關鍵.Kernel函數以線程網格（Grid）的形式組織，每個線程網格由若干個線程塊（Block）組成，而每個線程塊又由若干個線程（Thread）組成.在執行時，Kernel函數以線程塊為單位，各線程塊間并行執行，不同線程塊間只能通過全局顯存進行數據共享，同一線程塊內的線程之間可以通過共享內存通信，即在Kernel函數中存在2個層次的并行：在線程網格中的線程塊間并行和在線程塊中的線程間并行.

在Kernel函數的程序中，每個線程擁有自己的私有寄存器（register）和局部存儲器（local memory）；每個線程塊擁有一個共享存儲器（shared memory）；線程網格中所有線程都可以訪問全局存儲器（global memory）.還有2種程序中所有線程都可以訪問的只讀存儲器：常數存儲器（constant memory）和紋理存儲器（texture memory），它們分別為不同的應用進行優化.其中，全局存儲器、常數存儲器和紋理存儲器中的值在一個內核函數執行完成后被繼續保持，可以被同一程序中的其他內核函數調用.

2 有限元總剛矩陣

2.1 單元著色

為能夠并行計算單元剛度矩陣且不產生線程沖突，通過單元著色方法生成剛度矩陣.單元著色原則是：對于共享同一個節點的所有單元都著色為不同的顏色，即對于任一自由度中的單元沒有2個單元的顏色是相同的.

完成對結構中的所有單元著色后，再進行同一顏色的分組.

2.2 總剛矩陣的壓縮存儲

有限元法生成的總剛矩陣為大型稀疏矩陣，如果使用與稠密矩陣一樣的滿陣存儲法存儲，不僅存儲量和計算量大，而且會浪費很多不必要的內存空間，需采用壓縮存儲的方式存儲該稀疏矩陣，即只存儲稀疏矩陣中的非零元素.在進行單元剛度計算之前，先計算出單元中相互有貢獻的節點，對于相互沒有貢獻的節點，不給予存儲空間.壓縮存儲法中的行壓縮存儲（Compressed Sparse Row，CSR）法對矩陣的結構沒有要求，而且以CSR格式存儲的稀疏矩陣能夠更好地滿足GPU的并行計算.有限元法生成的總剛矩陣通常是無規則的稀疏矩陣，因此采用CSR格式存儲.

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

關鍵詞：GPU計算； 有限元法； 剛度矩陣； 預處理共軛梯度法

中圖分類號： TB115.7；TP311

文獻標志碼：B

0 引 言

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

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

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

本文利用GPU強大的并行計算能力和CPU的高效邏輯處理能力，將有限元法中計算量很大的單剛矩陣計算和總剛矩陣生成及線性方程組的求解交給GPU運算，CPU負責相應數據的前處理和后處理及整個分析過程中的邏輯關系.應用CUDA平臺完成基于GPU的有限元分析程序，在本文的計算平臺上運行，并用彈性力學平面問題和空間問題的有限元求解測試其運行效果.在總剛矩陣生成部分，為避免計算時線程沖突，采用文獻[9]和[11-12]中都提到的對結構單元進行著色的方法.在求解大規模線性方程組時充分利用現有資源，在完全調用現有庫函數的情況下實現基于GPU的預條件共軛梯度法.

1 CUDA編程模型

NVIDIA的CUDA并行計算模型使用C語言作為開發語言，同時也支持其他編程語言或應用程序接口.一個完整的CUDA程序包含連續運行在CPU端的主程序及用CUDA編寫的運行在并行GPU設備上的內核程序.由圖1可知，在CUDA架構下，應用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）組成.運行流程：在CPU端準備數據，然后傳到GPU端進行并行計算，最后將計算好的數據再傳回CPU端.

CPU端的程序主要負責數據的準備和一些邏輯運算，其中Kernel函數是整個應用程序的關鍵.Kernel函數以線程網格（Grid）的形式組織，每個線程網格由若干個線程塊（Block）組成，而每個線程塊又由若干個線程（Thread）組成.在執行時，Kernel函數以線程塊為單位，各線程塊間并行執行，不同線程塊間只能通過全局顯存進行數據共享，同一線程塊內的線程之間可以通過共享內存通信，即在Kernel函數中存在2個層次的并行：在線程網格中的線程塊間并行和在線程塊中的線程間并行.

在Kernel函數的程序中，每個線程擁有自己的私有寄存器（register）和局部存儲器（local memory）；每個線程塊擁有一個共享存儲器（shared memory）；線程網格中所有線程都可以訪問全局存儲器（global memory）.還有2種程序中所有線程都可以訪問的只讀存儲器：常數存儲器（constant memory）和紋理存儲器（texture memory），它們分別為不同的應用進行優化.其中，全局存儲器、常數存儲器和紋理存儲器中的值在一個內核函數執行完成后被繼續保持，可以被同一程序中的其他內核函數調用.

2 有限元總剛矩陣

2.1 單元著色

為能夠并行計算單元剛度矩陣且不產生線程沖突，通過單元著色方法生成剛度矩陣.單元著色原則是：對于共享同一個節點的所有單元都著色為不同的顏色，即對于任一自由度中的單元沒有2個單元的顏色是相同的.

完成對結構中的所有單元著色后，再進行同一顏色的分組.

2.2 總剛矩陣的壓縮存儲

有限元法生成的總剛矩陣為大型稀疏矩陣，如果使用與稠密矩陣一樣的滿陣存儲法存儲，不僅存儲量和計算量大，而且會浪費很多不必要的內存空間，需采用壓縮存儲的方式存儲該稀疏矩陣，即只存儲稀疏矩陣中的非零元素.在進行單元剛度計算之前，先計算出單元中相互有貢獻的節點，對于相互沒有貢獻的節點，不給予存儲空間.壓縮存儲法中的行壓縮存儲（Compressed Sparse Row，CSR）法對矩陣的結構沒有要求，而且以CSR格式存儲的稀疏矩陣能夠更好地滿足GPU的并行計算.有限元法生成的總剛矩陣通常是無規則的稀疏矩陣，因此采用CSR格式存儲.

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

關鍵詞：GPU計算； 有限元法； 剛度矩陣； 預處理共軛梯度法

中圖分類號： TB115.7；TP311

文獻標志碼：B

0 引 言

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

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

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

本文利用GPU強大的并行計算能力和CPU的高效邏輯處理能力，將有限元法中計算量很大的單剛矩陣計算和總剛矩陣生成及線性方程組的求解交給GPU運算，CPU負責相應數據的前處理和后處理及整個分析過程中的邏輯關系.應用CUDA平臺完成基于GPU的有限元分析程序，在本文的計算平臺上運行，并用彈性力學平面問題和空間問題的有限元求解測試其運行效果.在總剛矩陣生成部分，為避免計算時線程沖突，采用文獻[9]和[11-12]中都提到的對結構單元進行著色的方法.在求解大規模線性方程組時充分利用現有資源，在完全調用現有庫函數的情況下實現基于GPU的預條件共軛梯度法.

1 CUDA編程模型

NVIDIA的CUDA并行計算模型使用C語言作為開發語言，同時也支持其他編程語言或應用程序接口.一個完整的CUDA程序包含連續運行在CPU端的主程序及用CUDA編寫的運行在并行GPU設備上的內核程序.由圖1可知，在CUDA架構下，應用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）組成.運行流程：在CPU端準備數據，然后傳到GPU端進行并行計算，最后將計算好的數據再傳回CPU端.

CPU端的程序主要負責數據的準備和一些邏輯運算，其中Kernel函數是整個應用程序的關鍵.Kernel函數以線程網格（Grid）的形式組織，每個線程網格由若干個線程塊（Block）組成，而每個線程塊又由若干個線程（Thread）組成.在執行時，Kernel函數以線程塊為單位，各線程塊間并行執行，不同線程塊間只能通過全局顯存進行數據共享，同一線程塊內的線程之間可以通過共享內存通信，即在Kernel函數中存在2個層次的并行：在線程網格中的線程塊間并行和在線程塊中的線程間并行.

在Kernel函數的程序中，每個線程擁有自己的私有寄存器（register）和局部存儲器（local memory）；每個線程塊擁有一個共享存儲器（shared memory）；線程網格中所有線程都可以訪問全局存儲器（global memory）.還有2種程序中所有線程都可以訪問的只讀存儲器：常數存儲器（constant memory）和紋理存儲器（texture memory），它們分別為不同的應用進行優化.其中，全局存儲器、常數存儲器和紋理存儲器中的值在一個內核函數執行完成后被繼續保持，可以被同一程序中的其他內核函數調用.

2 有限元總剛矩陣

2.1 單元著色

為能夠并行計算單元剛度矩陣且不產生線程沖突，通過單元著色方法生成剛度矩陣.單元著色原則是：對于共享同一個節點的所有單元都著色為不同的顏色，即對于任一自由度中的單元沒有2個單元的顏色是相同的.

完成對結構中的所有單元著色后，再進行同一顏色的分組.

2.2 總剛矩陣的壓縮存儲

有限元法生成的總剛矩陣為大型稀疏矩陣，如果使用與稠密矩陣一樣的滿陣存儲法存儲，不僅存儲量和計算量大，而且會浪費很多不必要的內存空間，需采用壓縮存儲的方式存儲該稀疏矩陣，即只存儲稀疏矩陣中的非零元素.在進行單元剛度計算之前，先計算出單元中相互有貢獻的節點，對于相互沒有貢獻的節點，不給予存儲空間.壓縮存儲法中的行壓縮存儲（Compressed Sparse Row，CSR）法對矩陣的結構沒有要求，而且以CSR格式存儲的稀疏矩陣能夠更好地滿足GPU的并行計算.有限元法生成的總剛矩陣通常是無規則的稀疏矩陣，因此采用CSR格式存儲.