劳保强 安 涛 于 昂 郭绍光

(1 中国科学院上海天文台 上海 200030)

(2 中国科学院射电天文重点实验室 南京 210008)

(3 桂林电子科技大学信息与通信学院 桂林 541004)

1 引言

为了获得更高分辨率和灵敏度,综合孔径干涉阵列设备得到了飞速的发展.国际上,从上世纪70年代至90年代末,主要建成了厘米波、毫米波和亚毫米波的综合孔径干涉阵列.如美国的甚大阵(Very Large Array,VLA)[1]、澳大利亚望远镜致密阵(Australia Telescope Compact Array,ATCA)、英国多天线微波相连干涉仪网(Multi-Element Radio-Linked Interferometer Network,MERLIN)[2]、印度巨米波射电望远镜(Giant Metrewave Radio Telescope,GMRT)[3]等.进入21世纪,分米波、米波甚至更长波段的低频综合孔径干涉阵列重新成为了热点.如荷兰的低频阵列(Low Frequency Array,LOFAR)[4]、澳大利亚的默奇森宽视场阵列(Murchison Widefield Array,MWA)[5]、美国的长波长阵列(Long Wavelength Array,LWA)以及即将建造的平方公里阵列(Square Kilometre Array,SKA)[6]等.在国内,射电天文研究起步相对较晚,综合孔径干涉阵列屈指可数.目前的设备有新疆的21 cm阵(21CMA)、内蒙古的日像仪和新疆的天籁阵.从国际形式来看,射电天文的研究将进入以SKA为核心的新时代.

随着综合孔径干涉阵列设备的发展,原始数据的规模随之增大,同时数据处理也面临着巨大的挑战.SKA第1阶段(SKA1)的科学数据预处理所需要的计算能力为260 PFLOPS (即每秒260千万亿次浮点运算),相当于我国超级计算机天河二号(33 PFLOPS)的8倍,比目前世界排名第2的中国神威太湖之光的计算能力(93 PFLPOS)还高2倍,相当于约1亿台个人计算机的处理能力[7].近年,由于高性能计算和人工智能等技术的不断提高以及天文大数据处理的需求,使得利用这些技术来优化综合孔径干涉阵列数据处理方法成为了热门.

大视场成像是SKA等低频综合孔径干涉阵列数据处理中较为关键的技术之一.目前,比较有影响力并且将用于SKA的算法有:uv-faceting[8]、w-projection[9]、wstacking[10]、w-snapshots[11].其中,Romein[12]对w-projection算法进行了加速优化,研发了基于Compute Unified Device Architecture (CUDA)和Open Computing Language(OpenCL)的w-projection成像算法.该方法支持单个或多个GPU节点运行,并完成了相关性能测试,结果表明提出的加速方法比以往的方法在性能上提高了一个数量级.w-snapshots、w-projection和w-stacking成像算法均是一次性对整个视场的图像进行校正处理,这导致了优化后的算法只能有效地重建较低分辨率天空图像[13].相反,uv-faceting将视场分割成多个分面,每个分面采用传统的2维傅里叶变换方法重建天空图像,最后将所有分面的图像进行合并.在重建大像素天空图像上,有较好的拓展性,较适合用于SKA这种需要产生万兆像素图像项目中.此外,与其他方法相比,uv-faceting方法的另一个优点是可以在栅格化和成像之前直接对数据进行方向依赖效应的校正[14].因此,本文主要对uv-faceting算法进行并行优化改造.

首先,介绍了uv-faceting成像算法具体的工程实现方法,从而分析了该算法各个处理步骤的时间消耗情况,并确立对栅格化(gridding)和数据读取这两个部分进行并行优化.然后,提出了基于Message Passing Interface (MPI)+Open Multi-Processing(OpenMP)的uv-faceting成像算法和基于MPI+CUDA的uv-faceting成像算法.最后,对这两种算法进行正确性的验证与分析和性能测试与分析.

2 uv-faceting成像

本节主要从工程实现的角度描述uv-faceting成像算法,我们之前的研究给出了详细的uv-faceting算法理论推导过程,最终给出了第j个分面的亮度分布函数I与其可见度函数Vfacet的关系式[15−16]:

