基于CUDA的加速MATLAB計算研究

摘 要:介紹了NVIDIA公司新的編程框架CUDA的特點以及CUDA加速MATLAB的方法，測試了CUDA加速巖土工程中常用的算法如矩陣計算、快速傅里葉變換、支持向量機。隨后分析了數據規模、算法復雜性與加速效果的關系，指出了基于CUDA的MATLAB加速計算的應用前景。測試結果表明，CUDA方式相對傳統計算方式的最好加速效果分別達到了22.39倍、46.88倍、51.32倍，證明了CUDA加速計算的有效性。

關鍵詞:統一計算設備架構; MATLAB; 加速計算

中圖分類號:TP311文獻標志碼:A

文章編號:1001-3695(2010)06-2140-04

doi:10.3969/j.issn.1001-3695.2010.06.041

Model of accelerating MATLAB computation based on CUDA

LIU Shaobo， LIU Minggui， ZHANG Guohua

(State Key Laboratory of Geomechanics Geotechnical Engineering， Institute of Rock Soil Mechanics， Chinese Academy of Sciences， Wuhan 430071， China)

Abstract:The character of NVIDIA Corporation’s new program framework， described CUDA and its usage of accelerating MATLAB computation basedon CUDA firstly. Then tested the algorithms of matrix computation， fast fourier transform (FFT)， support vector machine (SVM) in detail while used those algorithms widely in geotechnical engineering. Analyzed the relationship between acceleration performance and data size or algorithm complexity. Finally， discussed the prospect of accelerating MATLAB computation basedon CUDA.

The test shows that these algorithms computed by CUDA are 22.39 times， 46.88 times and 51.32 times faster than those by traditional ways.

Key words:compute unified device architecture(CUDA); MATLAB; accelerating computation

由于邊坡、隧道等巖土工程的地質條件限制，人們很難深入了解巖土工程內部的復雜結構。對于這類問題的研究，常用有限元[1]、神經網絡[2~4]等方法，模擬或預測環境內部的一些參數變量如溫度、應力等。許多研究人員利用MATLAB軟件的強大功能，在上述研究中取得了不錯的成績，給科研工作帶來了很大便利。

MATLAB雖然在易用性等方面有著許多優點，但是其計算速度慢、效率低，尤其是在大規模計算方面常常需要耗費較多時間。為了解決MATLAB計算效率問題，目前研究人員大多數嘗試下述兩種方式:a)通過購買昂貴的服務器、工作站等設備，來加快MATLAB計算速度。這種方式雖然能夠很好地解決計算速度問題，但是給科研人員帶來了經濟上的負擔。b)移植MATLAB算法，通過高級語言如C++等重新編制程序，提供執行效率。這種方式雖然能夠顯著提高MATLAB計算速度，但是對研究人員的編程功底要求較高，而且有些算法本身很復雜，C++等高級語言沒有提供相應的庫函數，這需要研究人員從底層開始實現算法，費時費力，加大了科研工作量。

近幾年來，計算機顯卡核心單元GPU(graphics processing unit)在浮點運算性能方面取得了突飛猛進的發展(如NVIDIA公司推出的GeForce GTX200系列高端顯卡，浮點運算速度逼近每秒萬億次左右，相當于過去小型機的處理能力)，GPU以其遠遠超過CPU的浮點運算性能、高內存帶寬、高性價比等優點，越來越受到廣大科研人員的重視，并在分子生物模擬、地震數據處理、超級計算機、金融模型計算等[5~8]方面取得了不錯的加速效果，一些計算的加速效果達到100多倍。雖然GPU計算有如此好的加速效果和高速浮點運算性能等優點，但是在巖土工程的相關計算中還沒有引入該方法?；谶@一目的，本文介紹了NVIDIA公司的GPU通用計算框架CUDA，以及利用CUDA加速計算MATLAB程序的原理，通過規模矩陣計算、快速傅里葉變換(FFT)、支持向量機(SVM)等巖土工程計算中的常見算法，分析了CUDA加速MATLAB程序的情況。

