GPU加速范文
GPU加速范文（精选8篇）
GPU加速 第1篇
目标检测是计算机视觉领域中一项计算密集的工作。现实生活中的绝大多数应用,如智能监控中的行人检测、汽车安全中的车辆检测和图像检索中的特定类别目标检测等,都要求检测速度足够快以便达到实时性,同时也希望尽可能降低误报率。近年来,研究人员提出了一批优秀的卷积检测模型CDM(Convolutionbased detection model),并在诸如PASCAL VOC[1]、Image Net[2]等极具挑战性的标准评测数据集上取得了出色计算性能。这些模型以卷积运算为基础,用特定的模板与输入图像或特征图进行卷积得到响应图,进而通过最大响应值确定目标位置或将其作为特征图再次进行卷积运算,如图1所示。然而,尽管卷积检测模型简单有效,但卷积运算高强度的计算消耗阻碍了其在现实生活中的应用。
针对卷积运算耗时这一问题,一种有效的方法是利用卷积定理将空间域中的卷积运算通过傅里叶变换到频率域中更为高效的点乘运算,进而显著地降低计算复杂度[3]。尽管如此,在面对实际尺寸目标检测任务时,即使采用多线程并行技术,这种方法在CPU平台上的运行效率依然低下。而随着高性能计算技术的快速发展,GPU的通用科学计算能力得到了越来越多的重视。快速傅里叶变换和点乘运算中蕴含了丰富数据独立性,使得GPU上的计算并行化成为可能,可以进一步提高计算效率,节省处理时间。
本文从算法加速和硬件加速两个角度入手,首先将卷积定理推广至卷积检测模型,通过对比应用卷积定理前后的模型计算复杂度,从数学上论证了算法有效性。接下来,本文分析了算法的并行性,在此基础上提出了基于Open CL的GPU加速策略并进行实现。实验结果表明,本文算法能够有效利用GPU的计算能力,相对于多线程CPU实现,取得了2.13~4.31倍的加速比。
1 相关工作
一直以来,卷积操作在统计、信号处理、计算机视觉等领域都是最基础的数学算子之一,在空间域上通常以滑动窗口的形式实现。在许多早期的研究工作中,卷积算子常被用来调整信号的频率特性,如高斯平滑、平均平滑和Sobel滤波等。此后,模板匹配利用卷积算子在空间域上搜索特定的模式,这被认为是一种简单的检测模型。然而,这类模型难以在包含诸多类别的静态图片中准确地定位特定目标。为解决目标检测中存在的姿态和形态各异的问题,绝大多数常用的检测模型都采用了大量的模板和卷积运算,如形变部件模型DPM(Deformable partbased models)和卷积神经网络CNN(Convolutional neural networks)。尽管这类模型的检测效果显著提升,但其计算量也随卷积算子的增加而线性增长。
1.1 形变部件模型
Felzenszwalb等人[4]提出的形变部件模型采用HOG(Histograms of oriented gradient)描述子[5]作为底层特征,使用树状图模型[6]来表示目标结构,并据此将根模板和若干可形变部件模板组合到一起,形变部件用于反映目标局部形态特性。该模型通过滑动窗口的方式,在HOG特征金字塔上逐层卷积估计目标可能出现的位置。这种计算方式需要处理全部的搜索空间,从而运行效率非常低下。
为了解决此问题,Felzenszwalb等人[7]进一步提出了采用部分假设剪枝策略的级联实现。但由于级联中存在过多控制流和负载不均衡等问题,难以适用于GPU架构。Song等人[8]和Hirabayashi等人[9]进行了形变部件模型的CUDA实现,但同时也将其应用局限于Nvidia GPU平台。De Smedt等人[10,11]对构建HOG特征金字塔进行了Open CL实现,将卷积遗留在CPU端,以流水线在CPUs/GPUs异构体系下实现形变部件模型。相较于特征提取,形变部件模型中的卷积运算计算量更大,也更适于采用GPU并行计算架构。
1.2 卷积神经网络
卷积神经网络[12,13]源于对动物视觉皮层的研究,是多层感知器MLP(Multilayer perceptions)的变形。不同于形变部件模型采用固定的HOG特征,卷积神经网络将模板直接运用于原始图像上,通过大量卷积操作模拟视觉识别系统。卷积层和采样层是卷积神经网络的核心。卷积层主要用于产生不同层次的特征编码,在前馈和反向传播过程中都消耗了绝大部分的计算资源。采样层通过降采样操作,针对局部畸变和简单几何变换等问题实现了一定的不变形。
实际上,形变部件模型也是一种卷积神经网络。Girshick等人[14]简化了由Krizhevsky[15]等人提出的参数规模非常庞大的卷积神经网络,用于构建特征金字塔,以替代HOG特征。同时他们将距离变换推广为采样层,将形变部件模型的推断过程重构为一个卷积神经网络。
受限于网络深度和计算速度,卷积神经网络很难应用于实际尺寸的目标检测任务中。当前主流的卷积神经网络库都是使用CUDA实现GPU加速,如Caffe[16]等。Cecotti等人[17]将傅里叶变换引入到卷积神经网络,但其主要目的是为了便于脑电信号分析而非加速。
2 频率域算法加速
为加速线性目标检测系统的评估速度,Dubout等人[3]引入傅里叶变换,实现了6~8倍的加速比。我们可以将此方法进行推广,用于加速卷积检测模型。
算法加速的核心思想是傅里叶变换的卷积定理,即原始信号在空间域的卷积可由其对应的傅里叶变换乘积的反变换求得,可表示为:
算法加速的基本流程主要包括以下3个步骤:
(1)计算所有空间域特征图和模板的傅里叶变换;
(2)将频率域特征图和模板进行点乘和累加运算,得到频率域响应图;
(3)计算傅里叶逆变换得到空间域响应图。
为深入分析算法加速的有效性,假定卷积检测模型的输入为具有L维特征的J张尺寸为M×N的特征图和K个尺寸为P×Q的模板。同时,将特征图和模板的第l维分别标记为xl∈RM×N和yl∈RP×Q,则由特征图和模板卷积而得的响应图z∈R(M-P+1)×(N-Q+1)可表示为:
其中,表示翻转的模板。
卷积运算的计算复杂度为O(MNPQ)。由于特征图和模板的每个系数之间都需要进行一次乘法和一次加法运算,因此每一维所需的浮点数运算量为:
假定每张特征图都需要与所有K个模板进行卷积,则每张特征图所需的浮点数运算量共为:
其中,MN表示最终为得到1维响应图而进行的累加运算数。为简便起见,假定响应图的大小仍然保持为M×N。
为了在频率域中顺利进行点乘运算,首先需要在傅里叶变换之前将模板填补成与对应特征图同一大小。对于实数傅里叶变换,其变换结果满足Hermitian冗余,即其中一半的变换结果与另一半共轭。因此傅里叶变换后的特征图与模板大小均为M×(N+2),共包含M×(N/2+1)个复数系数。同样,简记其大小为M×N,共包含M×N/2个复数系数。
两维傅里叶变换的计算复杂度为O(M2N2),快速变换算法可将其降低至O(MNlog2MN)。根据文献[18,19],可得到每一维实数特征图傅里叶变换所需的浮点数运算量为:
对于频率域中的点乘和累加运算,特征图和模板对应位置上的复数系数之间需要各进行一次复数乘法(6次浮点数运算)和复数加法(2次浮点数运算),故总运算量为:
在实际系统中,通常只需要载入一次模板即可,因此模板的傅里叶变换可以离线进行。同时,根据傅里叶变换的线性特性,多维响应的累加运算可以在频域进行,故只需要计算一维的逆变换。因此,引入傅里叶变换后,对于每张特征图,共需要进行L次正变换、K次逆变换和KL次点乘和累加运算,所需的浮点数运算量共为:
对比空间域和频率域的浮点数运算量,则理论上的加速比为:
根据式(8),易知若模板数量越多,特征维度越高,即KL>>K+L,则能够获得越高的加速比。同时,对于更大的模板也能够获得更高的加速比。这两点都是当前主流卷积检测模型的共同特性。
3 GPU硬件加速
为了进一步加速卷积检测模型,本文在GPU上基于Open CL对前述算法进行了并行实现,使硬件加速不再受限于GPU平台类型。Open CL是一个面向异构系统,用于通用并行计算的免费开放标准,也是一个统一的编程环境,可广泛适用于各类不同架构的CPU、GPU、DSP及FPGA等并行处理器。
图2所示的算法加速的串行实现伪代码如算法1所示。
算法1频率域算法加速的串行实现
输入:特征图F、模板T
输出:响应图R
在算法1中,存在十分明显的并行性,主要有两个方面:
(1)数据并行性:无论是傅里叶变换、点乘,还是累加,这些基于像素的运算之间都不存在数据依赖性。因此,可以将输入数据映射到Open CL的NDRange索引空间,利用一个线程处理一个像素的运算过程。
(2)任务并行性:算法1中同类型的任务之间不存在数据依赖性,如所有模板的傅里叶变换。而同一FOR循环内的不同类型任务则需要顺序执行,如点乘、累加和响应图傅里叶逆变换。
通过调整FOR循环的顺序,可以使算法1更适宜GPU并行化,如算法2所示。将同一FOR循环中的任务依次发送至GPU,插入到不同的Open CL命令队列,可以使之并行执行。若存在数据依赖性,则将其插入同一命令队列,并通过Open CL事件进行通信,保证执行次序。不同FOR循环间通过同步命令保证内存的一致性。
算法2频率域算法加速的并行实现
输入:特征图F、模板T
输出:响应图R
3.1 快速傅里叶变换
N点一维离散傅里叶变换(DFT)可表示为:
其中,0≤k<n,ωn=exp(-2πi/n)。在诸多快速傅里叶(FFT)算法中,Cooley-Tukey算法[20]最为常用。其以分治法为策略,递归地将长度为n=n1n2的DFT分解成为两个长度分别为n1和n2的较小规模DFT以及旋转因子的复数乘法。式(9)可重写成:
其中,j=j1n2+j2,k=k1+k2n1,称作旋转因子,n1或n2称为基。计算示例如图3所示。
由于输出并非顺序,Cooley-Tukey算法中存在显式位反转操作,这并不适合在GPU上实现。Stockhan算法[21]通过在每次分解中穿插进行多维转置,可以避免这一问题而直接得到顺序输出,适于GPU实现[22]。
在利用GPU并行计算2维FFT时,通常的方法是使用1维FFT分别变换其所有行和列。在分解完成后,FFT核函数将数据从显存加载到寄存器中,利用GPU的局部内存存储中间数据,每个线程递归地计算一个基的FFT。
3.2 拼接策略和优化
前文曾提及需要在傅里叶变换之前将模板填补成与对应特征图同一大小。而简单的填补方法或将导致过高的内存开销,或需要消耗较多的计算资源,特别是对于部分需要构建特征金字塔的卷积检测模型。
对于上述问题,可以通过利用一种快速启发式左下装箱算法[23]拼接特征图来解决。其基本思路是将所有特征图组合成为若干特征拼图,后将模板填补成同一大小[3],如图4所示。特征拼图的大小通常等同于最大的特征图。
为避免GPU产生不必要的调度开销,不仅需要控制GPU上任务的数量,也需要合理设置每个Open CL核函数运行时所启动的线程数。因此,本文将特征图和模板的点乘运算和累加运算合并到一个核函数中,使用3维NDRange索引空间,每一个工作组计算一个像素上L维特征的相关运算。这样不仅可以避免过多的调度开销,还能够更为规律地进行内存合并访问。同时有计划地发送任务,适当的同步也可降低调度开销。
4 实验结果与分析
本文选择形变部件模型作为基础,进行Open CL实现,GPU上的FFT计算使用cl Math数学库中的cl FFT函数库。基准CPU多线程实现为Idiap研究所的FFLD[3],其使用CPU SIMD指令集和Open MP实现多线程并行计算,FFT计算使用FFTW函数库,矩阵运算使用Eigen函数库。
本文在两个不同硬件平台上进行了实验,一个配置为Intel Core i7 3770(3.4 GHz,4核)和Intel HD Graphics 4000(0.35 GHz,128个流处理器);另一个配置为Intel Core i5 3570K(超频4.2 GHz,4核)和Nvidia Ge Force GTX 670(超频1.17GHz,1344个流处理器)。
测试数据来自PASCAL VOC 2007数据集,检测模板来自基准CPU实现。以自行车为例,其检测模板共54个,模型构建的32维HOG特征金字塔共17层,组成6个特征拼图。根据算法2可将整个卷积过程(从输入特征图和模板到输出响应图)分为四部分:模板正变换、特征图正变换、点乘和累加以及响应图逆变换。表1给出了各部分总的计算时间,包括GPU计算时间、数据传输时间,以及各类在CPU上进行的预处理或后处理的计算时间,以便更客观地对比本文算法的加速效果。
在Intel GPU平台上,点乘和累加运算上取得了4.26倍的加速比,而其他三部分则更慢一些,模板正变换、特征图正变换和响应图逆变换的计算时间分别为基准的0.58、0.83和0.23倍。同样地,在Nvidia GPU平台上,本文实现了25.16倍的点乘和累加运算加速,而其他三部分则分别为基准的0.91、0.94和0.25倍。
实验结果表明,点乘和累加运算能够充分利用GPU中的大量流处理器,从而实现较高的加速比。而有关FFT运算的其他三部分较基准更慢的一个主要原因是其中数据传输以及CPU上的各类预处理或后处理操作所需的计算时间较高。例如,响应图逆变换速度远慢于基准的原因便是分解响应拼图操作中数据传输耗时较多。另一个主要原因在于所选用的FFT函数库。由于目前便于使用的以Open CL实现的FFT函数库较少,本文选用的cl FFT函数库仍处于开发初期,其功能和效率还有待提升。事实上,基准CPU实现所使用的FFTW函数库能够利用批处理操作,同时精确计算32维HOG特征图和模板的快速傅里叶变换。而cl FFT函数库在HD 4000上却只能同时精确计算5维,在GTX670上也不过增加至11维。对此,在具体实现中不得不先将特征图和模板进行拆分,计算FFT后再将其融合。
另外,GTX 670上的各部分计算时间均少于HD 4000,加速比分别为1.83、1.41、5.45和1.25倍。这主要得益于GTX 670拥有更高的工作频率和更多的流处理器,能够产生更多的并发线程。
由于模板的傅里叶变换可以离线进行,所以若从实际应用的角度出发,可将特征图正变换、点乘和累加,以及响应图逆变换这三部分计算时间的总和作为标准,来衡量卷积在线流程的性能。在HD 4000上,卷积在线流程的计算时间由2097 ms降至984 ms,提速2.13倍,其中点乘和累加运算的计算时间占比从80.88%降至40.45%,如图5所示,其中饼图的面积表示卷积在线流程所需要的计算时间。在GTX 670上,卷积在线流程的计算时间由2162 ms降至502 ms,提速4.31倍,而点乘和累加运算的占比更是明显下降,由84.97%降至14.54%,如图6所示。
综上所述,尽管在Open CL实现中所选用的第三方函数库的效率并不理想,从而导致部分运算时间较CPU实现有所增加。但从整体上来看,对于卷积检测模型本文的GPU实现在不同硬件平台上仍然能够取得不错的加速效果。
5 结语
本文提出了将傅里叶变换推广至卷积检测模型,通过将空间域中的卷积运算转换到频率域中的点乘运算以大幅降低计算复杂度,并从数学上论证了其有效性。同时,本文在GPU上使用Open CL实现了进一步的硬件加速。实验结果表明,尽管本文所选用的第三方函数库效率并不理想,但本文的GPU实现仍然能够在不同架构的硬件平台上取得不错的加速效果。相对于4核CPU平台下的多线程实现,本文的GPU实现能够达到2.13~4.31倍的加速比。
GPU加速 第2篇
说到高清大家印象里一定是蓝光、1080p等字眼，体积大、购买难让高清成为了镜中花水中月。我们手中的大显示器、液晶电视往往观看的都是DVD画质的标清视频。在等待蓝光普及的日子里，很多人为提前享受高清也想出了办法，其中倍线就是方法之一。
倍线，一个简单而又神秘的技术
倍线的原理实际上说起来很简单，大家都知道广告商宣传的399购买“千万级像素”数码相机是采用插值，用软件补偿的方法强行拉大图片。那么倍线也就很好解释了，它也一样是用软件把原本只有标清画质的DVD拉伸至高清的分辨率，分析前后数帧的画面，将其中的有效色彩信息提取出来补全到每一帧当中，实现提升色彩表现、清晰度和对比度的目的。
图1
●使用倍线后的效果很明显，无论是细节还是色彩都有一定的改善。
如何实现倍线
常见于数码影像产品附赠软件光盘中的ArcSoft公司的TotalMedia Theatre软件，便能够利用GPU运算能力，将低分辨率的DVD视频差值处理为高清视频实时播放出来。
第一步：下载TotalMedia Theatre软件（下载地址：http://work.newhua.com/cfan/200913/TotalMedia.rar），并安装。
第二步：打开TotalMedia Theatre，在播放窗口点击右键，选择设置或者按Ctrl+S键（见图2）。
图2
第三步：在打开的设置窗口中，选择左侧的视频，在取消启用硬件加速，并启用ArcSoft SimHD（见图3）。
图3
由于硬件解码会与倍线产生冲突，所以我们在使用倍线的时候要需要取消硬件加速。
小提示:看什么格式的视频可以用倍线
GPU加速 第3篇
近年来,随着计算机技术的发展,图形处理器(GPU)越来越强大,并且在计算能力上已经超过了通用CPU。使用GPU计算可以以低廉的价格获得巨大的计算性能,因此成为了科学计算领域的一个应用热点。自2007年NVIDIA公司推出了CUDA运算平台,并使用C语言为CUDA构架编写程序以来,GPU计算技术已经广泛应用于诸如信号处理、图像处理、信息安全等热门领域中。在数值线性代数中,GPU可以用于加速大规模的矩阵计算问题,包括矩阵的分解和求特征值以及特征向量等。许多组织和机构在GPU上实现了LAPACK库中的函数,并发布了相关软件包以供科研人员使用。比较著名的软件包有CULA和MAGMA等。
Cholesky分解是矩阵计算中的一个基本分解问题,常常用于求解系数矩阵为对称正定矩阵的线性方程组,由于其计算量比使用LU分解求解一般线性方程组的算法约减少一半。因此在科学计算中也有广泛的应用。此外,Cholesky分解也应用于Kalman滤波
CUBLAS是在GPU上实现的BLAS。利用CUBLAS库可以很方便地移植CPU代码到GPU上进行加速计算。自2008年Volkov等人提出CPU-GPU混合计算矩阵的LU、Cholesky和QR分解算法
1 使用CPU计算Cholesky分解
给定一个对称正定矩阵A,则存在主对角元素全为正数的下三角矩阵L满足:
通常称下三角矩阵L为矩阵A的Cholesky因子。下面给出Cholesky分解的直接算法
1.1 直接Cholesky分解
对于对称正定矩阵A直接由式(1)得到如下等式:
则对于第j个对角元素ajj,有:
对于位置为(i,j)(i>j)的元素,我们可得如下等式:
由此可得到求解Cholesky分解的公式如下:
若令j=1,2,…,n,循环使用式(5)、式(6)计算L的下三角元素,这样便得到计算Cholesky分解的直接算法。由于该直接算法主要是Level 2 BLAS操作,不适合移植到GPU上计算。
1.2 基于块的Cholesky分解
BLAS库内部实现了各种线性代数的基本运算。由于其代码的高效性和完整性,因此一直为科研人员广泛使用。其中Level 3 BLAS是依据现代计算机的内存等级进行了大量的优化处理,使得其计算性能远高于Level 2 BLAS。类似的结论也在GPU上成立。
为了能利用矩阵的子块进行计算,下面使用Cholesky分解的块递推更新算法。对矩阵A进行分块,如下:
则由式(1)有:
于是得到如下递推公式:
对矩阵按式(7)-式(9)依次递推更新下去,便可得到Cholesky分解。
具体的递推过程如图1所示。
图1 基于块的Cholesky分解
1.3 算法及分析
基于块的Cholesky分解算法如下:
Step1使用Intel MKL库中的dpotrf函数计算第一个子块A1的Cholesky分解;
Step2使用BLAS库中的Level 3函数trsm求解矩阵方程X×A1T=A2;
Step3使用BLAS库中的Level 3函数syrk更新尾部矩阵A3=A3-A2×A2T;
Step4对尾部矩阵A3重新分块并针对A3重复执行以上三步,直到计算完最后一个子块的Cholesky分解则结束。
在余部矩阵A3的阶数较大时,Step 3占据一次循环约80%的计算时间。而A1由于分块大小的固定,计算量不会改变,Step 1所占的计算时间非常小。当A2的行远大于分块值时,Step 2几乎占据了剩余20%的计算时间。
2 使用GPU计算Cholesky分解
2.1 CUBLAS库
CUBLAS是由NVIDIA公司提供的并在CUDA构架的GPU上实现BLAS的函数库。利用CUBLAS库可以很方便地移植矩阵计算的CPU代码到GPU上进行计算。由于矩阵算法中主要涉及矩阵与矩阵相乘等的Level 3操作以及矩阵和向量相乘等的Level 2操作,且BLAS和CUBLAS库中的函数是经过高度优化的,并被全世界广泛认可,因此在计算Cholesky分解中,可直接使用该库中的函数在CPU或GPU上实现矩阵的Level 2或Level 3操作。这样,可以将主要精力放在Cholesky分解算法的任务划分及任务调度上。
2.2 使用CUBLAS计算Cholesky分解
由于CUBLAS库中提供了trsm和syrk函数,但没有提供potrf函数。对此,可以在每次循环时先将需要计算Cholesky分解的子矩阵块传回到CPU上并利用MKL中的potrf函数计算其Cholesky分解,然后在将结果传输至GPU对应的位置。再依次使用CUBLAS库函数在GPU上更新矩阵的剩余部分。然而这样做的结果是当GPU(或CPU)执行计算任务时,CPU(或GPU)处于闲置状态,这导致了硬件资源的浪费。使用CUBLAS计算Cholesky分解的具体过程如图2所示。
图2 使用CUBLAS加速Cholesky分解
这样直接使用CUBLAS库函数加速实对称正定矩阵的Cholesky分解的效果也是十分明显的。首先由于在GPU上的操作全为Level 3的BLAS操作;其次,相对较小的分块值nb,如nb<512,在本实验环境下,使用MKL库函数计算512阶对称正定矩阵的Cholesky分解耗时约为4 ms,而来回传输512阶矩阵总共耗时约2 ms。相比较大的对称正定矩阵而言,直接使用CUBLAS库函数已经可以称之为充分利用GPU进行计算了。
2.3 Volkov的混合算法及分析
Volkov等人于2008年给出了CPU-GPU混合计算Cholesky分解的算法
图3 第1步更新
图4 第2步更新
图5 第3步更新
图6 第4步更新
Volkov的混合算法巧妙地移动了矩阵循环的更新顺序,这样余部矩阵的上次更新和本次子块的Cholesky分解及其所在列的更新可以同时进行。
考察Volkov算法中的混合并行部分,对GPU上的某次迭代更新后,有如图7所示的剩余矩阵:
其中矩阵A1、A2、A3为剩余待更新矩阵,子块B1、B2为已求出的Chol-esky分解部分。之后,在GPU上按如下公式对A1与A2进行更新:
然后将A1与A2异步传输到CPU。同时,在GPU上对余部矩阵按如下公式进行更新:
在CPU上再次对A1和A2更新:
其混合并行的示意如图8所示。
图8 原Volkov的混合调度策略
可以发现在GPU上更新矩阵A1、A2与更新A3是相互独立的。而且当剩余待更新矩阵阶数较大时,可以考虑将A1的首次更新移到CPU上进行。注意,对A2的两次更新不能同时移到CPU上进行,因为对A2的两次更新计算量都比较大,在CPU上计算时,耗时很多,极容易超过GPU上对A3的更新耗时,从而导致GPU的空闲。
此外,当剩余待更新的矩阵较小时,在GPU上更新A3的时间较短。此时不适合再将子矩阵A2转移到CPU上进行更新,否则会使得GPU出现长时间空闲。
整个混合计算过程中对数据传输采用异步方式,使得每次数据传输时间被GPU计算或CPU计算隐藏。
对此,我们综合考虑了计算的CPU以及GPU的计算性能,以Volkov的混合算法为基础,给出了一种新的调度策略。
2.4 新的混合调度算法
记剩余待更新的矩阵A3的阶数为ns。将新的混合调度算法分为以下三个阶段:
阶段1当ns>x时,更新子矩阵A3耗时较长,因此CPU可以做更多的工作以减少空闲时间。由于对子矩阵块A2的两次更新只能有一次移到CPU上进行,故有以下2种调度策略:
1)将A2的第一次更新放在GPU上执行,如图9所示。
图9 阶段1的混合调度策略1
2)将A2的第二次更新放在GPU上执行,如图10所示。
图1 0 阶段1的混合调度策略2
通过分析A2两次更新的计算量,可以发现A2的前次更新计算量为A2再次更新计算量的两倍。为防止GPU出现空闲,只有当ns>x1(x1>x)时,GPU上对A3的更新任务能完全隐藏CPU上的计算任务,采用混合调度策略2;否则,x<ns<x1时,将采用混合调度策略1,以防止GPU的空闲。这样相互选择调度是最优的结果,其x1、x2也需要根据计算机的CPU及GPU的计算性能大致估计出来。
阶段2当y<ns<x时,遵循Volkov的混合并行策略。混合调度策略参见图8所示。
阶段3当ns<y时,在GPU上首先更新A1,然后启动异步传输A1到CPU上进行再次更新操作。其余的更新操作在GPU上完成,如图11所示。
图1 1 阶段3的混合调度策略
其中x与y的值需要依据各自计算的CPU和GPU的性能来大致确定,以保证GPU不会出现空闲时间并且CPU的空闲时间最小。这样的混合调度策略则是非常合适的。
同时,要对Volkov混合算法的第一步和最后一步进行调整。在第一步中,对第一主列块的更新与余部矩阵传输到GPU上可以异步进行,这样其中的一个耗时会被隐藏。最后一步更新时,可以先将最后一个子块先传回至CPU,然后将更新最后一个子块与传回剩余矩阵数据异步执行。
对于分块大小nb,可以参考MAGMA软件包中关于分块值的设定。针对Kerpler构架的GPU,可将分块值设定如下:
3 数值实验及结果分析
3.1 实验环境
本文中描述的GPU计算过程均是在NVIDIA的Kerpler构架的GPU上实现的。有关数值试验的硬件和软件平台见表1所示。
表1 实验环境
3.2 实验结果与分析
实验测试1000~10 000阶对称正定矩阵,而且使用双精度数据类型,以保证计算结果的精确性。
首先分别给出本实验环境下MKL中执行dgemm以及dtrsm的耗时以及CUBLAS中一执行dgemm,dtrsm以及syrk的耗时的一个参考值,如表2所示。该数据用于估计划分阶段的x(包括x1)与y值。
表2 各函数测试结果
(单位:毫秒)
其中针对dgemm函数,m为表格的第一列,n=k=512。而dtrsm与dsyrk函数的n为表格第一列,k=512。
表3给出了数据矩阵在CPU与GPU之间的传输耗时的一个参考值。
表3 数据传输测试(单位:毫秒)
下面考虑x=6000时,使用阶段1的混合调度策略1是合适的。首先依据表3可得到,异步传输A1,B1(矩阵规模均为512×512)以及在CPU上更新A1的时间被GPU上更新A2隐藏。而异步传输A2(矩阵规模6000×512)耗时t1<20 ms,在CPU上再次更新A2耗时t2<80 ms,异步传回A1,A2耗时t3<20 ms,而GPU上对矩阵A2(规模6000阶)的更新耗时为t≈150 ms>(t1+t2+t3),故可取x=6000。类似,我们可以由此得到比较合理的x1,x以及y的值,分别如下:
最后,对1000~10000阶对称正定矩阵分别使用4种方法进行计算。表4给出了4种方法计算1000~10 000阶对称正定矩阵的Cholesky分解的结果。其中可以看到基于新调度策略的函数New_Sche明显优于使用CUBLAS直接计算Cholesky分解以及基于Volkov的混合算法的函数Volkov_S。
表4 各算法测试结果
(单位:秒)
4 结语
科学计算中使用CPU与GPU混合并行计算的例子层出不穷,然而最困难的在于设计一种负载均衡的混合调度方案,从而使得CPU的空闲时间尽可能小,同时又不能让GPU空闲。糟糕的混合调度方案会使得整体计算性能下降,或者使得计算不正确,从而失去了混合计算的目的,这都是需要避免的。合适的混合调度策略往往需要依据自身计算机的CPU以及GPU的计算性能来共同确定。首先细致分析每步操作在CPU以及GPU上的计算速度,然后合理划分计算任务到CPU和GPU上执行。再者,算法中会大量使用异步传输操作,使得这些数据传输耗时被CPU和GPU上的计算时间隐藏,这也是优化加速的一个重要部分。最后,本文中没有考虑分块值对混合调度策略的影响。由实验的分析可知,分块值与调度策略有直接关系,影响阶段值x1,x以及y的确定。如何确定最优的分块值以及相应的调度策略将是后面要进行的工作。
摘要：针对大型实对称正定矩阵的Cholesky分解问题,给出其在图形处理器(GPU)上的具体实现。详细分析了Volkov计算Cholesky分解的混合并行算法,并在此基础上依据自身计算机的CPU以及GPU的计算性能,给出一种更为合理的三阶段混合调度方案,进一步减少CPU的空闲时间以及避免GPU空闲情况的出现。数值实验表明,当矩阵阶数超过7000时,新的混合调度算法相比标准的MKL算法获得了超过5倍的加速比,同时对比原Volkov混合算法获得了显著的性能提升。
关键词：图形处理器,乔里斯基分解,加速比,混合算法
参考文献
[1]Chandrasekar J,Kim I S,Bernstein D S,et al.Reduced-Rank Unscented Kalman filtering using Cholesky-based decomposition[C]//American Control Conference,June,2008:1274-1279.
[2]Yu H,Chung C Y,Wong K P,et al.Probabilistic Load Flow Evaluation With Hybrid Latin Hypercube Sampling and Cholesky Decomposition[J].IEEE Transactions on Power Systems,2009,24(2):661-667.
[3]David S.Watkins.Fundamentals of Matrix Computations[M].New York:John Wiley and Sons,2013.
[4]Gene H Golub,Charles F,Van Loan.Matrix Computations[M].Baltimore:Johns Hopkins University Press,2013.
[5]Volkov V,Demmel J W.Benchmarking gpus to tune dense linear algebra[C]//Proceedings of the 2008 ACM/IEEE Conference on Supercomputing,Nov,2008:1-11.
[6]胡鹏飞,袁志勇,廖祥云,等.基于CPU-GPU混合加速的SPH流体仿真方法[J].计算机工程与科学,2014,36(7):1231-1237.
[7]张健,焦良葆,陈瑞.CPU-GPU混合平台上动态场景光线跟踪的研究[J].计算机工程与应用,2012,48(21):151-154.
[8]Yaohung M Tsai,Weichung Wang,RayB ing Chen.Tunning Block Size for QR Factorization on CPU-GPU Hybrid Systems[C]//Proceedings of the IEEE 6th International Symposium on Embedded Multicore Socs,Sept,2012:205-211.
[9]John Cheng,Max Grossman,Ty McK ercher.CUDA C Programming[M].Indianalpois:John Wiley&Sons,2014.
GPU加速 第4篇
近年来,图形处理器(Graphic Processing Unit,GPU)的计算能力远远超过其它通用平台,可编程性越来越灵活,将体绘制算法采用到GPU可编程单元进行加速处理,可大大提高系统的实时性,并具有很高的性价比。本文以传统的光线投射算法为基础,充分挖掘GPU的并行处理能力,实现了GPU加速的医学影像光线投射体绘制,在提高绘制图像质量的同时又保证了交互的实时性。
1 GPU加速的体绘制算法
基于GPU的体绘制技术主要可分为基于纹理映射的体绘制和光线投射的体绘制两种。
1.1 GPU加速的纹理映射体绘制
基于纹理映射的体绘制由Cullip和Neumann首次提出[3],经Cabral等进行了进一步研究[4]。这种方法将某种采样平面作为纹理,根据体数据空间坐标和纹理空间坐标的转换关系,充分挖掘GPU的纹理映射能力进行绘制。采样平面可以是基于视线空间的三维纹理[5],也可以是基于物体空间的二维纹理序列[6],按照所谓的代理几何体(proxy geometry)对纹理数据进行重采样,由从前到后或者从后到前的顺序对各个片层进行绘制。
这两种基于纹理映射的方法广泛应用于中等大小的数据集的体绘制,经过很多改进[7,8],在绘制图像质量要求不高的情况下可达到交互的帧率。但是,这两种方法由于其采样的片层是一个平面,所对应的光线是平行的,因此只适用于正交投影,无法进行透视投影,限制了其应用。
1.2 基于GPU的光线投射体绘制
与基于GPU纹理映射的体绘制算法相比,基于GPU的光线投射体绘制算法既可以进行正交投影,又可以进行透视投影,而且实现方式更为灵活,为空间跳跃等改进算法提供更多的方便,是一种能够兼顾速度和绘制质量的方法[9]。
国内外利用GPU对光线投射算法进行加速,有很多的研究。其中Simon Stegmaier等人在GeForce 6800GT平台上实现了加速算法,程序可实时加载不同类型的着色器,达到不同的绘制效果。但是,由于采用NVIDIA公司特定的OpenGL扩展,缺少通用性,无法在ATI显卡上运行[10]。Ruijters等人主要研究了GPU加速的优化算法,如分块、KD树和光线提前结束等方法[11]。Henning等人主要研究基于GPU的内窥镜重建,实现了基于GPU的透视投影重建[9],但是由于当时并不支持OpenGL的高级着色语言GLSL,着色器仍然采用低级的着色语言。
2 GPU加速的光线投射体绘制算法实现
在支持着色模式3.0的GPU出现以前,由于不支持动态流程分支,基于GPU的光线投射算法需要分成几个渲染通道来实现,例如利用深度测试和碰撞检测来模拟循环的执行[12]。而现在的GPU大部分支持动态流程分支和循环,这使光线投射算法的实现更为简单。本文提出的光线投射算法即采用单个渲染通道,着色语言选择GLSL,具有更强的通用性。而且,这种算法能够利用GPU如下的特性得到速度提升:
(1)硬件插值自动得到光线位置;
(2)内置的三维纹理快速三线性插值;
(3)几乎无需代价的全浮点合成;
(4)充分利用反射、点积等内置函数,计算光照效果;
(5)利用动态流程分支,实现光线提前结束。
2.1 像素算法原理
本文提出的GPU加速的光线投射算法,由两部分构成。一部分是由着色语言GLSL编写的顶点着色器和片元着色器[13],利用图形硬件的可编程性,使传统光线投射的6个步骤(光线生成、光线遍历、重采样、分类、光照和颜色合成),得以完全通过GPU实现;另一部分由C++和OpenGL编写,完成数据的准备(如加载纹理、计算变换矩阵和设置光源、视点等)实时编译和连接着色器等功能。GPU加速的光线投射算法主要步骤如下:
(1)将体数据作为三维纹理载入内存。如果选择实时计算梯度,则直接将体数据值载入三维纹理的alpha通道;如果选择预处理时计算梯度,则将计算好的每个体素的梯度值载入3维纹理的RGB颜色值,将体素的灰度值载入alpha通道。另外,分类得到的转换函数查找表也作为纹理载入显存。
(2)顶点着色器进行顶点变换,插值得到在裁剪空间中的片元位置。
(3)片元着色器部分是整个算法的核心。主要实现光线对体数据的遍历、重采样和计算。先计算光线的开始点和离开点,求得方向和循环次数,后从开始点进行遍历循环。利用硬件的自动3线性插值功能进行重采样,索引转换函数查找表得到每次步进后当前体素点的颜色值,计算光照效果,反复合成得到像素最终颜色和不透明度。
片元着色器的伪代码如下所示:
2.2 基于片元着色器的高质量光照
在传统的基于三维纹理映射的方法中,表示法向量的体数据场需要在数据预处理中预先计算,并做归一化处理,在绘制时实时加载入图形硬件中。当对体数据进行采样时,在对体数据插值的同时,也对预处理得到的法向量进行了三线性插值。插值后的梯度一般不是归一化的,这些非归一化的法向量用于计算光照,严重影响了光照效果。另一方面,传统的固定管线采用逐顶点的Phong光照模型,先在每个顶点上计算光照,然后插值得到每个像素上的光照量。显然,这种方法会导致图像质量的下降。
本文算法克服了以上的缺点,能够计算得到归一化的法向量,得到高质量的绘制效果。即采用逐像素的Phong光照模型,利用可编程片元处理器对每个像素进行光照计算。Phong光照模型的公式如下:
其中L为入射光向量,N为法向量,R为反射光向量,s是一个和物体有关的镜面反射指数,反映了物体表面的光滑程度,Ka、Kd、Ks分别是环境光系数、漫反射光系数和镜面反射系数。编程实现时可利用着色语言的内置函数求点积、反射光以及进行归一化操作,充分发挥图形卡的性能。
其中,光照计算所需的法向量N由规格化的梯度矢量近似得到。梯度在三维空间中是一个指向某个方向的向量,是一个定位信息,因而可以给出数据集合中结构的方位信息。本文的梯度采用一个三维中心差分梯度估算器得到:
梯度可以预先计算,作为三维纹理进行加载,以提高绘制速度,但是所需要的纹理内存是体数据的3倍,不利于处理大规模体数据。梯度也可以在片元着色器中,通过访问六邻域进行实时计算,以减少预处理时间和纹理内存的消耗,但会降低绘制速度。这两种方法均在后面所设计的工具包中编程实现,可供选择。
3 GPU加速的体绘制工具包设计
VTK(Visualization Toolkit)类库利用面向对象的技术,将算法封装成一系列定义清晰的、易于扩展的类,便于增加新的功能,并且源码开放[14],本文通过扩展VTK类库,设计了一个能够充分利用GPU计算能力的体绘制工具包,使系统开发人员无须了解硬件算法,也可使用体绘制的硬件加速功能,从而提高系统实时性和交互能力。
通过对VTK类库的研究知道,与光线投射算法相关的VTK类为vtkVolumeRayCastMapper和vtkVolumeRay CastFunction,其中vtkVolumeRayCast Function是纯软件的光线投射算法,即计算一条光线上的积分,得到最终图像上该像素点的颜色值;vtkVolumeRayCastMapper负责准备光线投射算法所需的参数、调用光线投射函数和显示。因此,我们可以对这两个类进行继承以得到新的GPU加速的光线投射类vtkGPUVolumeRayCasterMapper和vtkGPUVolumeRayCastFunction,即可实现GPU加速的体绘制再加速。新类vtkGPUVolumeRayr继承了标准类的参数准备和对光线投射函数的调用,增加了纹理加载和图像空间四边形绘制功能。而新的vtkGPUVolumeRayCastFunction调用了OpenGL API函数,加载了光线投射核心算法的着色器。
本文所设计的工具包特点是易于使用和扩展。开发人员使用GPU加速的体绘制工具包和使用标准的VTK软件工具包一样简单,并且工具包易于扩展,方便添加新的功能。工具包尽可能地保留了VTK的架构和设计,用户可以无需OpenGL和GPU加速的知识而直接调用工具包,实现硬件加速的功能。在数据处理方面,工具包支持8位和16位的数据,支持梯度预先计算和实时计算;在光照方面,运行时可以实时改变光照状态和光源位置。
4 实验结果和讨论
本文设计的工具包测试环境如下:CPU为Pentium 42.66GHZ,内存为1G DDR2,GPU为NVIDIA GeForce8800GT,显存为512MB,图形连接器为PCI EXPRESS16X。编程环境是Windows XP SP3下的Visual Studio6.0,图形API选择OpenGL 2.0,着色语言选择OpenGL的着色语言GLSL。
4.1 绘制效果
实验对VTK工具包中自带的切片数据集(64*64*93)进行体绘制,结果如图1所示,图1(a)是本文提出的GPU加速的光线投射体绘制结果图,图1(b)是VTK软件实现的光线投射体绘制结果图,图1(c)是VTK标准库中基于二维纹理映射体绘制的结果图,图1(d)是VTK标准库中基于三维纹理映射体绘制的结果图,几种算法采用同样的传递函数。可以看出,基于二维纹理映射的算法由于缺少层间插值,成像质量最低。其它三种算法,静态图像质量都比较高,其中图1(a)GPU加速的光线投射算法由于使用全浮点的硬件三维纹理插值和合成计算,以及逐像素的光照计算,因此光照效果最为满意。而图1(d)基于三维纹理映射的体绘制,放大和旋转时会产生很明显的伪影。通过对比得出,采用本文的GPU加速的光线投射算法成像质量最高。
4.2 绘制速度
对不同大小的医学切片体数据集采用不同的绘制算法,绘制窗口大小均为800*600,测得的绘制帧率如表1所示。由表1可以看出,GPU加速的光线投射算法绘制速度最快,比软件算法速度提升了300多倍,而且绘制速度受体数据大小影响最小。表2给出了实时计算梯度和预先计算梯度的速度比较。实时进行梯度计算,可加快数据载入,减少预处理时间,但是因为纹理访问带来的存储器延迟[15],造成了绘制帧率的下降,这种情况对于比较大的数据集的绘制尤为明显,而对小型数据集影响则并不明显。本文设计的工具包可以根据用户需要,在实时计算梯度或者预先计算梯度间进行计算。
5 结论
本文充分挖掘了GPU强大的计算能力,实现了基于图形硬件的光线投射算法,对医学影像三维体绘制过程进行了加速,并采用了一种逐像素光照的绘制方法,能在60帧/s的帧率下得到高质量的图像绘制效果,能满足大规模数据集实时显示的要求,与基于纹理映射的方法相比,具有一定优势。同时设计了一个简单易用的医学影像三维体绘制工具包,使GPU加速的光线投射算法具有一定的实用性和推广性。
摘要：为了提高体绘制的速度,充分利用了GPU的可编程性和并行计算能力,在兼容着色模式3.0的图形卡上实现了光线投射算法,并通过对可视化工具包VTK进行扩展,设计了一个可充分利用GPU加速功能的体绘制工具包。实验结果表明,GPU加速的算法显著提高了绘制速度,能在60帧/s的速度下得到高质量的图像绘制效果,并保证了交互的实时性。所设计的工具包具有一定的可用性和可扩展性。
GPU加速 第5篇
可满足性求解(satisfiability)是指,基于给定的一组布尔变量V={v1,v2,,vn},及由这些变量构成的短句集合C={c1,c2,,cm},求得满足该短句集合C的一组解V={v1,v2,,vn}的算法,简称为SAT[1]。SAT目前被广泛用于多个应用领域,如软件验证、理论证明、微处理器验证、模块验证等,其最重要最成功的应用之一为数字系统的测试与验证[1]。
SAT求解的工业实例问题规模很大,求解变量数目已达到百万数量级,约束子句数目达到千万数量级,巨大的求解时间开销难以满足需求。因此,设计高效的求解算法对于实际应用具有重要的意义。SAT用于工业实例求解时,需要对大量的变量进行迭代传播,每次传播需处理的数据量很大,且不同数据传播过程之间无需信息交互,为并行化提供了前提。
GPU具有大量的处理单元,可对大规模数据进行并行处理。GPU计算性能的提高和编程模型的发展,使得GPU越来越多地应用于通用计算中,如线性代数[2]、分子动力学[3]、统计学[4]、工程科学[5]等。本文针对于SAT求解特点及GPU特性,提出一种CPU与GPU协同求解的思想,基于GPU设计并实现了SAT算法中的最耗时(占算法开销80%～90%)的BCP过程的并行化。最后通过随机生成k-SAT测试用例进行测试,取得了5.4～10.3加速比。
1 相关工作
1.1 SAT
SAT分为完全SAT和不完全SAT两类。完全SAT对于整个求解空间进行搜索,如有解,则给出;如无解,可证明问题不可满足。不完全SAT由于求解中带有随机性,不能用于证明。本文基于完全SAT进行研究。现有的完全SAT算法中主要应用包括一下几流程:分支选择(branching),BCP,冲突学习(Conflict-driven clause learning)与回退(Backtracking)[1]。其中冲突学习与回退的出现对于SAT的发展具有里程碑的意义,并出现了一些经典的求解器,如GRASP[6],CHAFF[7],SATO[8]等。现有主流完全SAT求解器主要基于DPLL(Davis-Putnam-Logemann-Loveland)算法框架(如图1所示)来实现的[9]。
make_branch_decision() 通过启发式策略从现有的自由变量中选取一个作为决策变量,对其进行赋值指派,该过程也叫做决策(decision)过程;
deduce() 函数对新指派的变量以单文字蕴含(unit implication)规则基于约束集进行传播;运用单文字规则传播的全过程称为布尔约束传播BCP。如果一个变量被不同的子句通过单文字规则同时指派成TURE和FALSE,则说明当前的指派存在冲突,算法进入冲突分析(conflict analysis)过程;
conflict_and_backtrack() 分析当前冲突产生的原因,并表示成子句形式插入到现有的约束集中,这个过程叫做冲突学习(conflict-driven learning);通过冲突学习会撤销一部分变量的指派,这个过程称为回退(backtracking)。
1.2 并行SAT
伴随着多核体系结构和多机系统的出现,相继出现了一些并行SAT算法。如PMinisat[10]中实现了多线程之间对短子句的共享;PSatZ[11]设计了Master/Slaver模型,从线程分别与主线程通信,由主线程动态调整已实现动态的负载均衡;PSATO[12]中也采用了Master/Slaver的结构,并提出了将从进程的信息进行备份的容错措施。PaSAT[13]中设计了共享空间来实现对于学习到的子句的共享。这些并行的SAT算法都是在原有算法的基础上,通过对搜索空间的划分,将问题划分成若干个子问题由多个线程分别执行求解,在每个线程内部仍然是原有的串行SAT算法。
近年来出现了很多硬件加速器,包括图形处理器(GPU),元件可编程逻辑阵列(FPGAs)以及Cell B.E.,利用这些硬件加速器的很多研究都取得了比较好的效果[14,15]。GPU最初只用于图形处理,伴随其性能提高和编程模型日益成熟,如NVIDIA的CUDA,AMD的Brook+等,如今已广泛应用于通用计算领域[2,3,4,5]。通过提取算法的可并行部分,设计计算核心,将任务以合理粒度划分成多个线程,经过线程处理器调度分配到各个流处理单元执行,即可实现程序的并行化处理。
John D.Davis等在文献[16]中提出通过FPGA实现对SAT的加速,并融入了冲突学习过程,取得了很好的实验效果。基于GPU的SAT并行研究很少,主要因为受GPU每个处理单元功能的局限,其对于复杂问题的应用还不成熟。Luo Zhongwen等在文献[17]提出将3-SAT问题编码并存储到GPU中来实现。该文的主要工作是将3-SAT算法在GPU上实现,验证了GPU对3-SAT的适用性,但是并没有充分利用GPU的特点进行算法加速的研究,因此算法的性能提高并不理想。
2 基于GPU的SAT算法BCP过程加速方法
为提高性能而出现的许多启发式策略使得部分算法结构粘性很大,因此完全将现有的“类chaff”的SAT移植到GPU上实现难度很大,且对性能提高会有很大的影响。BCP过程是SAT算法求解过程的核心部分(见图1 DPLL算法),占算法总开销高达80%～90%。且其本质是以单文字规则迭代传播,具有处理数据量庞大、运算简单的特点,适用于GPU处理。因此,我们以BCP过程作为基于GPU加速研究的主要目标,对于SAT求解算法中的变量选择和冲突分析与学习过程仍交由CPU来执行(如图2所示)。
现有的SAT求解器各模块设计相对独立,且模块之间提供了很好的交互接口,使得将GPU加速部分集成到不同的求解器中变得容易可行;将经由GPU加速的BCP部分集成到现有的SAT算法(如MiniSat等)中,即可实现更高效的SAT求解器。
2.1 GP_BCP设计
现有SAT算法的BCP过程设计的很精巧,每次迭代过程只需几千个时钟周期即可,使得通过GPU来获得可观的加速比很具有挑战性。因此,我们要充分利用GPU并行处理的特点充分利用程序的并行度,并使每次迭代过程尽量与问题规模无关。基于以上两点,我们提出了基于GPU并行加速的BCP过程:GP_BCP。
2.1.1 GP_BCP总体结构
(1) CPU&GPU通信
该模块通过PCIe总线从CPU接收选取的变量及其指派值,并将传播过程的推导结果及当前求解状态传送回CPU。为充分利用总线带宽、减少通信开销,我们基于数据预取思想,每次选择一组变量并对其进行指派,以数据块模式传送给GPU。
设置缓冲区存放该组指派,每次从中选取一个变量交由推理引擎进行处理直至缓冲区为空,再请求CPU重新选取一组变量传输。在推理过程中如产生冲突,只需作废缓冲区中指派而不需其他额外开销。通过以这种数据块为单位的集中传输模式,可充分利用传输带宽,减少通信次数与通信开销。
(2) 推理引擎
为变量与约束子句集合建立关联结构,每次给定决策变量及其指派后,推理引擎通过对关联的约束子句进行计算,从而推导出新的变量的指派值。
(3) 推导结果
经推理引擎处理后会导致两种结果状态:产生新的推导指派和遇到冲突。如产生新的推导指派,将其放入缓冲区以进行迭代传播,同时将其传送回CPU端更新变量指派值。这一步可通过异步传输实现,从而隐藏了传输开销。如产生冲突,需要终止推导过程并将冲突状态(包括冲突变量、冲突的原因)传送回CPU以进行冲突分析与学习。
2.1.2 基于GPU的新型推理引擎设计
推理引擎是GP_BCP的核心部分,决定算法的最终加速比。本部分主要介绍其实现的总体思想。
实际求解中,问题规模往往很大,甚至达到百万数量级变量数和千万数量级子句数。为提高传播效率,我们为变量与约束子句建立关联结构。每次选择新的变量进行指派时,只需对关联结构进行传播,使得传播开销与问题规模无关,极大减少了传播的开销。同时,对于关联结构中的每个约束子句,分布到GPU不同处理单元来并行处理。因此,每次迭代传播过程开销相当于串行算法中对一条约束子句进行传播的开销。
对每个关联子句运用推导规则:如果子句中只有一个文字未被指派且其他文字均被赋为FALSE,则剩下的文字就会被赋为TRUE,从而推出新的指派。基于Brook+编程模型为其设计计算核心,其伪代码如图4所示。
如图4代码,对每个与当前指派相关的子句进行推理,遍历子句中的每个文字。如果存在文字j,使得j之前和之后的文字均被指派为false,且j为undefined(未被指派),则可推导出新的指派j为true;如果子句中所有文字均已被指派为false,则产生冲突,返回冲突状态。
2.2 GP_BCP优化
在2.1节(2)的推理引擎设计中,通过将变量相关的子句分配到GPU多个处理单元上传播而实现并行,但并行度不够充分。为提高数据的并行度,将每次指派后产生的新的指派分组,以组为单位进行并行化传播。
每次选择变量进行指派并根据约束子句集进行传播过程中,由单文字蕴含规则可推导出新的指派值,称为推导指派。一般来说,问题规模越大,每次传播产生的推导指派的数目越多。经过对不同的求解实例进行测试,我们发现每次传播产生的推导指派数目最多可达到103的数量级。由于由同一个变量指派经传播派生出的推导指派,具有相同的decision level(决策级别,由同一指派经传播推导出的变量决策级别相同),对其同时传播不会更改程序的推导规则,对于冲突分析过程及回退不会产生影响。因此可将多个指派的传播过程并行执行,充分利用GPU并行处理单元,既可增大程序的并行度获取更好的加速比,又减少了GPU与CPU通信次数及通信开销。
对于不同的求解问题,与变量关联的子句规模不同,决定着如何将推导指派进行分组并行传播,决定着GP_BCP优化过程的设计。对于这个问题,可在求解前通过CPU进行预处理得到与变量关联的子句规模,将结果以参数形式传入GPU作为优化尺度,以设计实现GP_BCP的并行优化。如图5给出优化后的计算核心伪代码。
kernel表达了优化后的程序语义。在具体实现时,基于编程模型特点,可通过数据并行展开而使循环展开,从而分布到GPU每个处理单元执行inference(assignment i, attached clause j),从而实现多个指派并行传播。
3 性能评测
测试平台为4G内存、3G主频的HP服务器与AMD RadeonTM4890GPU。我们通过编写测试样本生成程序,随机生成k-SAT测试用例(3≪k≪6)[18]来对优化前后的GP_BCP进行测试,并选取高效的串行SAT求解算法MiniSat进行对比测试结果和分析。
对随机生成的k-SAT测试样本,来对初始的GP_BCP过程进行测试,相对于MiniSat的BCP过程,不同的测试样本取得了不同的加速比(见表1所示)。
由表1可见,对于不同的测试用例,GP_BCP所取得加速比并不相同。随着每个子句所含文字数的增加,加速效果呈增长趋势。这是由于随着测试用例中平均文字数目增加,每个文字相关的子句集增大,从而使每次选取变量指派并进行常量传播时程序的并行度增加,使得BCP过程的加速比更高。
以相同的测试用例对经过优化的GP_BCP进行测试,结果见表2所示。
由表2可知,经过优化GP_BCP过程更充分的利用了GPU并行处理能力,使GP_BCP数据处理并行度增大,相对于原有的GP_BCP,优化后的GP_BCP加速比明显增加。
图6所示给出GP_BCP过程优化前和优化后的性能提升情况的对比结果。
对比优化前后GP_BCP过程获得的加速比(图6),我们发现,经过优化后,3-SAT的GP_BCP性能提升幅度最大。而随着测试用例子句平均文字数的增加,性能提升幅度呈下降趋势。这是因为3-SAT平均文字数最小而使得优化前的GP_BCP数据并行度小,所以经优化后3-SAT的并行度大幅提高;而随着平均文字数增加,GP_BCP并行度增加,优化幅度减小。因此优化后的GP_BCP性能提高幅度呈下降趋势。
4 结 语
本文基于GPU实现了SAT算法的并行化加速研究,将SAT算法中BCP过程移植到GPU并行,实现了基于GPU的并行BCP过程GP_BCP,最后进行了测试。对于BCP过程并行化后,该部分取得了较好的加速效果,针对于不同的测试用例,取得了5.4～10.3的加速比。SAT算法本身及应用的特点决定了传播过程与冲突学习过程的频繁交互,而使得基于GPU的SAT算法整体性能提升有限。在以后的工作中,伴随着GPU单个处理单元处理能力的增强,期望能够将冲突学习及回退过程移植到GPU上处理,实现完全基于GPU的SAT算法,从而使效率有更大的提升。
GPU加速 第6篇
Ising模型[1]是研究铁磁相变的经典模型。一维情况下的Ising模型是没有相变的。1944年,Onsager得出了无外场下二维正方晶格Ising模型的相变临界点的解析解[2]。Ising模型具有丰富的内涵,使得它不仅在物理学领域有重要的应用,而且在生态物理,社会科学,经济学等领域也有着极为广泛的应用。
在过去的40多年里,得益于数值模拟、分析方法的发展,以及计算机性能的提高,数值模拟被广泛地应用于统计物理的研究中。在系统的临界点附近,传统的Ising模型Monte Carlo(MC)模拟算法的效率会急剧下降,出现所谓的“临界慢化”现象。随后出现了一些改进算法,例如Swendsen-Wang类团簇算法[3]通过同事处理一个集团内所有自旋大大地提高了算法的效率。虫子(worm)算法[4]则通过同事处理虫子内的所有自旋以提高算法效率。近年来,随着图形处理器(GPU)以及编译技术的迅猛发展,尤其是NVIDIA公司的计算统一设备架构(Compute Unified Device Architecture,CUDA)[5]的出现,使得“CPU+GPU”的通用计算模式得到越来越广泛的关注。
本文针对三角晶格Ising模型设计了一种基于GPU的MC算法,由CPU负责自旋的创建,GPU负责自旋的模拟,将自旋模拟中所要用到的能量计算、随机数产生以及状态更新等耗时运算从CPU转移到GPU中,降低了CPU的负荷,减少全系统自旋模拟所需的时间,提高了算法效率。
1 GPU及CUDA架构
GPU是针对高密度和并行计算设计的,不同于CPU架构,GPU由一组多处理器(Streaming Multiprocessors)和全局内存(Global Memory)组成。每个多处理器都含有8个处理核心,一组32位寄存器,只读的常量内存(Constant Memory)和纹理内存(Texture Memory)。在CUDA环境下,CPU被看作主机(Host),GPU为设备(Device),GPU作为CPU的协处理器(Coprocessor)使用。GPU中执行运算的最小单元是线程(Thread),一定数目的线程组成一个线程块(Block),多个具有相同线程数的线程块再构成一个网格(Grid)。在GPU中所有线程块以及线程块中的线程都是同时运行的,对线程的创建、管理、调度和执行指令由SIMT单元以warp为单位完成。
GPU中所有线程同时运行,执行相同代码但操作不同数据。GPU的存储单元以及GPU与CPU的数据传输带宽对程序效率有较大的影响。GPU中共享内存的访问速度较全局内存快,而GPU与CPU之间的数据传输又较GPU内部的数据传输慢,为了使程序更为优化,应最小化设备与主机之间的数据传输,同时应尽可能多地在共享内存进行数据的访问,以减少访问数据所需的时间。本文所用的GPU的主要参数见表1。
2 Ising模型MC算法的GPU实现
二维三角晶格Ising模型如图1所示,定义系统尺寸为L,即有LL个自旋粒子,处于二维三角晶格格点位置,每个粒子的自旋只能取向上或者向下两个态,并只考虑近邻自旋间的相互作用,模型的哈密顿量表示式(1)。
式(1)中J为与交换积分成正比的耦合常数,J>0代表铁磁体,J<0代表反铁磁体。第i个格点位置的自旋si的取值为+1或者-1,分别对应于自旋向上和向下;<i,j>表示最近邻自旋对。式(1)中右边第一项为自旋的相互作用,第二项代表外场的作用,H代表外磁场。以下仅讨论J>0,H=0的情形,系统的哈密顿量变为
根据Metropolis算法[6],自旋si翻转的概率为
P(si-si)=min[1,e-βΔE] (3)
式(3)中β=1/(kBT),T为绝对温度,kB为波尔兹曼常数,ΔE为自旋翻转前后系统的能量差。
传统的Metropolis算法对全系统执行一次扫描的方法如下:(1) 选取1个自旋计算其自旋翻转前后的能量差ΔE并计算翻转概率P;(2) 通过随机数发生器生成随机数r,若rP则将自旋翻转,否则保持自旋不变;(3) 重复(1)和(2)直到所有自旋均被扫描过。
如图1所示,对于三角晶格Ising模型,自旋si是否翻转仅取决于si及其6个最近邻的状态,这表明了同时模拟多个自旋的可行性。而并行运算正是GPU的强项,因此采用GPU对系统进行模拟。在对系统的自旋进行并行模拟时,若同时模拟所有格点上的自旋,每个格点上的自旋状态会被其邻居读取(即读内存)。与此同时,该自旋又通过读取邻居自旋状态以更新自己的状态(即写内存),这样就会发生同时读写内存的冲突。为了避免这种情况,可以将晶格分割成类似于棋盘的周期分布,如图2所示,同时模拟相同编号的自旋,这样既能得到较高的运算效率又不会产生读写冲突。
在GPU程序中,对于LL的系统,首先将晶格分成L/2个2L的晶格块,分别指派给L/2个线程块。再将晶格块中的2L个格点分成L/2个22的子格并指派给线程块中的L/2个线程, 即每个线程需模拟一个22子格上的自旋。为了方便描述,22子格中的自旋依次编号为0,1,2,3,各线程同时模拟对应的小块中相同编号的自旋,避免了对自旋同时读写的问题。
以88三角晶格Ising模型为例,如图2所示,CUDA程序的具体步骤如下:
(1) 执行内核,指派4个独立的线程块,每个线程块中运行4个线程;
(2) 在16个线程并行地生成随机数储存在共享内存中;
(3) 各线程并行地从全局内存中读入对应子格中0号位置的自旋及其6个最近邻的自旋状态、计算ΔE、执行Metropolis算法;
(4) 退出所有线程以确保所有小块中的0号位置自旋的模拟均已完成;
(5) 依次对1,2,3号自旋执行步骤(1)～(4)。
3GPU算法的效率及有效性
为了验证GPU算法的效率,需要对GPU算法所用时间与CPU算法进行比较。本文在Intel Core i7 920和Tesla C1060的平台上分别模拟了L=16～1 024的系统,为了比较的准确性,没有计算统计系统能量、磁化率等基本物理量所用的时间。将全系统执行一次MC模拟所需时间除以系统的自旋总数作为单个自旋的平均翻转时间,结果表明,仅使用CPU模拟时,单个自旋翻转时间与系统尺寸无关,基本稳定在50 ns;而在GPU模拟中,随着系统尺寸L的增大,GPU同时处理的自旋数也相应地增加,从而使单个自旋的翻转时间随着L的增加迅速减小,在L=1 024时仅为0.9 ns。可以直观地用加速比(speedup)[7]作为GPU加速效果的衡量标准,其定义如下:
图3给出了不同尺寸的三角晶格下,GPU和单CPU算法对所有自旋进行1 000次模拟所需时间以及GPU的加速比。可以看到,对于小系统,由于不能有效地使用GPU的线程数,其加速比甚至小于1。但是,随着尺寸的增加,GPU的效率显著上升,尤其是对于L=1 024的系统,GPU可以实现约69倍的加速。对于更大的系统,由于GPU核心数的限制,当GPU中的线程利用率达到饱和时,GPU的加速性能将会达到瓶颈。
在看到GPU所能带来的时间成本节约的同时,还需验证GPU算法的有效性。仅比较GPU和CPU这两种算法下得到的系统的能量E和磁化率m是远远不够的,较为合理的办法是比较这两种算法下系统的临界点和临界指数。
Binder累积量[8]是寻找临界点的一个重要工具,其定义如式(5)。
式(5)中m为系统的磁化强度,L为系统尺寸,<>表示平均值。QL在临界点的值与系统尺寸无关, 因此可以使用式(6)来确定临界点:
式(6)中L′和L是两个不同的系统尺寸。图4所示为不同尺寸系统的QL随β的变化,图中曲线的交点即为三角晶格Ising模型的临界点。
为了更准确地获得三角晶格的临界点及临界指数,将获得的数据使用式(7)进行拟合。
QL=Q0+a1(β-βc)Lyt+a2(β-βc)L2yt+b1Ly1+b2Ly2+c1(β-βc)Lyt+y1 (7)
式(7)中Q0,ai,bi,ci,βc均为未知数。在对QL的拟合中取y1=-1.75,y2=-2.0,得到βc=0.274 66(1),Q0 =0.86(3),yt=1.01(2),βc与理论值ln(3)/4高度吻合。采用模拟获得的参数,QL的数据经过有限尺度修正后的结果如图5所示。
为了获得yh,还需对磁化率χ进行了拟合,磁化率χ的定义为:
χ=β(<m2>-<m>2) (8)
所用拟合公式如下:
χ=χ0+L2yh-2[a0+a1(β-βc)Lyt+a2(β-βc)2L2yt+b1Ly1+b2Ly2+c1(β-βc)Lyt+y1] (9)
得出yh = 1.875 6(3),与理论值15/8的结果吻合。在确定的β=0.274 653上拟合了式(9)。在临界点附近,χ的数据如图6所示,图6中近似的直线说明χ的行为满足2yh-2的标度率。
4结论
本文在CUDA环境下,使用GPU作为CPU的协处理器,实现了对二维三角晶格Ising模型的模拟。结果表明,在小尺寸情况下,由于GPU线程的利用率过低,GPU的运行效率甚至不如单CPU的效率。但是,随着系统尺寸的增加,GPU的运行效率快速增长,尤其是对于10241024的系统,使用Tesla C1060,可以获得Intel(R) Core(TM) i7 CPU (2.67 GHz) 69倍的加速比。同时,应用优化后的GPU程序,对三角晶格Ising模型的临界行为进行了探讨,获得了高精度的临界点βc=0.274 66(1)和临界指数yt=1.01(2),yh=1.875 6(3),结果与理论值βc=ln(3)/4,yt =1,yh =15/8吻合,验证了GPU加速三角晶格Ising模型MC模拟的有效性。
摘要：在分析传统Monte Carlo算法的基础上,针对三角晶格Ising模型提出了一种基于GPU的并行模拟方法,大大提高了算法的效率。对1 024×1 024的模型,实现了69倍的加速比。通过该算法所得数据分析模型的临界行为,获得了高精度的临界点βc=0.274 66(1)和临界指数yt=1.01(2),yh=1.875 6(3)。
关键词：GPU,CUDA,Ising模型,Monte Carlo模拟,临界行为
参考文献
[1] Ising E.Beitrag zur theorie des ferromagnetismus. Z Phys, 1925; 31(1):253—258
[2] Onsager L.Crystal statistics.I:a two-dimensional model with an or-der-disorder transition.Phys Rev,1944;65(3—4):117—149
[3] Swendsen R H, Wang J S. Nonuniversal critical dynamics in Monte Carlo simulations. Phys Rev Lett, 1987; 58:86—88
[4] Deng Y J, Garoni T M, Sokal A D. Dynamic critical behavior of the worm algorithm for the Ising model. Phys Rev Lett, 2007; 99:110601
[5] NVIDIA. NVIDIA CUDA C Programming Guide. http://developer.download.nvidia.com/compute/cuda/3_1/toolkit/docs/NVIDIA_CUDA_C_ProgrammingGuide_3.1.pdf. 2010
[6] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equations of state calculations by fast computing machines. J Chem Phys, 1953; 21(6):1087—1092
[7] Preis T,Virnau P,Paul W,et al.GPU accelerated Monte Carlo sim-ulation of the 2D and 3D Ising model.J Comput Phys,2009;228(12):4468—4477
GPU加速 第7篇
关键词：图像融合,GPU加速,几何校正
1 引言
图像融合 (Image Fusion) 是将数张有重叠部分的图像融合成一幅大型的无缝高分辨率无重影的图像的技术。是指将多源信道所采集到的关于同一目标的图像数据经过图像处理, 最大限度的提取各自信道中的有利信息, 最后综合成高质量的图像。它是数据融合的一个重要分支, 随着应用技术的发展, 图像融合开始引起人们的关注, 特别是在最近二十几年里, 国内外在图像融合的不同层次上开展了大量的研究, 取得了很多宝贵的研究成果, 涌现出大量的文献, 如基于小波包和边缘特征的遥感图像融合算法[1], 基于边缘统计特征的遥感图像融合改进方法[2], 一种基于相似度的HIS变换融合算法[3], 节点内无冗余图像合成方法[4], 基于小波变换和邻域特征的多聚焦图像融合算法[5]等。在众多算法中, 处理平滑过渡, 提高图像质量的较多, 对视频图像融合尚不能做到差强人意, 仍有很多问题有待解决, 比如, 在如何改善图像融合算法的运算效率、稳定性、可靠性等方面, 都值得深入研究, 还有在多投影显示系统应用中, 需要很多图像并行计算操作, 而现在微机处理器的并行计算能力有限, 特别是在处理超高清图像播放时, 算法的运算效率的问题仍然没有得到很好的解决。
图形处理器 (Graphics Processing Unit, GPU) 是显卡的计算核心, 是一种专用的图形渲染设备, 随着GPU技术, 特别是可编程图形处理器技术 (PGPU) 的飞速发展, GPU高度并行化架构和可编程GPU处理单元, 使得GPU适用于处理高度并行的通用科学计算, 并在图形图像中的应用越来越广泛[6], 使用GPU进行通用计算和图形计算已成为一种发展趋势;GPU在并行性能方面远高于传统的多核CPU并行计算, 在图形计算方面, 基于GPU的计算获得了更加广泛的应用, 例如大规模视频拼接、快速光线投射、图像压缩与绘制等。因此, GPU在图形图像研究领域的应用也越来越多。本篇文章利用GPU的高度并行计算能力, 解决了多投影系统算法的运算效率的问题, 实现了超高清图像的多投影实时播放系统。
2 基于GPU加速的多投影融合系统
本投影系统实现是基于局域网的集群式显示绘制系统, 采用Client-Server模型, 系统结构如图1所示, 主要包括:主控单元、校正核心模块、显示控制模块
图形图像处理及服务器模块。其中主控单元客户端应用程序向用户提供投影阵列布置、多投影显示效果参数设置、投影区域设定、切分纹理设定、几何校正和颜色校正实时编辑等功能, 它实时调用运行在服务器端的校正核心模块完成几何校正与颜色校正的功能。对图像进行相应的几何变形和颜色参数进行实时修改和显示, 将测试结果保存到各个显示节点。测试结果将应用在多投影实时播放器, 对视频图像的颜色与几何实时校正, 完成图像显示工作。本多投影融合系统最大的特点是在核心几何校正与颜色校正整合过程中采用了GPU加速技术, 提高图像处理的效率, 完成了超高清图像的实时投影播放。
3 基于GPU加速半自动几何校正技术
多投影拼接是多投影融合显示系统中核心技术, 它是将数张有重叠的图像拼成一幅大型的无缝高分辨率的图像, 其中图像配准与图像融合是其关键, 多投影播放图像的质量主要取决于这两项技术。图像配准又称为几何校正, 是投影融合的必要前提条件, 融合的清晰度, 图像有没有重影和融合带受图像几何校正精度影响较大。在投影参数配置阶段的几何校正操作是根据校正图案由人工操作来完成的, 校正图案是根据投影墙的特点, 预先使每个投影仪投射一种易于识别和定位的特定模式图形。典型的校正图案有黑白格、网格线, 校正图案的好坏和人为操作的因素直接影响几何校正的精确度。鉴于以上问题, 本文提出一种基于GPU加速的半自动几何校正技术, 它具有稳定性强, 精度高, 运算快等优点, 并且它具有对配置后对由外界引起投影仪的位置移动进行几何校正参数自动补偿的功能。该算法的具体步骤如3.1节所示, 以两台投影仪 (L, R) 几何校正为例:
3.1 半自动几何校正
1.校正投影仪投射校正图案;
2.根据校正图案进行人为校正, 并计算出相应的投影几何变换矩阵参数HL与HR;
3.投影仪L单独投校正图像, 同时C启动拍下校正图像IL;
4.投影仪R单独投校正图像, 同时C启动拍下校正图像IR;
5.根据后面的图像配准算法, 精配准IL与IR;
6.根据配准结果校正HL。
3.2 基于GPU加速的图像配准算法
由于在多投影融合系统中要配准的图像的来源相同, 所以在配准过程中选用的相似度量函数是相关系数, 其定义如式1所示。
为了便于公式 (1) 的快速计算, 它可以简化为公式 (2)
其中g为模板图像 (大小为mn) , 大图像为S (大小为MN) , 用Sx, y表示S中以 (x, y) 为左上角点与g大小相同的子块.将两个相同尺寸的图像的相关系数记为ρ (x, y) .
分析式 (2) 可知, 其中.S与g相乘并求和的部分的算法复杂度依然是O (Mx Nxmxn) , 记Sum (x, y) 为Sx, y与g相乘并求和的结果, 共要计算MN次Sum (x, y) , 它成为整个匹配算法计算的瓶颈。本文采用GPU并行运算, 利用CUDA编程加以优化加速。
CUDA是NVIDIA推出的运算平台, CUDA™是一种由NVIDIA推出的通用并行计算架构, 该架构使GPU能够解决复杂的计算问题, 它使用C语言为基础.在CUDA架构下, 显示芯片执行时的最小单位是thread, 数个thread可以组成一个block。
4 基于GPU加速的图像融合算法
多投影显示系统图像显示质量主要决定于融合带处理的质量, 对融合带中的每个画面的显示都要经过几何校正与亮度校正的渲染过程, 其中几何校正主要是根据第3部分计算得到的几何变换参数确定每台投影仪所投图像在全局画布中相对位置。多投影颜色校正工具对投影仪颜色的校正过程包括对各个投影仪的色彩平衡、亮度区间、融合曲线、Gamma曲线进行调试。首先要调整各投影仪之间的融合亮带, 然后通过指定亮度衰减曲线对颜色全通道 (或红、绿、蓝三通道分别调整) 进行Alpha亮度衰减调整。根据多投影显示墙的画面效果, 选择最好的一组调试参数进行保存, 从而确定图像重合带的融合所需的全部参数值。整个显示画面当中, 特别是融合带及其周围选定扩展区域内的每一个像素的渲染过程都要进行以上两个过程的校对工作。
本系统中将几何校正与亮度校正渲染计算过程合并并采用GPU进行加速, 从而完成整个显示画面的渲染工作。具体实现过程分为GPU显存访问的优化与融合算法的CUDA并行实现。
4.1 GPU显存访问的优化
在GPU中计算, 必须先把内存中的数据拷贝到显存中, 然后访问显存进行计算, 而访问显存需要消耗400~600个时钟周期时钟周期, 而所以设计时应该尽量减少对显存的访问次数.CUDA在显卡端提供了一些可以高速缓存的寄存器, 寄存器 (register) 等快速存储部件访存一次只需要4个左右的时钟周期。优化时充分利用GPU显存中的多层次存储部件, 发挥快速存储部件读取数据优势, 最大化提高执行性能。
4.2 CUDA并行的融合算法原理
将图像显示渲染过程中几何校正与颜色校正的总的计算, 根据所用显卡硬件的每个线程块的最多拥有的纯种数, 比如每个流处理单元SP最多可同时执行256个线程, 就block设置为1616二维形式, 从而使得处理每个block时, 使每个SP单元能够满负载执行, 最大限度地提高资源利用率, 然后根据总的显示画面 (大小为widthheight) , 将整个渲染任务分为宽为 (width+16-1) /16, 高为 (height+16-1) /16的任务节点网格, 每个网格节点为一个block, 当然整个任务所分的任务节点个数, 与相应的显卡参数有关, 如果每流处理单元SP最多可同时执行block.yblock.x线程, 则总任务可分为宽为 (width+block.x-1) /block.x, 高为 (height+block.y-1) /block.y的任务节点网格。
5 实验
下面是通过几何校正和颜色校正获得的图像融合效果, 分别设置了红、绿、蓝三种颜色, 如下图:
6 总结
定量地评价图像融合性能是一项重要而复杂的工作。实验表明, 通过几何校正和颜色校正获得的图像融合效果, 其颜色过度更加平稳, 中间的色带几乎完全消隐掉, 其算法适用于不同颜色, 且融合带均消隐掉, 充分显示了此算法的优越性。
7 结论
多源图像的融合层次多种多样, 融合技术千差万别, 然而, 基于GPU加速的图像融合技术很大程度上解决了传统融合算法在视频播放低的局限性, 本算法在提取原来算法优越性的基础上, 加入了GPU技术, 克服了原来算法的不足, 使视频融合技术真正适用于任何普清及高清视频的播放, 实验表明, 该算法充分利用GPU的多发核并发计算的特性, 有效的提升播放速率, 当然, 我们下一步除了要提高播放速率外, 还要进一步优化算法性能, 使其达到更理想的融合效果。
参考文献
[1]曾宇燕, 何建农.基于小波包和边缘特征的遥感图像融合算法[J].计算机应用, 2011.31 (10) .
[2]曾宇燕, 何建农.基于边缘统计特征的遥感图像融合改进方法[2][J].计算机工程与应用, 2013.49 (3) .
[3]王晓艳, 刘勇.一种基于相似度的HIS变换融合算法[J].遥感技术与应用, 2011.26 (5) .
[4]刘华海, 王攀, 李思昆.节点内无冗余图像合成方法[J].计算辅助设计与图形学学报, 2013.25 (5) .
[5]基于小波变换和邻域特征的多聚焦图像融合算法[J].西北工业大学学报, 2011.29 (3) .
[6]宋暨, 周松涛.基于GPU并行计算的图像快速匹配[J].湖北民族学院学报 (自然科学版) , 2011.29 (3) .
[7]詹洪陈, 袁杰.图像处理的GPU加速技术研究[J].现代电子技术, 2012.35 (20) .
GPU加速 第8篇
蚂蚁在觅食路径上会释放一种特殊的分泌物—“信息素”,随着时间流逝,信息素会挥发,后面的蚂蚁根据路径上的信息素浓度,选择信息素多的路径去觅食,这样便形成了一个正反馈机制。在整个寻径过程中,虽然单只蚂蚁的选择能力有限,但它们的行为具有非常高的自组织性,相互之间交换路径,最终寻找到最优路径。受到蚁群系统信息共享机制的启发,意大利学者Dorigo M于1992年首次系统提出了蚁群算法,并成功地将该方法应用到求解TSP问题中[1]。该算法是启发式搜算算法的一种,采用了分布式并行计算机制,易于与其他方法结合,具有强的鲁棒性。同时,相对于其他搜算法,对初始路线要求不高,在搜索过程中不需要人工调整。研究表明,蚁群算法是解决TSP问题有效的算法之一,因此也成为近年来的研究热点。
近年来,基于GPU(图像处理器)的大规模通用并行计算,大大提高了计算机图形处理的效率[2]。GPU的高速浮点运算能力、并行计算和可编程功能也为通用计算提供了良好的并行计算平台,同时也为蚁群算法的高速并行实现提供了可能。NVIDIA公司的统一计算设备架构(CUDA),为研究人员利用CPU进行数据并行处理提供了便捷的手段[3]。通过GPU加速并行算法,是蚁群算法并行化、提高算法执行速度的有效途径之一。
1 蚁群算法基本原理
蚁群算法受蚁群行为和反应特征的启发而得来的。觅食的蚂蚁在走过的路上释放信息素,其他觅食蚂蚁将沿着留有信息素的路径在相同位置找到食物。因此,信息素可用来实现间接通信的目的[4]。基于信息素的城市探知转移概率公式:
这里τij(t)是城市i与j之间的信息素浓度,ηij(t)是城市i与j之间的可见距离,dij是城市i与j之间的距离,表述为公式(2):
Pkij(t)是选择路径时的转移概率,α是信息素系数,β是绝对距离系数,分别用来控制τij(t)和ηij(t),t表示时间,当α设置的合适,有利于选择较好的行动路径;β值越大,说明探知能力越强。
蚂蚁在觅食过程中,将在通过的路径上释放信息素去记录路径,同时其他蚂蚁将优先选择该条路径,通过的蚂蚁越多,该条路径上的信息素越浓。同时,信息素也会因环境因素而挥发。因此,信息素更新公式可以用式(3)表示:
其中,ρ表示信息素挥发系数,其值介于0和1之间。如果ρ太小(接近于0),则结果不会收敛,因为信息素消失并且不会随时间积累。如果ρ太大(接近于1),则结果会过早收敛。Δτij是城市i与j之间的信息素,由式4表示。
Q是蚂蚁隐含的信息素。当蚂蚁经过路径时,τij(t+1)等于ρτij(t)加上Q。如果路径上没有蚂蚁通过,τij(t+1)等于ρτij(t),随着时间的流逝,τij(t)将变为0。
蚁群算法的流程如图1所示。首先,设置蚂蚁的基础参数和起始位置,其次,计算选择城市的转移概率并遍历城市。当第一次遍历结束,计算该次遍历经过的总路径长度,同时更新路径上的信息素,再次将蚂蚁放置起点位置,进行下一轮搜索食物的旅程,最后每次旅程将进行比较,并找出最短路径,输出结果。
计算路径转移概率、选择转移路径时,通过如下方式进行遍历城市。路径示意如图2所示,城市S是原点,城市A到G为允许转移的城市。PSA、PSB、PSC、PSD、PSE(t)、PSF、PSG分别为城市S到城市A~G的转移概率。给出随机数R、求和数sum、转移概率P,其中,sum=PSA+PSB+PSC+PSD+PSE+PSF+PSG。首先,P=PSA/sum,比较P和R,若P大于R,则选择路径A;否则,P=(PSA+PSB)/sum,比较P和R,若P大于R,则选择路径B;否则,P=(PSA+PSB+PSC)/sum,比较P和R,若P大于R,则选择路径C;依次类推,通过这种方式选择下一个遍历城市。
2 基于GPU并行的蚁群算法
通过对CUDA和GPU并行程序设计的研究,将基于蚁群算法的TSP求解优化路径问题转化为GPU中的多线程的并行处理过程,充分利用GPU的高速浮点运算和并行计算的特性,以提高算法的速度。
在蚂蚁遍历城市的过程中,每只蚂蚁独立地建立自己的遍历路径,对蚂蚁来说,自然地存在并行性。根据CUDA的编程模型[5],我们很自然地会将每只蚂蚁个体映射到CUDA的一个线程上,用线程来模拟你蚂蚁个体,每个线程完成一只蚂蚁的城市周游回路。
设蚂蚁个数为M,城市个数为N,CUDA中Block个数为4,蚁群算法的GPU并行化模型中,总线程个数为M,每个Block中线程个数为M/4,如图3所示。在模型中,将路径状态初始化、路径转移概率计算、路径长度计算、更新信息素定义为CUDA函数。首先,M个线程被分配在4个Block中,同时启动,计算各自的城市转移概率,寻找转移城市;其次,由N个线程分成N个Block,计算路径上信息素增量;再次,用一个线程计算路径长度和最优路径;最后N个线程分成N个Block,更信路径信息素。
3 实验与分析
本文采用TSPLIB中的算例,实验环境为Inte(lR)Core(TM)i7-5820K@3.30GHz CPU,NVIDIA Ge Fore GTX 960显卡的计算机,实验软件环境为Win7操作系统、CUDA7.0版本,在VS2010环境下编程。实验数据和相关参数如表1所示。
分别采用CPU和GPU算法,对每个问题进行多次实验。实验数据表明本文算法获得了较高的加速比,表2数据显示了对不同问题进行300次迭代时,本文算法对基本蚁群算法的加速比。
实验结果表明,采用GPU并行化蚁群算法取得了良好的加速效果,当蚁群规模由256增大到1024、城市规模由21增大到76时,加速倍数由3.0增大到10.75。数据显示,问题规模越大,蚁群个体和城市规模越大,加速效果越好,对于更大规模的复杂问题,基于GPU的并行化蚁群算法将取得更高的加速比。
4 结束语
本文研究了基本蚁群算法的原理,并基于GPU和CUDA高速并行计算模型,利用GPU在CUDA7.0环境下,对蚁群算法进行化加速设计,实验结果表明,该方法取得了良好的加速效果,当蚁群规模增大时,加速倍大幅度提高。数据显示,蚁群个体和城市规模越大,加速效果越好。
摘要：蚁群算法是求解旅行商问题的有效方法之一,但是随着蚁群规模和城市规模的增大,算法的运行速度将大大降低,本文利用GPU在CUDA7.0环境下,对蚁群算法进行化加速设计,实验结果表明,该方法取得了良好的加速效果,当蚁群规模增大时,加速倍大幅度提高。数据显示,蚁群个体和城市规模越大,加速效果越好。
关键词：蚁群算法,GPU,CUDA
参考文献
[1]宗德才,王康康,丁勇.蚁群算法求解旅行商问题综述[J].计算机与数字工程,2014(11).
[2]占正锋,李戈,张学贺,伊旭悦.基于CUDA的图像预处理并行化研究[J].智能工程,2014.
[3]李建明,胡祥培,庞占龙,钱昆明.一种基于GPU加速的细粒度并行蚁群算法[J].控制与决策,2009(8).
[4]Shi-An Li,Min-Hao Yang,Chung-Wei Weng,Yi-Hong Chen,Chia-Hung Lo,Ching-Chang Wong.Ant Colony Optimizationalgorithm design and its FPGA implemention[J].IEEE International Symposium on Intelligent Signal Processing and Commu-nication Systems,2012.
GPU加速范文
声明：除非特别标注，否则均为本站原创文章，转载时请以链接形式注明文章出处。如若本站内容侵犯了原著者的合法权益，可联系本站删除。