上式中,第j个分面的可见度函数可以通过改变原始可见度函数V的相位得到:

其中,lj、mj和nj是第j个分面中心的方向余弦,e表示指数函数,∆lj=l−lj和∆mj=m−mj,l和m是原始观测方向的方向余弦,u′和v′由原始的u、v和w坐标计算得到:

这里,图像坐标l和m在工程实现中可以利用其与相位中心的关系得到[17],因此lj、mj和nj可以通过下面的公式计算得到:

其中,nra和ndec是分面相位中心的赤经和赤纬,odec和ora是原始相位中心的赤经和赤纬.

从上述的分析,我们很容易得出理论上只要给出需要求解的分面的相位中心,就可以利用传统的2维傅里叶变换方法求出该分面的亮度分布函数.因为原始的可见度数据、相位中心和u、v和w坐标可以在观测数据中获取.

图1是实际观测数据进行uv-faceting成像的流程图.这里我们使用的输入数据必须为校准后的数据(calibrated data),格式是Measurement Set (即ms文件)1http://www.astron.nl/casacore/trunk/casacore/doc/notes/229.html,它是目前比较通用的数据存储格式.该文件格式目前版本号是V2.0,SKA科学数据处理器团队正在将其升级为V3.0,并将作为SKA1数据的文件存储格式.图1中uv-faceting成像主要由3层嵌套循环来完成,最外层循环是按照分面的个数执行,使用facet-id变量作分面的循环计数,num-facets是分面的总数.次外层循环按照可见度数据的行数来执行,这里使用变量row作行数的循环计数,可见度数据的总行数为row-count.内层循环按照频率通道数来执行,变量c作为频率通道数的循环计数,总的频率通道数为channel-count.1次循环即完成1个分面1行可见数据中1个频率通道的栅格化(gridding)处理任务.gridding的实质是可见度函数与卷积核函数(卷积核滤波器)进行卷积,这里的gridding算法我们将采用casacore grider的设计思路实现.循环结束后,对每个分面得到的gridding结果进行2维快速傅里叶逆变换(IFFT),得到各个分面的图像结果.最后,马赛克(mosaic)处理所有分面图像得到最终的目标图像.其中,步骤phase rotation利用(2)式计算完成;步骤baseline rotation由(3)式和(4)式联立计算得到;步骤scaleuvcoverage根据傅里叶变换的展缩(尺度变换)特性获得.例如当想获得的图像像素大小为nx×ny,并且像元大小(cellsize)设定为xcellsize和ycellsize时,u、v坐标需要分别乘以一个尺度常量,分别为:

图1 uv-faceting成像流程图.图中N表示否,Y表示是Fig.1 The flow chart of uv-faceting imaging.N indicates no,and Y indicates yes

3 uv-faceting成像并行算法

3.1 基于MPI+OpenMP的uv-faceting成像并行算法

首先,我们探讨一下uv-faceting成像算法处理过程的耗时情况.从传统的射电成像中得知,在成像的所有处理步骤中gridding是耗时最多且是较为关键的一步.另外,随着输入数据的增大,数据的读取所需的时间也会急剧增大.因为ms文件存储底层的数据输入输出(IO)使用的是casacore库中的Casacore Table Data System (CTDS)[18],目前CTDS只支持数据的串行读写.所以当数据量非常大的时候,数据读取消耗的时间将是巨大的,这将严重影响后续的科学数据处理与分析任务的开展.如图2是uv-faceting成像关键步骤随着数据量增大所耗时间的变化曲线,关键步骤依次为卷积核滤波器的建立(Conv-filter Create)、数据加载(Data Loading)、栅格化(gridding)和傅里叶逆变换(Fourier Inversion).数据大小分别是0.36、3.7、36和349,单位GB.图中也充分证明了数据读取和gridding是uv-faceting成像算法最耗时间的部分,其中gridding耗时最长,并且这个增大呈线性性.基于此,本文将重点对gridding和数据读取这两个部分进行并行化改造.