1 CUDA框架特點

1.1 CUDA簡介

為了充分發揮GPU的高速浮點運算等方面性能，NVIDIA公司研究人員針對該公司的GPU設計出了一種新的計算框架CUDA。該框架是一個軟硬件協同的完整并行計算解決方案，提供了硬件的直接訪問接口以及高性能計算指令，可以使GPU與CPU協同工作，合理管理GPU設備、存儲器、流、事件等，充分利用GPU高速內存帶寬，從而實現復雜的問題加速計算，尤其是大規模的并行計算問題。CUDA開發工具主要包括CUDA驅動程序、CUDA運行時庫、CUDA庫三個組件。CUDA軟件的核心層次結構如圖1所示[9]。

1.2 MATLAB在CUDA框架下的的使用

由于MATLAB的M語言編程時與外部環境的數據和程序的交互非常有限，MATLAB提供了功能更強大的MEX腳本文件。利用MEX文件，可以調用C、Fortran等語言、輸入或輸出數據、與其他軟件建立客戶/服務器通信、直接控制硬件等功能[10]。

NVIDIA公司針對MATLAB開發了一個新的插件，該插件包含一個名為nvmex腳本和一個簡化CUDA配置的可選配置文件。通過nvmex可以將CUDA源程序編譯為MATLAB的MEX腳本文件，使MATLAB軟件方便地調用CUDA函數庫，方便數據計算[11，12]。

一個基本的CUDA mex文件數據計算操作包含以下幾個步驟:a)在GPU中分配合適的數據存儲空間。b)將數據從主機內存移動到GPU內存中，GPU開始分配線程塊等準備工作，并根據CUDA代碼調用相關函數庫，做相應的數據計算操作。GPU計算完成后，將數據返回主機內存中。c)釋放GPU中的數據存儲空間，完成整個計算流程[10]。其流程如圖2所示。

2 MATLAB加速計算在巖土工程中的應用

在巖土工程中，常常使用矩陣計算、快速傅里葉變換(FFT)、支持向量機等算法來解決巖土工程中的計算問題。例如，矩陣計算可以解決巖土力學[13]、滲流有限元分析[14]等方面問題，快速傅里葉變換模型可以解決動荷載作用下土的動力響應[15，16]問題，支持向量機模型可以解決巖體工程分級[17]、地下工程可靠性分析[18]、邊坡穩定性[19]等方面問題。

下面將從測試平臺、測試算法、測試結果方面，詳細探討MATLAB對上述算法的加速效果。

2.1 測試平臺環境選擇

本文軟件測試平臺選用的是Windows XP SP2操作系統，支持CUDA插件的MATLAB 7.3.0 (R2006b)，CUDA 2.1版本庫。硬件測試平臺是雙核的Intel Core2Duo E4400 2.0 GHz CPU，影馳9600GT顯卡(其核心是NVIDIA GeForce 9600 GT)，DDR2 800金士頓1 GB內存。

2.2 測試算法

CUDA加速MATLAB的主要原理是將流程控制的操作如串行代碼放在CPU中操作，而一些需要大規模并行計算的操作代碼放在GPU中處理，從而縮短了整個計算操作的時間。MATLAB中數據的計時方式采用tic、toc函數來標記起止時間。CUDA計算數據的流程如圖2所示。下面詳細介紹這些算法的實現。

2.2.1 矩陣計算

矩陣計算模型采用的是常見的矩陣乘法，通過不同維度的矩陣相乘，比較出在不同維度下的加速效果。為了研究的方便和矩陣數據一致性，這里的矩陣都采用方陣，矩陣數據是由生成的隨機數填充，并保存在txt文件中。

傳統MATLAB的實現方式是先在主機上分配好矩陣的存儲空間，將txt文件數據讀入到數組中，然后調用矩陣乘法命令，將計算的結果保存到指定的結果數組中，完成整個計算。

