基于CUDA的加速MATLAB计算研究
近几年来,计算机显卡核心单元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 矩阵计算
#p#page_title#e#
矩阵测试的数据采用行列数相同的方阵,分别采用维度为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
#p#page_title#e#
从这里可以看出,矩阵计算和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