图2 uv-faceting成像关键步骤耗时曲线.Conv-filter表示卷积核滤波器Fig.2 The consuming time curve of the key step of uv-faceting imaging.Convfilter represents the convolution kernel filter

在gridding的并行改造方面,从图1流程可以看出,uv-faceting算法实际上提供了一种粗粒度并行的方法.这种方法只需对最外层循环进行并行化,无需考虑任何内层循环的处理过程和同步锁问题,最为简单方便.这个粗粒度并行可以使用OpenMP编程实现,本文则是在最外层循环外加“#pragma omp parallel for schedule (static)”语句实现多线程并行.这相当于将分面的处理过程分为多个线程并行实现.例如有n个分面需要处理,并使用t个CPU核数或线程处理时,那么将给每个线程静态分配n/t个循环.当使用的CPU核数或线程数接近于分面数,该方法所消耗的时间最短.

在数据读取方面,大致有两种方法作并行优化:底层存储数据IO并行化和应用层数据读取并行化.第1种并行化方法是利用并行IO软件包,如Adaptable IO System(ADIOS)等,开发基于该并行IO的数据管理系统.它可以作为新的管理系统与传统的串行数据管理系统并存,同时将复杂的并行化代码封装起来并提供一种简洁的接口给用户使用,因此使用时避开了繁琐的编程实现过程.这种方法能够从根本上解决海量数据读写上的瓶颈,但是目前仍处于验证性实验阶段,还没得到广泛应用[19].第2种并行优化方法是保留存储底层数据IO的读写方式,在用户应用层通过编写并行程序对数据进行分割同时读取.这种方式需要对数据格式和内容有充分的了解,且具备较强的编程功底.本节的并行算法将使用第2种方法.

图3是本文提出的基于MPI+OpenMP的uv-faceting成像并行算法.首先,根据数据的存储格式特点将输入的校准后的数据利用MPI进程按行分割成多个分块,分割数等于MPI进程数,图3使用4个进程或4个节点为例;其次,使用MPI将分割之后的数据分配到不同的CPU节点上进行并行处理,每个进程或CPU节点将完成gridding等数据处理.同时在每个进程使用OpenMP多线程并行完成多个分面的gridding处理;然后,将所有进程的gridding结果进行求和并发送至主节点;最后,在主节点对最终的gridding结果进行2维傅里叶逆变换得到目标多个分面的图像并输出相应的fits格式文件,再将多个分面的图像进行马赛克处理得到最终的目标图像.本文中马赛克的处理软件采用Swarp[20].

图3 基于MPI+OpenMP的uv-faceting成像并行策略.k表示观测时间戳的开始索引,m表示观测时间编号,pol表示偏振Fig.3 The uv-faceting imaging parallel strategy based on the MPI and OpenMP.k is the start index of observation time stamp,m is the numbering of observation time,and pol is the polarization

3.2 基于MPI+CUDA的uv-faceting成像并行算法

就单节点而言,CPU的核数是非常有限的,而GPU核处理器则通常有成千上万个处理核.因此,上述的gridding处理步骤,使用GPU加速会更加有效.本文主要采用CUDA编程实现gridding加速,利用了由相同基线进行连续测量之间的空间局部性.由于连续测量之间的积分时间足够短,相同基线的连续测量数据将会被gridding到相同的网格点上.我们把相同基线下连续测量的多个时间戳数据的gridding结果进行累加并缓存在寄存器中,当处理下一个基线数据时,再将寄存器的数据写入全局存储器(global memory).当global memory的写访问是原子性(atomic)的,多线程之间的不存在竞争,多个基线就可以实现并行gridding.在gridding处理中,我们给每个卷积核滤波器的抽头分配一个线程,一共需要CFull support-size×CFull-support-size个线程,其中CFull-support-size是卷积核滤波器的全支撑大小.因此,基于该方法的uv-faceting成像算法需要的最小线程数量为Nfacets×NBaselines×CFull-support size×CFull-support-size,其中Nfacets是分面总数,NBaselines是基线总数.