基于CUDA的實現方式主要是利用CUDA自帶的CUBLAS庫進行矩陣計算。其步驟是主機先將txt文件數據讀入主機內存;然后GPU調用cublasInit函數初始化GPU，并調用cublasAlloc函數分配矩陣的GPU存儲空間，將主機內存中的矩陣數據讀入GPU中;最后GPU調用cublasSgemm函數進行矩陣乘法計算，并通過cublasGetVector函數將計算結果返回主機，調用cublasFree函數釋放GPU上的矩陣存儲空間，完成整個計算工作[11]。

2.2.2 快速傅里葉變換(FFT)

FFT計算方式與矩陣計算類似，也是先生成隨機的FFT原始數據，并將數據保存在txt文件中，保證測試數據的統一性。

傳統MATLAB的實現方式是將txt文件的數據讀入到指定的數組中，然后調用fft函數，并返回指定結果，完成整個計算。

基于CUDA的實現方式是利用CUDA自帶的CUFFT庫。實現步驟是先將txt文件的數據讀入主機內存;然后GPU調用cudaMalloc分配GPU上的存儲空間，調用cudaMemcpy將數據從主機內存復制到GPU內存，再調用cufftPlan1d建立一維FFT plan;最后調用cufftExecC2C函數完成FFT變換，將計算結果返回主機內存，釋放plan以及GPU上的數據存儲空間，完成整個計算工作[11]。

2.2.3 支持向量機(SVM)

SVM算法分為訓練、預測兩部分。其計算方式與上述算法類似，先生成一定數量的訓練樣本和測試樣本數據，然后保存在txt文件中，由兩種算法分別對測試樣本進行訓練、預測。

傳統的MATLAB的實現方式是將txt文件的數據讀入指定的數組中，然后設定好訓練參數，調用svmtrain函數進行訓練并得到訓練模型，根據訓練模型和測試數據調用svmclassify函數進行預測，并返回相應的預測結果。

基于CUDA的實現方式主要是利用參考文獻[20]的cuSVM庫進行實驗。實現步驟是先讀取txt文件數據，然后設置好訓練參數，調用cuSVMTrain函數對訓練樣本進行訓練并生產訓練模型，根據訓練模型和測試數據調用cuSVMPredict函數進行預測，返回相應預測結果，完成整個計算。

2.3 測試結果

2.3.1 矩陣計算

矩陣測試的數據采用行列數相同的方陣，分別采用維度為500×500、800×800、1 000×1 000、3 000×3 000、5 000×5 000的矩陣(方陣)相乘，計算時間如表1所示。表1中的加速倍數是由傳統方式計算時間與CUDA方式計算時間的比值得到的，本文中的其他算法表格的加速倍數均采用類似方式計算。從表1中可以看出，當矩陣規模為500×500時，CUDA方式時間為0.026 3 s，而傳統方式時間為0.015 6 s，可以看出CUDA方式不僅沒有加速作用，反而落后于傳統的計算方式約1倍的時間(0.59倍);當矩陣規模增加到800×800、1 000×1 000時，CUDA方式略有增加，但是效果不是很明顯，只有1.28倍和1.58倍加速;當矩陣規模增到3 000×3 000時，CUDA方式取得了明顯的加速效果，達到了5.16倍，比傳統方式節約了29.734-5.761=23.973 s的時間;當矩陣規模達到5 000×5 000時，CUDA方式取得了不錯的加速效果，達到了22.39倍，遠小于傳統方式計算所花的時間。

表1 矩陣計算結果

矩陣規模傳統方式時間/sCUDA方式時間/s加速倍數

500×5000.015 6000.026 3000.59

800×8000.578 0000.453 0001.28

1 000×1 0001.188 0000.752 0001.58

3 000×3 00029.734 0005.761 0005.16

5000×5000137.844 0006.156 00022.39

2.3.2 快速傅里葉變換(FFT)

