苏 辉a*,邱夏青b,马文鹏b

(信阳师范学院 a.网络信息与计算中心;b. 计算机与信息技术学院,河南 信阳 )

摘 要:基于Matlab平台,采用有限元方法实现了对二维拉普拉斯(Laplace)方程在GPU平台上的加速.通过对物理问题的分析与物理模型的构建,完成总体CSR格式存储的刚度矩阵的生成;使用Matlab和CUDA混合编程,在Matlab平台上实现该有限元问题的并行加速;并结合CuBlas数值计算库采用PCG算法求解装配后的大型线性稀疏方程组,从而高效地迭代出各格点的速度势.该算法既充分发挥了Matlab在数值计算方面的高效性,又充分发挥了GPU在细粒度并行加速方面的优势.

关键词:CUDA程序设计;有限元方法;GPU;预处理共轭梯度算法

对于很多大规模数值模拟来说,稀疏线性系统的求解占用了大部分的计算时间和资源.随着科学与工程问题的规模和复杂程度的提高,对稀疏性系统的求解规模和速度提出了更高的要求,利用高性能计算技术是当前解决大规模稀疏线性系统求解的重要途径[1].在图形处理器(GPU)硬件架构的迅猛发展和NVIDIA CUDA[2]并行计算平台的不断完善之下,GPU的可编程能力也在不断增强.随着CUDA运算平台的推出,GPU的并行计算已经渗透到各个学科,如在生物科学研究领域中蛋白质与基因序列串的匹配等;在计算金融方面利用GPU加速债券定价、风险分析以及算法交易[3];在化学计算方面进行量子化学的计算加速;在数据挖掘、分析学方面的信息高效排序等,其本质就是利用GPU多线程的特点来辅助CPU计算,以此来提高计算机整体的运算速度.

目前,针对GPU的有限元分析在求解大规模稀疏线性方程组的研究较多,主要是由于在求解线性方程组时数据量大并且计算比较集中,因此主要研究稀疏矩阵的数据结构和稀疏矩阵向量乘(Sparse Matrix Vector Multiplication, SpMV)的算法[4].同时,很多研究学者也研究了针对GPU的基于CUDA平台的共轭梯度(Conjugate Gradient Method, CG)算法及各种改进算法,如BICG、BiCGSTAB算法.

本文所述算法是采用行压缩存储(Compressed Row Storage, CSR)模式存储大规模的稀疏矩阵,在稀疏矩阵向量乘(SpMV)的基础上采用预处理共轭梯度(Preconditioned Conjugate Gradient Method, PCG)算法,对大规模稀疏矩阵组成的线性方程组迭代求解.采用CUDA和Matlab混合编程,实现有限元方程的数值求解,通过多线程求解来实现GPU加速,算法性能有明显的优越性.

有限元方法(Finite Element Analysis, FEA)是一种高效常用的数值计算方法,由于有限元方法计算精度较高,具有适应大多数复杂情况的能力,使其成为力学工程分析问题的有效解决手段.化整为零、积整为零是有限元方法的基本思想,这种解决问题的思想与并行计算中的分治思想相协调[5].因此,对于规模较大的物理工程问题,可以借助于有限元方法将大规模问题分解,使用GPU多线程结构将运算效率低下的串行方法改进成为并行分析方法,使得计算机的执行效率变得更高.

矩阵的类型有很多,按照存储数据的稠密程度分类,可以分为稠密矩阵和稀疏矩阵.在稀疏矩阵中,其存储的零元素较多,而其存储的有效值较少;而稠密矩阵所存储的有效值较多,非有效值较少.因此本文采用一种新型的存储方式,即将矩阵信息压缩存储,也是将矩阵中的零元素都剔除,只存储矩阵中的有效值,对于内存空间的利用将会更加充分.

SpMV算法是用来求解矩阵和向量相乘的函数,通过SpMV算法,可以求得矩阵与向量相乘的结果,此算法主要是为后面的预处理共轭梯度(PCG)算法做铺垫、以供给其调用来求解线性方程组.