图4是基于CUDA的单个GPU节点uv-faceting成像算法流程.图中主机先进行输入参数的初始化,同时对GPU device进行初始化,包括device的重置和device的选择,因为有些GPU节点会挂载多个GPU加速卡.随后,主机开始读数据,再将读入的数据进行重新排序和打包,最后将重排后的数据传送至GPU的global memory.GPU device经过初始化后,将会事先给从主机传送来的数据开辟相应的数组memory.数据传送到GPU global memory后,在GPU device中完成数据的gridding处理任务.首先,获取CUDA线程索引(thread-id),这里线程总数设置为上述提到的最小线程数量,线程块大小为256;其次,根据这个线程索引、设置的线程总数和线程块数量,计算出当前的基线起始索引、卷积核滤波器支撑索引(u-support和v-support)、每个基线的时间戳数量(num-t)和分面ID号(facet-id);然后,开始完成当前基线和分面下,多个时间戳的gridding并累加,累加的结果缓存到CUDA寄存器中,其中,这个gridding累加只计算当前卷积核支撑索引下的卷积结果.当需要进行下一个基线数据处理前,当前基线gridding结果将从寄存器拷贝到global memory;最后,把所有基线和分面计算的gridding结果传送给主机.主机完成IFFT和mosaic,得到最终的图像.

在数据读取方面,我们同样采用3.1节中的第2种方法,不同的是分割数据时按照观测时间的顺序来进行切割.如图5所示,是基于MPI与CUDA的uv-faceting成像并行策略.首先,将输入数据利用MPI进程按观测时间分割成多个ms文件,文件数等于MPI进程数,图5中同样是以4个GPU节点为例;其次,使用MPI将分割之后的文件数据分配到不同的GPU节点上进行并行处理,每个GPU节点将完成数据的重排序和gridding等数据处理,其中,gridding采用上述的CUDA编程进行加速处理;然后,将所有GPU节点的gridding结果进行求和并发送至主节点(GPU node0);最后,在主节点完成分面成像与马赛克处理获得最终的图像.

从上述的描述知,使用本文提出的并行算法的前提是观测数据以基线为一组并按照时间排序方式存储,即数据按照基线索引来排序的.然而,ms格式的数据一般不是按照这种方式存储的,因此在数据处理之前需要进行数据重新打包和排序.此外,需要注意的是:所有基线的观测时步的数量(时间戳的个数)并不是都相等,有些基线的某些时步数据有可能在标记和校准的时候就从Measurement Set中分离出来,所以需要一个可变长度的数组来存储每个基线的数据.如此,这个可变长度的数组想要传输到GPU上进行计算就非常困难.本文则使用一个一维数组存储每个基线第1个时步位置(即基线的起始索引)来替代.然而,Measurement Set没有存储基线索引的信息,只存储了天线的ID号.因此本文给出了一种计算基线索引的方法.首先,按照基线对数据进行分组,每组对应一个编号.这个编号可以使用天线的ID号来计算得到:

其中,B表示天线1 (ant1)和天线2 (ant2)组成的基线;S表示ant1和ant2中最小ID号;Nant表示天线的个数;Iant1和Iant2分别表示ant1和ant2的ID号.根据(6)式计算出基线的编号,然后统计相同基线编号下时步的个数(相同基线的采样数据个数),并存储到数组C中.最后,利用前缀求和方法计算出每个基线的起始索引Bstart-index,其数学表达式如下:

其中,ks表示相同基线下时步的编号,范围为0 ∼IB−1.

图4 基于CUDA的单GPU节点uv-faceting成像算法流程图Fig.4 The flow chart of the uv-faceting imaging algorithm based on the CUDA on a single GPU node

图5 基于MPI+CUDA的uv-faceting成像并行策略.k表示观测时间戳的开始索引,m表示观测时间编号,pol表示偏振Fig.5 The uv-faceting imaging parallel strategy based on the MPI and CUDA.k is the start index of observation time stamp,m is the numbering of observation time,and pol is the polarization

结合上述计算得到的基线索引起始值和相同基线索引下的时步个数,我们就可以将读入的数据按照基线索引进行重新排序.

4 验证性实验与分析

