基于CUDA的汇流分析并行算法的研究与实现

　　摘 要:针对基于数字高程模型(dem)生成流域等流时线的快速运算问题,提出了一种基于统一设备计算架构(cuda)平台同时可发挥图形处理器(gpu)并行运算特性的汇流分析的快速并行算法。采用改进后的归并排序算法进行数据排序及新的内存分配策略和改进的并行算法进行汇流分析。用该并行算法和cpu上的串行算法, 对生成基于dem的等流时线运算时间和矩阵乘法运算时间进行分析验证。实验结果表明,基于cuda的汇流分析并行算法能提高系统的计算效率,具有较好的效果。
　　关键词:并行计算;图形处理器;统一设备计算架构;汇流分析;数字高程模型
　　　　doi:.1001-3695.2010.07.011
　　
　　research and realization in parallel algorithm ofconfluence analysis based on cuda
　　
　　zhao xiang-hui1,2, miao qing 1,2, fu zhong-liang1,2, su chang1,2, li xin1,2
　　(u institute of computer applications, chinese academy of sciences, chengdu610041, china;te school, chinese academy of sciences, beijing 100049, china)
　　
　　abstract:aiming at the fast parallel computing of generating isochrones of watersheds that based on digital elevation model(dem),this paper proposed a fast parallel algorithm of confluence analysis based on compute unified device architecture (cuda) platform that could use parallel computing of graphic processing unit(gpu).carried data sorting out by using the improved merge sorting algorithm,adopted the new memory allocation strategy,carried confluence analysis out by using the improved parallel computing algorithm. using the presented parallel algorithm and the serial algorithm based on cpu to analyze and verify the time consumed when generating isochrones of watersheds based on dem and executing matrix multiplication. the experiments results illustrate that this parallel algorithm of confluence analysis based on cuda can improve the computational efficiency of the system and have a better effect.
　　key words:parallel computing; graphic processing unit(gpu); compute unified device architecture(cuda); confluence analysis; digital elevation model(dem)
　　cuda是nvidia公司提供的gpu用于通用计算的开发环境,是一个全新的软硬件架构,可以将gpu视为一个并行数据计算的设备,对进行的计算进行分配和管理。WWw.lw881.com图形显示卡,尤其是高端的显卡,现在几乎是单机必配的重要硬件。其核心部件gpu采用 simd体系结构,大量的计算部件使得它很适合处理高数据量的并行运算[1~3]。基于cuda平台,利用gpu强大的并行处理能力来加速通用科学计算及普通应用程序是非常好的应用趋势,已经成功应用到诸如金融风险管理、媒体图像以及科学研究[4]等领域。
　　本文结合科研项目“基于3s技术的自然灾害智能预测平台的研发与应用”中相关的关键信息技术来开展基于cuda的汇流分析并行算法的研究与实现。在自然灾害的山洪汇流分析中,项目中要绘制的等流时线是基于dem高程图绘制的[5,6]。若采用传统的基于cpu的串行算法,高分辨率的dem图在带来更好的计算效果的同时,也大大地增加了计算量,使得计算时间大幅度增加,降低了计算效率。例如作为四川省试点地区的彭州甘溪沟流域,其流域面积11.1 km2,沟道长度为6.83 km,若使用80 m作为单位栅格边长来划分流域,则流域可被划分为1 805个有效栅格单位,若采用高精度的划分方法,以20 m作为栅格边长来划分流域,则流域被划分为28 900个有效栅格,在能获得更高精度dem图的情况下,栅格的数量会更大。这样,在获得更精确结果的同时,大大牺牲了运算效率。传统的方法难以满足高数据量的计算需求,本文算法的计算核心是矩阵相乘计算,需要花费大量的时间,gpu在计算矩阵相乘时,比cpu要快得多,因此可巧妙地利用基于cuda平台的gpu强大的并行处理能力来加速计算速度[7,8]。本文提出了一种基于cuda平台可发挥gpu并行运算特性的快速并行算法。

