基于Matlab平台有限元方法的GPU加速

2018-10-17邱夏青马文鹏

信阳师范学院学报（自然科学版） 2018年4期

苏辉，邱夏青，马文鹏

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

0 引言

对于很多大规模数值模拟来说，稀疏线性系统的求解占用了大部分的计算时间和资源.随着科学与工程问题的规模和复杂程度的提高，对稀疏性系统的求解规模和速度提出了更高的要求，利用高性能计算技术是当前解决大规模稀疏线性系统求解的重要途径[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加速，算法性能有明显的优越性.

1 求解向量方程组的方法与技术

1.1 有限元方法

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

1.2 稀疏矩阵行压缩存储

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

1.3 稀疏矩阵向量乘算法

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

1.4 预处理共轭梯度算法

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

1.5 Matlab软件的优势

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

2 二维拉普拉斯方程的有限元求解

2.1 求解应用方程

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

(1)

(2)

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

2.2 插值函数及相关计算

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

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

(3)

其中

(4)

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

(5)

2.3 加权余量方程

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

(6)

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

(7)

2.4 单元方程的建立和总体方程的组合

单元方程的形式如下：

(8)

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

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

2.5 求解线性方程组

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

3 基于GPU的二维拉普拉斯方程并行加速

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

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

3.1 以格点为主序装配非零元素

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

3.2 基于Matlab平台的GPU加速

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平台实现对大规模线性方程组的并行求解操作.

4 计算结果及分析

4.1 测试环境与测试算例

本次实验是在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万.