本次实验的数据采用甚大阵VLA B和C构型组成的望远镜阵列在74 MHz波段观测后发座星系团(Coma cluster)得到的数据,观测时间为:1998-11-22UTC11:19:57.5—1998-12-02UTC22:59:55.0.该观测的时间总长为905998 s,分为3组观测.

Coma cluster是一个超大的星系团,它包含超过1000个已经被识别的星系,视场范围较大,因此非常适用于本文算法的验证测试.

实验的测试平台采用广州超算中心的天河二号超级计算机,由于测试需要用到GPU,因此实验最终在GPU节点上进行.天河二号目前共有25个GPU节点可用,每个GPU节点各指标参数如表1所示.基于MPI+CUDA的算法测试采用10个节点进行测试,每个节点使用1个GPU加速卡.基于MPI+OpenMP的算法测试同样采用10个节点进行测试,但每个节点只使用CPU核,每个节点共使用20个CPU核(即20个线程).

图6是Coma cluster数据进行传统的2维傅里叶变换成像得到的结果.图7是使用通用天文应用软件(Common Astronomy Software Applications,CASA)进行faceting成像得到的结果,其中图像的像素大小设为1800 × 1800,赤经和赤纬方向的图像单元大小(cellsize)均设为30′′(arcsecond),分面总数为25个分面(赤经和赤纬方向各为5个分面).图8是本文基于MPI+CUDA的uv-faceting成像得到的结果,图9是本文基于MPI+OpenMP的uv-faceting成像得到的结果,这两个测试的cellsize大小和分面数量的设置均与CASA软件的设置相同.从视觉效果可以明显看出,远离相位中心的星系或射电源在图7、图8和图9均能够得到较好重构.这说明本文提出的这两种并行算法能够完成射电干涉阵列大视场成像,结果与广泛使用的CASA软件的结果基本一致.此外,图8和图9在相同的流量范围内,重构的图像更加明亮.说明我们提出的方法,在提高了运行速度的同时也提高了图像的信噪比.

表1 天河二号超级计算机GPU节点各指标参数Table 1 The indicator parameters of the GPU node of the Tianhe-2 supercomputer

图6 2维傅里叶变换成像结果Fig.6 The results of two-dimensional fourier transform imaging

图7 CASA软件的成像结果Fig.7 The imaging results of CASA software

图8 基于MPI+CUDA的uv-faceting成像结果Fig.8 The results of uv-faceting imaging based on the MPI and CUDA

图9 基于MPI+OpenMP的uv-faceting成像结果Fig.9 The results of uv-faceting imaging based on the MPI and OpenMP

为了更进一步验证我们提出的算法的正确性和有效性,我们使用源查找软件对图7、图8和图9进行源查找,对比分析它们的图像结果的一致性.这里我们使用的源查找软件是由Paul Hancolk开发的Aegean软件[21],该软件已经广泛应用于默奇森广角阵列MWA的数据处理当中.源查找的对比结果如图10、图11和图12所示,这里我们采用的波束大小(beam size)均为193.20′′× 142.96′′,位置角为−64.46◦,这些值由原始观测数据计算得到.表2为详细的分析对比结果.从表中可知,以CASA得到的结果为基准,MPI+CUDA方法在1个波束大小的误差范围内与CASA得到的结果位置相同的射电源有140颗,相同率约为94.6%,比CASA的结果多恢复了11颗源.MPI+OpenMP的方法得到的结果有141颗射电源与CASA的结果位置相同,相同率为95.27%,比CASA多恢复了42颗射电源.这表明本文提出的这两种方法基本正确.存在差异的原因:我们的方法是将数据分割成多份,每份数据单独进行栅格化,再进行加权求和得到最终的栅格化结果.在加权求和时,边界位置部分网格点出现了重复累加,导致差异的产生.这些差异可以通过以下方法进行修正:以CASA图像结果的流量大小为基准,对我们提出方法的图像结果进行加权重,修正它们的流量大小,从而获得与CASA完全一致的图像结果.这项工作是我们下一步的工作内容.另一方面,MPI+OpenMP与MPI+CUDA的结果有159颗射电源的位置相同,说明MPI+OpenMP的方法多恢复了31颗源.表明在正确性上,MPI+CUDA的方法比MPI+OpenMP的方法更加准确.