FFT測試的數據分別采用維度為1 024×8、1 024×32、1 024×128、1 024×512、1 024×1 024的一維FFT變換。這里FFT數據規模均采用2的整數次冪，主要是為了方便FFT計算。計算時間如表2所示。

表2 FFT計算結果

FFT規模傳統方式時間/sCUDA方式時間/s加速倍數

1 024×80.016 0000.0121.33

1 024×320.094 0000.0234.09

1 024×1280.453 0000.04111.78

1 024×5122.015 0000.05735.36

1 024×1 0243.844 0000.08246.88

從表2可看出，當FFT規模為1 024×8時，CUDA方式與傳統方式計算相差不大，時間相差0.016-0.012=0.004 s，加速倍數為1.33;當FFT規模為1 024×32、1 024×128時，CUDA方式取得了明顯的加速效果，分別達到了4.09倍、11.78倍;當FFT規模為1 024×512、1 024×1 024時，CUDA方式取得了驚人的加速效果，分別達到了35.36倍、46.88倍，遠遠優于傳統計算方式。

2.3.3 支持向量機(SVM)

SVM采用的樣本特征值數為100，樣本數分別3 000、5 000、8 000、10 000、30 000，訓練參數采用的是CSVC核心函數，C值為100，γ值為0.5。預測計算采用的測試數據是訓練樣本中的部分數據，根據訓練模型進行預測。訓練結果和預測結果分別如表3和4所示。

表3 SVM訓練結果

樣本規模傳統方式時間/sCUDA方式時間/s加速倍數

3 00065.728 00018.735 0003.51

5 000112.997 00020.178 0005.60

8 000220.542 00023.972 0009.20

10 000357.984 00029.832 00012.01

30 000523.538 00032.614 00016.05

表4 SVM預測結果

測試規模

傳統方式

時間/s精度/%

CUDA方式

時間/s精度/%

加速倍數

2 0005.796 00092.740.386 00092.8315.02

4 0008.652 00093.270.412 00093.2121.00

6 00013.095 00093.860.485 00093.8427.01

8 00020.026 00094.760.572 00094.7935.01

10 00031.721 00095.590.618 00095.5651.32

從表3可以看出，當樣本規模為3 000、5 000、8 000時，CUDA方式雖然有一定的加速效果，但是效果不是很突出;當樣本規模達到10 000、30 000時，CUDA方式取得了明顯的加速效果。從表4可以看出，兩種計算方式的預測精度比較接近，但是計算時間差別較大。當測試規模為2 000時，CUDA方式就取得了15.02倍的明顯加速效果;當測試規模不斷加大時，CUDA方式加速效果越來越明顯，在測試規模達到10 000時，加速效果達到了驚人的51.32倍!

2.3.4 測試結論

分析上述測試結果，可以得出以下結論:

a)數據規模對加速效果的影響。從表1~4的測試結果簡單分析中可以看出，同樣的算法，當計算規模較小時，CUDA方式加速效果不太明顯，甚至出現“減速”效果，如矩陣計算規模為500×500，傳統計算時間為0.015 6 s，CUDA方式計算時間為0.026 3 s，落后于傳統方式計算;當規模較大時，CUDA加速效果顯著，產生了十幾倍到幾十倍的加速，而且規模越大，加速效果越明顯。

產生該現象的原因是，主機(host)內存數據與GPU內存數據交互(讀入、返回)需要一定的時間開銷，規模不大時這部分開銷對最終的計算時間有很大影響，只有當數據規模很大時，數據交互時間所占比例較小，GPU計算時間足以抵過這部分開銷，CUDA加速的效果才比較突出。