1 cuda并行计算的引入
　　cuda表示统一计算设备架构,使gpu能够解决复杂的计算问题。它是一个并行编程模型和一个软件编程环境,它主要是为了帮助广大的程序员来更好地开发平滑扩展的并行程序。当前主流gpu产品在浮点运算性能上都超过了1 tflops[9],而intel公司最新的四核处理器core2 i7也仅达到70 gflops[10]。cuda模型中,gpu相当于cpu的协处理器,它能并行执行非常多的线程,处理大规模的并行运算。整体来说,gpu芯片类似于流处理器,适合一次进行大量相同的工作;cpu则比较有弹性,能同时处理变化较多的工作[3]。图1表明,gpu具备了超高的计算能力,在多个处理核心和高存储器带宽的配合下,最新的gpu成为了图形和非图形处理的超强工具。cuda的软件堆栈[11]采取如图2所示的多层组成方式,以此实现高性能计算。它允许定义一种叫kernel的函数扩展,当一个kernel函数被调用时,会有n个不同的cuda线程在gpu中并行地执行相同的程序,并行执行的线程数量用扩展的“〈〈〈〉〉〉”来指定。例如functionongpu〈〈〈4,256〉〉〉(…),该函数一旦成功调用,便共有4×256=1 024个独立的并行线程来执行该函数的内容。其中4为块(block)的数量,256为每块内的线程(thread)的数量。
2 基于dem的流域等流时线绘制算法
　　
　　为了解决山洪灾害预测问题,通过建立山洪形成过程的数学模型,实现对山洪预测的数值模拟。等流时线是一种经典的流域汇流曲线[6],根据对基于数字高程模型的流域等流时线绘制算法的描述[12],在整个计算和绘制的过程中,主要经过五个大步骤的计算:a)对dem图作填洼处理;b)根据高程对栅格进行排序;c)判断栅格水流方向;d)计算单个栅格汇流时间;e)计算流域汇流时间。其中:步骤a)填洼处理可以由gis软件,本文使用supermap deskpro2008软件来完成;步骤b)对栅格按高程排序,一般排序算法的复杂度为o(n log 2n ),对于大量元素的排序,较为耗费系统计算资源;步骤c)判断栅格水流方向,需要与栅格周围的八个栅格进行比较;步骤d)计算单个栅格汇流时间;步骤e)从dem值最小的栅格开始计算流域汇流时间。若采用传统的基于cpu的串行算法, 这些运算都需要对栅格进行多次的遍历,在栅格划分较多的高精度dem图上,就会成为整个算法的瓶颈。
　　观察分析以上几步可发现,在栅格运算的过程中,排序可以使用优化的并行排序算法来提高效率,步骤c)和d)的运算,由于其栅格之间的相关性很低,计算水流方向时,只需要考虑与之相关的八个栅格;计算单个栅格的汇流时间时几乎是独立运算。这两步运算符合simd(single instruction multiple data,单指令多数据流)指令系统,宜于利用并行运算将其实现。按照规格网格法,把系统所需的dem图表示成高程矩阵,在计算机中以二维矩阵的方式存放。流域外的空白地区的高程值被设为-999。所以对dem图的操作,实际上就是对高程矩阵的运算。cpu的流式指令体系结构在提高矩阵运算速度上有着很大的限制。gpu(图形处理器)采用单指令多数据流体系结构,大量的计算部件使得它很适合处理矩阵运算。因此,本文将以上几步的算法进行优化,使其适合并行计算,并基于cuda平台利用gpu的通用计算功能,使用改进的并行算法在单机上进行并行实现,最大化利用电脑硬件的计算能力。
　　3 基于cuda汇流分析的并行算法
　　3.1 算法实现流程
　　本文使用改进后的归并排序算法对栅格进行排序。将第2章中介绍的相关算法进行优化,使其适合并行计算。该算法在gpu上的具体实现(以四川彭州甘溪沟流域的dem图为例来对上述算法并行部分进行验证),主要流程如下:
　　a)将dem转换为高程矩阵。
　　b)使用改进后的归并排序算法,按照高程值对栅格进行排序。对高程矩阵按照数据量将数据分块并排序,将结果送至共享内存并分配空间,同步线程后,对奇偶块进行merge运算。
　　c)计算栅格水流方向和经流时间。
　　d)绘制等流时线进行汇流展示。
　　3.2 算法及其应用分析
　　3.2.1 按照高程值对栅格进行排序
　　加利福尼亚州大学的satish等人[13]提出了一种在cuda上进行归并排序的算法,并提出通过预先抽取数据采样,对待排序序列根据采样值进行分块来提高归并算法中merge运算的效率。同时,采样以256个待排元素为间隔,来保证划分出的每个子块最多含有256个元素。这样使在merge运算时划分出的子块能够快速地在gpu上基于其共享内存来运算。
　　
　　在计算等流时线前需要对栅格进行排序,以获得高程值最低的栅格作为整个流域汇流的出口。在排序算法中,由于多路归并排序,将数据分组排序后再合并,适合并行计算,在此处使用改进后的归并排序。首先将dem对应的高程矩阵看做一维数组,然后按照数据量将数据分成k个块,每块分别送进cuda的一个block,使用bitonic对其进行排序,最后将结果送至共享内存中,同步线程后,对奇偶块进行merge运算。从国内外已有的研究来看,排序运算元素之间并不能达到高的独立运算,元素之间仍有一定的相关性,并行时加速比不是很高。因此,并行排序运算虽然有一定的效率提高,但并不是整个运算加速的重点。
　　
　　3.2.2 计算栅格水流方向和经流时间
　　假设高程矩阵是一个m×n的矩阵,根据汇流分析的相关知识,求取栅格水流方向,只需计算该栅格与相邻的八个栅格的距离权落差,在计算单个栅格的流经时间时,需要使用公式:
　　ti=nl/[kv(si)b](1)
　　式中各个参数的意义请参见文献[12],只有参数si与该栅格的流向栅格有关,这两个运算之间的相关性很高,可以在一个运算过程中完成。在这两个运算过程中,栅格与栅格之间的相关度较低,可以方便地实现高并行运算。
　　1)高程矩阵的分块和cuda中内存的分配
　　该算法的整个并行计算过程中,尽量降低内存带宽的使用是提高性能的重要方式。而该算法的特点决定了即使是并行计算,在计算过程中仍可能存在重复读取的现象。例如,分块时把栅格a分入block(0,0)中运算,但是栅格a的相邻栅格b可能在分块时被分入block(0,1)中计算,这样a和b栅格就会在运算的过程中被重复地读取。为提高并行算法在gpu上的效率,如何分块尤为重要。为了防止在计算过程中对内存的反复读取,所以分块时,将需要参与计算的元素全部分入一个block中。例如需要block(0,1)计算k×k个栅格的流向,就以这k×k的矩阵为中心,向四周相邻的方向各多取一个元素,即取以这个矩阵为中心的(k+2)×(k+2)阶矩阵。这样在计算该k×k个栅格时将需要的所有元素完全加载到存取速度很快的共享内存中来,在计算的过程中,该部分栅格的运算就不需要再存取任何外部的内存。
　　针对绘制等流时线的运算,需要创建三个与运算栅格相同大小的辅助矩阵,即流向矩阵、坡度矩阵和单位栅格流经时间矩阵。如果想取得更高的运算速度,最好能将这三个矩阵都存放在存取速度较快的共享内存中。同时,在cuda平台上,一个线程块最大只支持512个线程,而且根据内存节距对齐的原则,gpu的内存控制器从某个固定倍数的地址开始读取才会有最高的效率[14]。根据其支持的线程数,最好的办法就是将高程矩阵分解为14×14的小矩阵,这样,每次共享内存中就会加载16×16阶的矩阵,并在共享内存创建三个14×14阶的辅助矩阵。为了方便进行线程块的计算,针对高程矩阵,可为每个线程块分配14×14个线程,再建立(m/14)×(n/14)个线程块。这就充分利用了gpu多线程的特点,从而大大地提高运算速度。当然,分解高程矩阵时,有可能小矩阵块是边沿矩阵块,这就需要在分解时对小矩阵块作出是否是边界元素的判断,如果是边界元素,就需要用-999作为值填充在空白的相邻栅格的位置,创建出16×16阶的矩阵。