在运行速度方面,MPI+CUDA的方法消耗的总运行时间最短,CASA消耗的总运行时间最长,相比CASA、MPI+CUDA的方法提高了将近13倍.虽然MPI+OpenMP与MPI+CUDA消耗的总运行时间差不多,它们只相差不到3 s.但在gridding方面,MPI+CUDA的方法消耗的时间为1.16 s,MPI+OpenMP消耗的时间为3.92 s,说明MPI+OpenMP是MPI+CUDA的3.38倍左右.综上分析,表明本文提出的算法在运行速度上有很大的提高,不管是在运行速度还是正确性上MPI+CUDA方法更优,但也有上升的空间.

图10 图7与图8源搜寻结果对比图Fig.10 Comparison of source finding results of Fig.7 and Fig.8

图11 图7与图9源搜寻结果对比图Fig.11 Comparison of source finding results of Fig.7 and Fig.9

图12 图8与图9源搜寻结果对比图Fig.12 Comparison of source finding results of Fig.8 and Fig.9

表2 源查找结果对比Table 2 Comparison of source finding results

5 性能测试与分析

5.1 单节点

单节点的测试平台,我们采用上海天文台区域中心原理样机的GPU节点、澳大利亚Pawsey超算中心Galaxy超级计算机的GPU节点和上述天河二号的GPU节点进行测试.上海天文台GPU节点的各项指标参数见表3,Galaxy超级计算机的GPU节点各项指标参数见表4.测试数据采用第4节中Coma cluster的观测数据,分面总数设置为49个分面,即赤经和赤纬方向各7个分面,其余的参数设置与第4节中的测试相同.我们主要通过改变环境变量OMP-NUM-THREADS来使用不同的线程数量测试基于MPI+OpenMP的uv-faceting成像算法的性能.然后,通过使用不同型号的GPU来测试基于MPI+CUDA的uv-faceting成像算法的性能.

如图13所示,在56个线程范围内,随着线程数量的增加,gridding所消耗的时间随之递减.但没有呈现良好的线性性,不管是单精度还是双精度计算.极有可能是分面的数量没有均匀分配给CPU核,导致了各个线程之间的负载不均衡.当线程的数量接近分面总数时,gridding所消耗的时间最少;不管是单精度还是双精度情况下,相比没有进行并行改造的运行速度均增加了20倍左右.表明本文提出的基于MPI+OpenMP的uvfaceting成像算法的并行优化是有效的.另外,我们还对比测试了CPU和GPU之间以及不同GPU之间的性能.随着GPU加速卡性能的提升,gridding消耗的时间变短.如图13特斯拉(Tesla)V100所消耗的时间最短,无论是单精度还是双精度相比没有进行并行改造的运行速度增加了60倍左右.单精度情况下,Tesla K80 gridding所消耗的时间与CPU使用48个线程所消耗的时间相当,Tesla K20X gridding所消耗的时间与CPU使用40个线程的相当.双精度情况下,Tesla K80 gridding所消耗的时间与CPU使用8个线程所消耗的时间相当,Tesla K20X gridding所消耗的时间与CPU使用4个线程的相当.表明本文提出的基于MPI+CUDA的uv-faceting成像算法是有效的.同时也表明,在一定的线程范围内,基于MPI+CUDA的算法比基于MPI+OpenMP的算法性能要好,并且随着GPU加速卡的提高,这个范围的上限值会随之增大.

表3 上海天文台区域中心原理样机GPU节点各指标参数Table 3 The indicator parameters of the GPU node of the Shanghai Astronomical Observatory Regional Center Prototype

表4 Galaxy超级计算机GPU节点各指标参数Table 4 The indicator parameters of the GPU node of the Galaxy supercomputer

5.2 多节点

多节点的测试数据我们采用扩展甚大阵望远镜(Expanded Very Large Array,EVLA)的观测数据,选择了Sanjay Bhatnagar的超新星遗迹G55.7+3.4 L波段的观测数据.该观测于2010-08-23UTC01:00:25开始,2010-08-23UTC08:25:16.5结束,历时约7小时25分.首先我们使用CASA软件对G55.7+3.4源的原始观测数据进行校准,并导出校准后的数据作为本次测试的输入数据.