b)算法復雜度的影響。表5總結了各個算法的加速效果對比。從加速倍數的變化效果上看，矩陣規模從500×500擴大到5 000×5 000時，加速倍數從0.59擴大到22.39，即加速效果擴大了22.39/0.59≈37.95倍;FFT規模從1024×8擴大到1024×1024時，加速倍數從1.33擴大到46.88，即加速效果擴大了46.88/1.33≈35.25倍;SVM訓練算法規模從3 000擴大到30 000時，加速倍數從3.51擴大到16.05，即加速效果擴大了16.05/3.51≈4.57倍;SVM預測規模從2 000擴大到10 000時，加速倍數從15.02擴大到51.32，即加速效果擴大了51.32/15.02≈3.42倍。

表5 各算法加速效果

矩陣規模加速倍數FFT規模加速倍數SVM樣本規模加速倍數SVM測試規模加速倍數

500×5000.591024×81.3330003.51200015.02

800×8001.281024×324.0950005.60400021.00

1000×10001.581024×12811.7880009.20600027.01

3000×30005.161024×51235.361000012.01800035.01

5000×500022.391024×102446.883000016.051000051.32

從這里可以看出，矩陣計算和FFT計算加速效果非常明顯，而SVM訓練計算和預測計算效果略差。產生該現象的原因是矩陣乘法和FFT算法相對簡單，沒有復雜的控制;而SVM訓練算法和預測算法比較復雜，尤其是SVM訓練算法有很多分支控制，而且GPU在分支控制等指令方面能力遠不如CPU，大大削弱了GPU的大規模并行處理能力，導致加速效果變化不明顯。

c)CUDA加速的擴展性問題。圖3~6列出了上述算法的時間對比曲線圖。

從圖3可以看出，當矩陣規模在1 000×1 000以內時，傳統方式和CUDA方式的曲線斜率變化都不大;當矩陣規模從1 000×1 000擴展到3 000×3 000時，傳統方式曲線斜率變化非常明顯，而CUDA方式曲線斜率仍變化不大;當矩陣規模從3 000×3 000擴展到5 000×5 000時，傳統方式曲線斜率急劇增加，而CUDA方式曲線變化仍不明顯。同樣地，從圖4~6中可以看出，隨著計算規模的加大，傳統方式曲線斜率的變化都大于CUDA曲線斜率的變化。

曲線斜率的大小，在一定程度上可以反應出計算規模與計算時間的關系，即在相同的計算規模下，曲線斜率越大，說明該計算方式耗費的時間越多。當曲線斜率大到一定程度時，計算規模稍微擴大，也能導致計算時間的急劇增加，這時計算規模與耗費時間的性價比極低，形成了計算瓶頸，擴展性較差。

根據上述分析，可以看出傳統方式的擴展性不如CUDA方式，傳統方式更容易形成計算瓶頸。在巖土工程計算規模和復雜度不斷擴大的趨勢下，傳統方式無法應對這一問題，而CUDA方式則更有優勢。

從上面的分析可以看出，CUDA方式加速MATLAB計算比傳統方式更有優勢，可以有效地應對巖土工程計算規模日益擴大的問題，同時為了充分發揮CUDA的加速效果，應盡量選擇數據計算規模較大、沒有復雜控制的算法才能取得較好的加速效果，而且規模越大，算法越簡單，加速效果越理想。

3 結束語

隨著科學技術的不斷提高，巖土工程方向的建設規模越來越大，一些相關計算工作也越來越復雜，MATLAB的計算效率嚴重影響了科研工作的進展。如何加快MATLAB計算速度，是科研人員研究的一項重要內容。為了加速MATLAB計算速度，本文利用NVIDIA公司開發的GPU通用計算框架CUDA，并結合巖土工程計算領域中常用的矩陣計算、快速傅里葉變換、支持向量機算法進行測試。測試結果顯示，這些算法的最好加速效果分別達到了22.39倍、46.88倍、51.32倍，取得了不錯的加速效果，證明了CUDA加速MATLAB計算模式的有效性以及可行性。本文提出了基于CUDA的加速MATLAB計算模式，充分利用GPU強大的浮點運算功能和并行處理能力，能夠有效地加速MATLAB計算，該方案的低成本、高效率的優勢，可以有效地解決過去MATLAB加速計算問題。隨著巖土工程研究工作的不斷深入，科研計算工作的不斷增加，基于CUDA的加速MATLAB計算模型將會發揮越來越重要的作用。