据科学统计,预条件共轭梯度(PCG)算法的迭代速度随着方程的阶数的增加而下降,从而可以推断出:对于高阶稀疏矩阵的线性方程组,无论其阶数有多高,PCG算法迭代的次数都是可以被接受的[6],这种方法使求解大型或者超大型方程组成为可能.同时,对于PCG算法,还允许在计算的过程中保留矩阵的稀疏性,允许矩阵采取更加高效的方式来存储它.此特点使得计算机在处理拓扑形状复杂的稀疏矩阵的时候节省大量的内存空间.同时,PCG算法适合于求解大型对称的矩阵,而本次实验所测试的矩阵正好符合这一特点.

Matlab是一款通常用于数据可视化、数据计算、算法开发以及数据分析数学软件[4].同时,Matlab可以将数值计算、科学数据可视化、矩阵运算以及图形图像仿真都集成在一个视图中,使得工程设计和科学研究中的有效数值计算问题能够得到更加有效的解决.Matlab最大的优点就是能够使数值计算结果和程序实现过程可视化,其强大的图形图像处理能力使得静态的数据转化为生动的图像,使得数学规律更加直观明了,使人从枯燥无味的数学计算中解脱出来.其次,在Matlab平台可以调用CUDA函数,实现复杂计算调用GPU中的核函数完成,从而实现整体运算加速[7].

理想流体流动的连续方程为:

其中:表示x方向的速度,其大小等于速度势φ在x方向的导数,表示y方向的速度,其大小等于速度势φ在y方向的导数,将其代入式(1)中得:

对其可描述为理想流体流动的速度势方程.该方程的形式与拉普拉斯方程的形式一致.

三角形线性单元插值函数为:

Φi=ai+bix+ciy,i=1,2,3,

在计算公式(3)和(4)时,当i=1时,j=2,k=3;当i=2时,j=3,k=1;当i=3时,j=1,k=2.其中Δ为单元面积,可计算如下:

在伽辽金有限元方法中,权函数等于插值函数.式(2)与插值函数相乘后再整个计算区域内积分,可以得到:

式(6)的左边用格林公式进一步推算和展开,得到系统(2)的加权余量方程:

单元方程的形式如下:

但是,由于存在第二类边界条件,则在建立单元方程时,应该分为内部单元和边界单元,这里不再叙述.

总体方程的组合分两步完成,首先装配CSR格式的刚度矩阵;最终装配目标刚度矩阵和装配好边界向量.

在装配完刚度矩阵和求解出边界向量之后,就可以使用预处理共轭梯度(PCG)求解线性方程组了.在求解线性方程组时,需要先求解出CSR存储格式的系数矩阵和由边界条件计算的向量,然后再调用PCG算法实现线性方程组的求解.由此计算出来的时间可以进行与后面的并行效率对比分析.

网格是GPU线程的组织方式,若干线程块组成了一个网格,最多512个线程能组成一个线程块,属于同一个线程块中线程的指令地址相同,这些线程具有很强的可并行性,同时还能够通过Share memory和barrier实现块内通信,可以形成细粒度并行[8].在线程中通过粗细粒度的双层并行可以实现GPU运算的多方面加速.

由于GPU多线程的特点,使其在大规模并行计算方面具有很大的优势,所以对于此次设计中的装配刚度矩阵和预处理共轭梯度算法来说,采用在GPU端细粒度并行方式可以实现计算提速[9].

在本次刚度矩阵的填充中,填充的刚度矩阵格式为CSR模式,直接将每行的数值填入value数组,填充位置的确定需要根据index数组确定每行非零元素在value数组中的起始地址和每行的非零元素的序号来确定.最终的value装配形式是按照每行的数据段从前到后顺序执行.

3.2.1 并行装配刚度矩阵

在Matlab中调用MEX文件极为方便,因为其调用方式和Matlab的内建函数几乎完全相同,只需要在命令窗口输入函数名和对应的参数即可.因此,本文在Matlab平台上编写一个接口函数mexFunction,其格式如下:

void mexFunction(int nlhs,mxArray*plhs[],int nrhs,const mxArray*prhs[])

nlhs:输出参数的数目

plhs:指向输出参数的指针

nrhs:输入参数的数目

prhs:指向输入参数的指针

mexFunction函数只是提供了一个Matlab和CUDA的接口,其实现的功能就是获取参数指针,并且调用FillValueArray函数来启动核函数,实现多线程的计算[9].

在FillValueArray.cu文件中,面临的主要任务就是用C语言实现对value数组的按格点填充.对于每个格点分配一个线程,使其并行完成对value数组的填充.在GPU的每个thread中,执行的操作都相同,区别是每个线程执行时的数据不同.

实现CUDA和Matlab的混合编程,把.cu文件用CUDA的编译器nvcc编译成.obj文件,作为一个库函数供.cpp调用,然后.cpp与.obj文件共同用MEX接口编译成供Matlab调用的函数.在此部分最核心的部分就是用CUDA语言编写.cu文件,实现对CSR存储格式的系数矩阵的在GPU上并行填充.完成在device端分配显存、数值计算、写回host端和显存释放等操作.

3.2.2 并行求解线性方程组

在核函数中,首先,为要使用的变量开辟临时显存空间,声明调cusparseDcsrmv函数中的各个参数类型;其次,使用CUBLAS库中的一些固有的函数进行向量的操作;最后,分配线程执行函数,将device端计算的结果写回host端,释放在显存中开辟的临时空间.

对于CUDA中的可执行文件,主要是获取PCG函数的参数,获取变量的指针,在显存中开辟显存空间,然后再调用CUDA编写的.cu函数,完成并行操作.

对于PCG算法的并行,也是在Matlab平台用CUDA语言编写.cu文件与.h文件用nvcc命令编写成.obj文件,执行在GPU上[10].然后.cpp与.obj文件共同用mex接口编译成供Matlab调用的函数.在Matlab端对生成的函数调用,实现在Matlab平台实现对大规模线性方程组的并行求解操作.

本次实验是在Windows7,64位操作系统上运行,处理器的型号为Intel(R)Core(TM)i7-4790,8核,主频为3.60 GHz,安装内存(RAM)为8.00 GB,GPU型号为NVIDIA GeForce GT 705.本次实验运行的软件运行平台是Matlab 2015a、Visual Stdio 2010和CUDA8.0版本.

本次实验的测试算例是以流水通过两平板之间的速度势变化为物理背景.测试的初始条件为宽H为5 m,高L为10 m,水流速度v为2 m/s,其中水平方向分割段数Nx为500,垂直方向分割段数Ny为400,则通过计算可以得出划分的网格数为40万,网格格点数目为20万.

通过计算大规模线性方程组为载体,来对比通过GPU加速前后程序的执行时间,计算GPU在大规模工程运算中的提速效率.下面通过测试各个程序段的时间来对比使用GPU加速的效果.

1)实验的并行提速效果

装配刚度矩阵(fillvalue)提速:3.222÷0.401=8.035倍;

求解线性方程组(PCG)提速:2.810÷0.411=6.84倍;

算法总体提速:6.957÷1.737=4.01倍.

2)并行加速效率分析

通过上述数据,可以看出在装配刚度矩阵中使用GPU中的线程计算可以提速8倍,在求解线性方程组时使用GPU中线程可以提速7倍,对于求解整体的方程提速为4倍,对于刚度矩阵的填充和线性方程组的求解来说,效果比较理想.

本文是在Matlab平台上实现Matlab和CUDA的混合编程,借助于GPU的多线程结构,使得GPU发挥细粒度方面的并行优势,实现整体算法的并行加速.综上所述,本次设计中通过在GPU中的两次提速,在大规模数值计算问题中,能节省较多的计算时间,实验效果比较理想,达到在大规模数值计算的GPU提速.

参考文献:

[1] 张健飞,沈德飞.基于GPU的稀疏线性系统的预条件共轭梯度法[J].计算机应用,2013,33(3):825-829.

ZHANG Jianfei, SHEN Defei.GPU-based preconditioned conjugate gradient method for solving sparse linear systems[J]. Journal of Computer Applications,2013,33(3):825-829.

[2] 张健飞,沈德飞.有限元GPU加速计算的实现方法[J].计算机辅助工程,2014,23(2):41-45.

ZHANG Jianfei, SHEN Defei. The realization method of finite element GPU acceleration calculation[J].Computer Aided Engineering,2014,23(2):41-45.

[3] 胡国耀.基于GPU的有限元方法研究[D].武汉:华中科技大学,2011.

HU Guoyao.Research on finite element method based on GPU[D].Wuhan: Huazhong University Science and Technology, 2011.

[4] 原思聪.Matlab语言与应用技术[M].北京:国防工业出版社,2011.

YUAN Sichong.Matlab language and application technology[M]. Beijing: National Defence Industry Press, 2011.

[5] SUH J W, KIM Y. 2-Configurations for Matlab and CUDA[J]. Accelerating Matlab with GPU Computing,2014:19-44.

[6] 韩琪,蔡勇.基于GPU的大规模拓扑优化问题并行计算方法[J].计算机仿真,2015,32(4):221-226.

HAN Qi,CAI Yong.Parallel computing method for large-scale topology optimization problem based on GPU[J]. Computer Simulation,2015,32(4):221-226.

[7] MARTNEZ-FRUTOS J, MARTNEZ-CASTEJN P J,HERRERO-PREZ D. Fine-grained GPU implementation of assembly-free iterative solver for finite element problems[J]. Computers and Structures, 2015,157:9-18.

[8] KARATARAKIS A, METSIS P, PAPADRAKAKIS M. GPU-acceleration of stiffness matrix calculation and efficient initialization of EFG meshless methods[J]. Computer Methods in Applied Mechanics and Engineering,2013,258(5):63-80.

[9] TIAN Jin, GONG Li, SHI Xiaowei, et al. GPU-accelerated fem solver for three dimensional electromagnetic analysis[J]. Journal of Electronics,2011,28(4):615-622.

[10] GAO J, WANG Y, WANG J, et al. Adaptive optimization modeling of preconditioned conjugate gradient on multi-GPUs[J]. Acm Transactions on Parallel Computing,2016,3(3):16.

SU Huia*,QIU Xiaqingb, MA Wenpengb

(a. Network Information and Calculating Center;b. College of Computer and Information Technology, Xinyang Normal University, Xinyang , China)

Abstract:The acceleration of the two-dimensional Laplace equation on the GPU platform was realized by using the finite element method based on Matlab platform. Through the analysis of the physical problem and the construction of the physical model, the generation of the stiffness matrix of the whole CSR format was completed. Parallel programming of the finite element problem was implemented on the Matlab platform using the mixed programming of Matlab and CUDA. Combining with the CuBlas numerical calculation library, the PCG algorithm was used to solve the large linear sparse system of the assembly, and the velocity potential of each lattice can be iterated efficiently. This algorithm not only gives full play to the efficiency of Matlab in numerical calculation, but also gives full play to the advantages of GPU in fine-grained parallel acceleration.

Key words:CUDA programming; finite element method; GPU; preconditioned conjugate gradient algorithm

中图分类号:TP391

文献标志码:A

文章编号:1003-0972(2018)04-0677-04

DOI:10.3969/j.issn.1003-0972.2018.04.031

收稿日期:2018-01-24;

修订日期:2018-07-19;*.

基金项目:国家自然科学基金项目(,61501393);河南省高等学校重点科研项目(17B520034);河南省教育厅教改项目(2017-JSJYYB-055)

作者简介:苏辉(1971—),男,河南平舆人,副教授,硕士,主要从事模式识别和图像信息处理等方向研究.

责任编辑:郭红建