本次测试的平台采用天河二号的GPU节点,如表1所示.测试数据的参数设置:像素大小设为1024×1024,赤经和赤纬方向的图像单元大小(cellsize)均设为8′′,分面总数为16个分面(赤经和赤纬方向各为4个分面).测试方法:计算节点从1到10自然数递增,每个节点数测试10次,获取数据加载的时间和gridding的时间,并取平均值和标准差.基于MPI+CUDA的方法每个节点使用1个GPU加速卡,基于MPI+OpenMP的方法每个节点使用20个CPU核.

图13 单节点性能测试结果Fig.13 The results of performance tests on a single node

数据读取(或加载)的耗时情况如图14所示.在GPU上运行时,在1–4个节点范围内,无论是单精度情况下还是双精度情况下,随着节点数的增加所消耗的时间均接近线性递减.当节点数超过4个之后,随着节点数的增加,所消耗的时间减少得比较缓慢,在8个节点数之后达到平衡.在CPU上运行时,只有在2个节点范围内接近线性递减.当节点数超过2个之后,变化缓慢,同样在8个节点数之后达到平衡.表明本文提出的基于MPI+CUDA方法和基于MPI+OpenMP方法在一定计算节点范围内具有良好的可拓展性,基于MPI+CUDA方法的范围更大,拓展性更优.另一方面,在CPU上运行时,所消耗的数据加载时间相对较少,单精度和双精度所消耗的时间基本相同.此外,CPU上运行的标准差较GPU上运行的小,运行相对稳定.在GPU上运行,所消耗的时间较大一些,且不够稳定.表明在数据加载方面,基于MPI+OpenMP的方法比基于MPI+CUDA的方法稳定性强且耗时短.

gridding的耗时情况如图15所示.从图中可以看出,无论是CPU还是GPU上运行,gridding所消耗的时间均随着节点数的增加逐渐减少,但是均没有呈现明显的线性性,在8个节点数之后达到平衡.1个节点数所消耗的时间是10个节点数的9倍左右.表明本文提出的两种方法在整体上具有良好的拓展性.在单精度情况下,在GPU上运行所消耗的时间较短,标准差几乎为零,相对较稳定.在双精度情况下,在CPU上的运行时间较短,但GPU的运行稳定性更强.表明在gridding方面,基于MPI+CUDA的方法比基于MPI+OpenMP的方法更加稳定.

图14 数据加载耗时曲线Fig.14 Data loading consuming time curve

图15 gridding耗时曲线Fig.15 Gridding consuming time curve

6 总结与展望

本文主要介绍了uv-faceting的工程实现流程,进而分析了uv-faceting主要步骤的耗时情况,明确了数据读取和gridding这两个步骤是uv-faceting中最耗时的部分.基于此,本文提出了两种并行优化算法:基于MPI+OpenMP的uv-faceting成像算法和基于MPI+CUDA的uv-faceting成像算法.验证性实验的结果表明本文提出的两种算法均能够正确实现uv-faceting成像,结果与目前广泛使用的CASA软件得到的结果基本一致.但还存在略微差异,这是我们将来工作需要完善的地方之一.性能测试方面,单节点的测试表明本文提出的两种算法是有效的,在一定线程范围内,基于MPI+CUDA的方法比基于MPI+OpenMP的方法性能更好.其中,反映了基于MPI+OpenMP的方法线程之间负载均衡不是很好,多线程的拓展性有待优化和提高.多节点的测试表明,在数据读取方面,本文的两种方法在一定的节点数范围内具有较好的可扩展性,基于MPI+OpenMP的方法比基于MPI+CUDA的方法更加稳定.在gridding方面,表明本文提出的两种方法在整体上有良好的可拓展性,但在局部上均没有展现出明显的线性性,说明节点之间的负载均衡不是很好.因此,在后续的工作中,需要完善多节点之间的可拓展性性能.