參考文獻:

[1]葛修潤， 豐定祥. 巖體工程中流變問題的有限元分析[J]. 巖土力學， 1980(3): 15-20.

[2]馮夏庭， 刁心宏. 智能巖石力學(1)——導論[J]. 巖石力學與工程學報， 1999，18(2): 222-226.

[3]馮夏庭. 智能巖石力學導論[M]. 北京: 科學出版社，2000.

[4]劉明貴， 岳向紅， 楊永波，等. 基于Sym小波和BP神經網絡的基樁缺陷智能化識別[J]. 巖石力學與工程學報， 2007，26(增1):3484-3488.

[5]UIUC VMD Research Group. GPU acceleration of molecular modeling applications[EB/OL].(2007-09-11)[2009-11-02]. http://www.ks.uiuc.edu/Research/gpu/.

[6]Headwave. NVIDIA unveils CUDATM: the GPU computing revolution begins[EB/OL]. (2006-11-08)[2009-11-02]. http://www.nvidia.cn/object/headwave_geophysical_data_cn.html.

[7]ENDO T， MATSUOKA S. Massive supercomputing coping with heterogeneity of modern accelerators[C]//Proc of IEEE International Symposium on Parallel and Distributed Processing Symposium. Miami: IEEE Press，2008: 1-10.

[8]SCHMERKEN I. Wall street accelerates options analysis with GPU technology[EB/OL]. (2008-11-07)[2009-11-02]. http://wallstreetandtech.com/technologyriskmanagement/showArticle.jhtml?articleID=215801950.

[9]NVIDIA. CUDA reference manual[EB/OL]. (2007-06-23)[2009-11-02]. http://developer.download.nvidia.com/compute/cuda/2_1/toolkit/docs/CudaReferenceManual_2.1.pdf.

[10]董長虹. MATLAB接口技術與應用[M]. 北京: 國防工業出版社，2004:105-121.

[11]NVIDIA. MATLAB Acceleration using CUDAenabled GPUs[EB/OL]. (2007-07-05)[2009-11-02].http://www.nvidia.com/object/MATLAB_acceleration.html.

[12]NVIDIA Corporation. Accelerating MATLAB with CUDATM Using mEX files[EB/OL]. (2007-09-11)[2009-11-02].http://developer.download.nvidia.com/compute/cuda/1_0/Accelerating%

20MATLAB%20with%20CUDA.pdf.

[13]趙明華，李微哲，楊明輝，等. 成層地基中傾斜偏心荷載下單樁計算分析[J]. 巖土力學， 2007，28(4): 670-674.

[14]董平川，韓德金，牛彥良，等. 油藏多相滲流的面向對象有限元程序設計[J]. 巖土力學， 2009，30(4): 1115-1121，1130.

[15]姚海林，盧正，羅海寧，等. 交通荷載作用下Kelvin地基上不平整路面的動力響應分析[J]. 巖土力學， 2009，30(4):890-896.

[16]徐斌，陸建飛，王建華，等. 移動荷載作用下層狀飽和土的動力響應[J]. 巖土力學， 2008，29(12): 3186-3192.

[17]趙洪波，馮夏庭，尹順德.基于支持向量機的巖體工程分級[J]. 巖土力學，2002，23(6):698-701.

[18]趙洪波，茹忠亮，張士科. SVM在地下工程可靠性分析中的應用[J]. 巖土力學， 2009，30(2):526-530.

[19]馬文濤. 基于PSO和LSSVM的邊坡穩定性評價方法[J]. 巖土力學， 2009，30(3):845-848.

[20]CARPENTER A. CUSVM: a CUDA implementation of support vector classification and regression[EB/OL].(2009-06-17)[2009-11-02].http://patternsonascreen.net/cuSVMDesc.pdf.