由于参加运算的高程矩阵的阶很可能不是16的倍数,为了使gpu更有效率地工作,在开始运算前为分解的分块小矩阵开辟内存空间时,就直接把内存大小配置成16的倍数,并在复制矩阵到显卡内存之前将其清零。这充分利用了gpu的内存读取特点,且符合gpu内存读取高效率的原则。使用cuda提供的cudamallocpitch()函数就可以满足该要求,用于提高共享内存的访问速度[15]。
　　分配空间的伪代码如下:
　　#define block_size 16
　　……
　　int block_num = ((n + (block_size – 1)-1) / (block_size-1))*(block_size-1);
cudamallocpitch((void**)&demc, &pitch_dem, sizeof(float)*block_num, block_num);
　　cudamallocpitch((void**)&dirc &pitch_dir, sizeof(float)*block_num, block_num);
　　cudamallocpitch((void**)&gradc,&pitch_grad, sizeof(float)*block_num, block_num);
　　cudamallocpitch((void**)&timec, &pitch_time, sizeof(float)*block_num, block_num);
　　cudamemset(demc, 0, pitch_dem*block_num) ;
　　2)运算核函数
　　有了上述的高程矩阵分块方法和内存分配机制,真正需要在gpu上执行的核函数的伪代码如下:
　　__global__ static void gendlsx(const float*dem, size_t lddem,const float*dir,size_t ldir,float*grad, size_t ldgrad,float*time,size_t ldtime, int n)
　　{ 分配共享内存空间
　　dem[block_size][block_size];
　　将该线程需计算的分块矩阵的对应行列载入共享内存对应的位置;
　　__synthreads();
　　//线程同步,确保各个线程需要数据完全装载在共享内存中;
　　for(int i=1;i　　//只计算中间的14×14阶矩阵
　　{
　　计算出水流方向表,依据栅格水流方向,将其写入dir矩阵中;
　　根据水流方向,计算出水流方向的长度和坡度的比值,放入矩阵grad中;
　　依据式(1),计算出单个栅格的汇流时间,存入time矩阵;
　　}
　　__synthreads();//线程同步,确保以上各个线程计算均完成;
　　将结果写入dir和time矩阵对应的位置;
　　}
　　4 实验分析
　　实验平台为2.60 ghz的intel e5300双核处理器,2 gb内存,windows 7×64操作系统。gpu为nvidia geforce 9600gt,显存256 mb,使用2.2版本的cuda toolkit及对应sdk。编程环境为visual studio 2008。
　　
　　考虑到传统的计算方法难以满足高数据量的计算需求,本文算法的计算核心是矩阵相乘计算,需要花费大量的时间,gpu在计算矩阵相乘时,比cpu要快得多。因此,本文采用改进后的归并排序算法进行数据排序及新的内存分配策略和改进的并行算法进行汇流分析,可巧妙地利用基于cuda平台的gpu强大的并行处理能力来加速计算速度。
　　实验以四川彭州甘溪沟流域的dem图为例,对上述算法并行部分进行验证(参与运算的数据类型均定义为float类型)。彭州甘溪沟流域dem图采用16位grid数据集,划分为2 511×2 550个栅格,计算机中表达为二维矩阵形式。其中最大值为4 810.000 0,最小值为0.000 0,空白边界值为-9 999。实验结果如表1所示。由于数据量仅达到百万级的运算量,且其中包含加速效果不明显的排序运算,速度的提升不是很明显,仅达到4倍多。借助 gpu强大的并行计算能力,如果数据量进一步增多,将会大大提高运算速度。实验运算过程中的方阵相乘使用cuda运算的实验数据如表2所示。若矩阵阶数较小,由于gpu 运算必须经过pci-e总线传输及显存分配等必要耗时,加速比不明显;若矩阵阶数较大,如在本实验条件下,矩阵阶数达到450,加速比就开始增大。实验结果表明,基于cuda的并行算法比采用cpu上的串行算法具有更好的计算能力,能提高系统计算效率。
　　表1 生成基于dem的等流时线的运算时间比较
　　ms
　　
　　cpu上串行算法基于cuda的并行算法
　　3 207650
　　
　　表2 矩阵乘法运算时间比较ms
　　
　　矩阵阶数cpu上串行算法基于cuda并行算法加速比矩阵阶数cpu上串行算法基于cuda并行算法加速比
　　10<541无700657102x6.44
　　100<1043无10001880232x8.10
　　4508479x1.062000153081386x11.04
　　5 结束语
　　本文算法能提高系统的计算效率,具有较好的效果。实验表明, 对相同的数据,通过比较分析生成基于dem的等流时线的运算时间和矩阵乘法的运算时间,gpu比cpu的计算处理速度要快很多。基于cuda平台可在一台主机上运行多个gpu,通过cpu的线程来管理多个gpu,或建立 gpu集群, 使其计算性能得到更大的提升。可以将基于cuda的汇流分析的并行算法的设计思想应用于涉及到需求转换后能应用快速并行计算的各种项目系统中,以提高系统的计算效率。本文的算法设计思想和相关技术解决方法,为相关功能的实现提出了新的思路和解决办法,将对构建模型化、定量化、直观现代的新型自然灾害监测预警分析系统等提供重要的借鉴与参考。
　　参考文献:
　　[1]
　　吴恩华.图形处理器用于通用计算技术、现状及其挑战[j].软件学报,2004,15(10):1493-1540.
　　[2]张舒,褚艳利.gpu高性能运算之cuda[m].北京: