有限元GPU加速计算的实现方法

张健飞　沈德飞

摘要：研究基于GPU的有限元求解中的总刚矩阵生成和线性方程组求解问题.通过对单元着色和分组完成总刚矩阵的生成，并以行压缩存储（Compressed Sparse Row，CSR）格式存储，用预处理共轭梯度法求解所生成的大规模线性稀疏方程组.在CUDA（Compute Unified Device Architecture）平台上完成程序设计，并用GT430 GPU对弹性力学的平面问题和空间问题进行试验.结果表明，总刚矩阵生成和方程组求解分别得到最高11.7和8的计算加速比.

关键词：GPU计算； 有限元法； 刚度矩阵； 预处理共轭梯度法

中图分类号： TB115.7；TP311

文献标志码：B

0 引 言

作为一种求解微分方程或积分方程的微分方法，有限元法[1]以其高度的适应性，成为现代工程设计和结构分析的重要方法之一，并在土木、水利、汽车、机械、航空航天、核工业和大地勘测等众多领域应用广泛.随着科学技术的不断发展，工程问题的规模和复杂程度相应提高，也对有限元计算提出更大规模、更快速度的要求.有限元法的基本思想是“化整为零，积零为整”，与并行计算技术的“分而治之”的基本原则相协调，因此，对于大规模的有限元结构分析，可将效率低下的串行有限元分析改进为并行有限元分析.

自从2006年NVIDIA正式发布用于通用计算的统一计算架构CUDA（Compute Unified Device Architecture）平台[2]后，图形处理器的体系架构得到迅速发展和完善，性能得到大幅提高，使得GPU不仅能高效应用于计算机图形处理，而且其强大的计算能力也能很好地适用于高性能计算[3]，大大推动基于GPU通用计算的研究，并广泛应用于医学成像、生物信息学、计算结构力学、计算流体力学、计算金融、地震勘探、地理信息系统以及电影和动画制作等领域.商业有限元分析软件ANSYS和Abaqus中也应用GPU技术.[4]

目前，基于GPU的有限元分析在求解由有限元生成的稀疏线性方程组的研究中比较多，主要是由于在求解方程组时数据量大而且计算比较集中，便于并行运算.NATHAN等[5]讨论稀疏矩阵的数据结构，并探讨几种有效的基于CUDA的高效的稀疏矩阵与向量相乘的方法；AIL等[6]探讨针对多GPU的基于CUDA的快速共轭梯度法，并探讨共轭梯度法中最耗时的稀疏矩阵与向量相乘的操作；JOLDES等[7]研究基于GPU的以六面体单元为主的混合网格及在有限元中寻求稳定解的问题，并通过CUDA实现基于非线性力学模型的自动模拟神经外科的过程；PAWEL等[8]探讨基于GPU的三维有限元数值积分算法和计算方面的内容；CECKA等[9]探讨基于GPU的有限元法刚度矩阵组装方法，评估每种方法的优缺点；李熙铭[10]验证基于CUDA的复电阻率问题，并详细研究复共轭梯度法；胡耀国[11]运用单元分组的方式计算得出有限元法中的总刚矩阵，并研究基于GPU的共轭梯度法，因当时的资源限制，其在总刚矩阵生成部分的加速效果并不明显.

本文利用GPU强大的并行计算能力和CPU的高效逻辑处理能力，将有限元法中计算量很大的单刚矩阵计算和总刚矩阵生成及线性方程组的求解交给GPU运算，CPU负责相应数据的前处理和后处理及整个分析过程中的逻辑关系.应用CUDA平台完成基于GPU的有限元分析程序，在本文的计算平台上运行，并用弹性力学平面问题和空间问题的有限元求解测试其运行效果.在总刚矩阵生成部分，为避免计算时线程冲突，采用文献[9]和[11-12]中都提到的对结构单元进行着色的方法.在求解大规模线性方程组时充分利用现有资源，在完全调用现有库函数的情况下实现基于GPU的预条件共轭梯度法.

1 CUDA编程模型

NVIDIA的CUDA并行计算模型使用C语言作为开发语言，同时也支持其他编程语言或应用程序接口.一个完整的CUDA程序包含连续运行在CPU端的主程序及用CUDA编写的运行在并行GPU设备上的内核程序.由图1可知，在CUDA架构下，应用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）组成.运行流程：在CPU端准备数据，然后传到GPU端进行并行计算，最后将计算好的数据再传回CPU端.

CPU端的程序主要负责数据的准备和一些逻辑运算，其中Kernel函数是整个应用程序的关键.Kernel函数以线程网格（Grid）的形式组织，每个线程网格由若干个线程块（Block）组成，而每个线程块又由若干个线程（Thread）组成.在执行时，Kernel函数以线程块为单位，各线程块间并行执行，不同线程块间只能通过全局显存进行数据共享，同一线程块内的线程之间可以通过共享内存通信，即在Kernel函数中存在2个层次的并行：在线程网格中的线程块间并行和在线程块中的线程间并行.

在Kernel函数的程序中，每个线程拥有自己的私有寄存器（register）和局部存储器（local memory）；每个线程块拥有一个共享存储器（shared memory）；线程网格中所有线程都可以访问全局存储器（global memory）.还有2种程序中所有线程都可以访问的只读存储器：常数存储器（constant memory）和纹理存储器（texture memory），它们分别为不同的应用进行优化.其中，全局存储器、常数存储器和纹理存储器中的值在一个内核函数执行完成后被继续保持，可以被同一程序中的其他内核函数调用.

2 有限元总刚矩阵

2.1 单元着色

为能够并行计算单元刚度矩阵且不产生线程冲突，通过单元着色方法生成刚度矩阵.单元着色原则是：对于共享同一个节点的所有单元都着色为不同的颜色，即对于任一自由度中的单元没有2个单元的颜色是相同的.

完成对结构中的所有单元着色后，再进行同一颜色的分组.

2.2 总刚矩阵的压缩存储

有限元法生成的总刚矩阵为大型稀疏矩阵，如果使用与稠密矩阵一样的满阵存储法存储，不仅存储量和计算量大，而且会浪费很多不必要的内存空间，需采用压缩存储的方式存储该稀疏矩阵，即只存储稀疏矩阵中的非零元素.在进行单元刚度计算之前，先计算出单元中相互有贡献的节点，对于相互没有贡献的节点，不给予存储空间.压缩存储法中的行压缩存储（Compressed Sparse Row，CSR）法对矩阵的结构没有要求，而且以CSR格式存储的稀疏矩阵能够更好地满足GPU的并行计算.有限元法生成的总刚矩阵通常是无规则的稀疏矩阵，因此采用CSR格式存储.

摘要：研究基于GPU的有限元求解中的总刚矩阵生成和线性方程组求解问题.通过对单元着色和分组完成总刚矩阵的生成，并以行压缩存储（Compressed Sparse Row，CSR）格式存储，用预处理共轭梯度法求解所生成的大规模线性稀疏方程组.在CUDA（Compute Unified Device Architecture）平台上完成程序设计，并用GT430 GPU对弹性力学的平面问题和空间问题进行试验.结果表明，总刚矩阵生成和方程组求解分别得到最高11.7和8的计算加速比.

关键词：GPU计算； 有限元法； 刚度矩阵； 预处理共轭梯度法

中图分类号： TB115.7；TP311

文献标志码：B

0 引 言

作为一种求解微分方程或积分方程的微分方法，有限元法[1]以其高度的适应性，成为现代工程设计和结构分析的重要方法之一，并在土木、水利、汽车、机械、航空航天、核工业和大地勘测等众多领域应用广泛.随着科学技术的不断发展，工程问题的规模和复杂程度相应提高，也对有限元计算提出更大规模、更快速度的要求.有限元法的基本思想是“化整为零，积零为整”，与并行计算技术的“分而治之”的基本原则相协调，因此，对于大规模的有限元结构分析，可将效率低下的串行有限元分析改进为并行有限元分析.

自从2006年NVIDIA正式发布用于通用计算的统一计算架构CUDA（Compute Unified Device Architecture）平台[2]后，图形处理器的体系架构得到迅速发展和完善，性能得到大幅提高，使得GPU不仅能高效应用于计算机图形处理，而且其强大的计算能力也能很好地适用于高性能计算[3]，大大推动基于GPU通用计算的研究，并广泛应用于医学成像、生物信息学、计算结构力学、计算流体力学、计算金融、地震勘探、地理信息系统以及电影和动画制作等领域.商业有限元分析软件ANSYS和Abaqus中也应用GPU技术.[4]

目前，基于GPU的有限元分析在求解由有限元生成的稀疏线性方程组的研究中比较多，主要是由于在求解方程组时数据量大而且计算比较集中，便于并行运算.NATHAN等[5]讨论稀疏矩阵的数据结构，并探讨几种有效的基于CUDA的高效的稀疏矩阵与向量相乘的方法；AIL等[6]探讨针对多GPU的基于CUDA的快速共轭梯度法，并探讨共轭梯度法中最耗时的稀疏矩阵与向量相乘的操作；JOLDES等[7]研究基于GPU的以六面体单元为主的混合网格及在有限元中寻求稳定解的问题，并通过CUDA实现基于非线性力学模型的自动模拟神经外科的过程；PAWEL等[8]探讨基于GPU的三维有限元数值积分算法和计算方面的内容；CECKA等[9]探讨基于GPU的有限元法刚度矩阵组装方法，评估每种方法的优缺点；李熙铭[10]验证基于CUDA的复电阻率问题，并详细研究复共轭梯度法；胡耀国[11]运用单元分组的方式计算得出有限元法中的总刚矩阵，并研究基于GPU的共轭梯度法，因当时的资源限制，其在总刚矩阵生成部分的加速效果并不明显.

本文利用GPU强大的并行计算能力和CPU的高效逻辑处理能力，将有限元法中计算量很大的单刚矩阵计算和总刚矩阵生成及线性方程组的求解交给GPU运算，CPU负责相应数据的前处理和后处理及整个分析过程中的逻辑关系.应用CUDA平台完成基于GPU的有限元分析程序，在本文的计算平台上运行，并用弹性力学平面问题和空间问题的有限元求解测试其运行效果.在总刚矩阵生成部分，为避免计算时线程冲突，采用文献[9]和[11-12]中都提到的对结构单元进行着色的方法.在求解大规模线性方程组时充分利用现有资源，在完全调用现有库函数的情况下实现基于GPU的预条件共轭梯度法.

1 CUDA编程模型

NVIDIA的CUDA并行计算模型使用C语言作为开发语言，同时也支持其他编程语言或应用程序接口.一个完整的CUDA程序包含连续运行在CPU端的主程序及用CUDA编写的运行在并行GPU设备上的内核程序.由图1可知，在CUDA架构下，应用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）组成.运行流程：在CPU端准备数据，然后传到GPU端进行并行计算，最后将计算好的数据再传回CPU端.

CPU端的程序主要负责数据的准备和一些逻辑运算，其中Kernel函数是整个应用程序的关键.Kernel函数以线程网格（Grid）的形式组织，每个线程网格由若干个线程块（Block）组成，而每个线程块又由若干个线程（Thread）组成.在执行时，Kernel函数以线程块为单位，各线程块间并行执行，不同线程块间只能通过全局显存进行数据共享，同一线程块内的线程之间可以通过共享内存通信，即在Kernel函数中存在2个层次的并行：在线程网格中的线程块间并行和在线程块中的线程间并行.

在Kernel函数的程序中，每个线程拥有自己的私有寄存器（register）和局部存储器（local memory）；每个线程块拥有一个共享存储器（shared memory）；线程网格中所有线程都可以访问全局存储器（global memory）.还有2种程序中所有线程都可以访问的只读存储器：常数存储器（constant memory）和纹理存储器（texture memory），它们分别为不同的应用进行优化.其中，全局存储器、常数存储器和纹理存储器中的值在一个内核函数执行完成后被继续保持，可以被同一程序中的其他内核函数调用.

2 有限元总刚矩阵

2.1 单元着色

为能够并行计算单元刚度矩阵且不产生线程冲突，通过单元着色方法生成刚度矩阵.单元着色原则是：对于共享同一个节点的所有单元都着色为不同的颜色，即对于任一自由度中的单元没有2个单元的颜色是相同的.

完成对结构中的所有单元着色后，再进行同一颜色的分组.

2.2 总刚矩阵的压缩存储

有限元法生成的总刚矩阵为大型稀疏矩阵，如果使用与稠密矩阵一样的满阵存储法存储，不仅存储量和计算量大，而且会浪费很多不必要的内存空间，需采用压缩存储的方式存储该稀疏矩阵，即只存储稀疏矩阵中的非零元素.在进行单元刚度计算之前，先计算出单元中相互有贡献的节点，对于相互没有贡献的节点，不给予存储空间.压缩存储法中的行压缩存储（Compressed Sparse Row，CSR）法对矩阵的结构没有要求，而且以CSR格式存储的稀疏矩阵能够更好地满足GPU的并行计算.有限元法生成的总刚矩阵通常是无规则的稀疏矩阵，因此采用CSR格式存储.

摘要：研究基于GPU的有限元求解中的总刚矩阵生成和线性方程组求解问题.通过对单元着色和分组完成总刚矩阵的生成，并以行压缩存储（Compressed Sparse Row，CSR）格式存储，用预处理共轭梯度法求解所生成的大规模线性稀疏方程组.在CUDA（Compute Unified Device Architecture）平台上完成程序设计，并用GT430 GPU对弹性力学的平面问题和空间问题进行试验.结果表明，总刚矩阵生成和方程组求解分别得到最高11.7和8的计算加速比.

关键词：GPU计算； 有限元法； 刚度矩阵； 预处理共轭梯度法

中图分类号： TB115.7；TP311

文献标志码：B

0 引 言

作为一种求解微分方程或积分方程的微分方法，有限元法[1]以其高度的适应性，成为现代工程设计和结构分析的重要方法之一，并在土木、水利、汽车、机械、航空航天、核工业和大地勘测等众多领域应用广泛.随着科学技术的不断发展，工程问题的规模和复杂程度相应提高，也对有限元计算提出更大规模、更快速度的要求.有限元法的基本思想是“化整为零，积零为整”，与并行计算技术的“分而治之”的基本原则相协调，因此，对于大规模的有限元结构分析，可将效率低下的串行有限元分析改进为并行有限元分析.

自从2006年NVIDIA正式发布用于通用计算的统一计算架构CUDA（Compute Unified Device Architecture）平台[2]后，图形处理器的体系架构得到迅速发展和完善，性能得到大幅提高，使得GPU不仅能高效应用于计算机图形处理，而且其强大的计算能力也能很好地适用于高性能计算[3]，大大推动基于GPU通用计算的研究，并广泛应用于医学成像、生物信息学、计算结构力学、计算流体力学、计算金融、地震勘探、地理信息系统以及电影和动画制作等领域.商业有限元分析软件ANSYS和Abaqus中也应用GPU技术.[4]

目前，基于GPU的有限元分析在求解由有限元生成的稀疏线性方程组的研究中比较多，主要是由于在求解方程组时数据量大而且计算比较集中，便于并行运算.NATHAN等[5]讨论稀疏矩阵的数据结构，并探讨几种有效的基于CUDA的高效的稀疏矩阵与向量相乘的方法；AIL等[6]探讨针对多GPU的基于CUDA的快速共轭梯度法，并探讨共轭梯度法中最耗时的稀疏矩阵与向量相乘的操作；JOLDES等[7]研究基于GPU的以六面体单元为主的混合网格及在有限元中寻求稳定解的问题，并通过CUDA实现基于非线性力学模型的自动模拟神经外科的过程；PAWEL等[8]探讨基于GPU的三维有限元数值积分算法和计算方面的内容；CECKA等[9]探讨基于GPU的有限元法刚度矩阵组装方法，评估每种方法的优缺点；李熙铭[10]验证基于CUDA的复电阻率问题，并详细研究复共轭梯度法；胡耀国[11]运用单元分组的方式计算得出有限元法中的总刚矩阵，并研究基于GPU的共轭梯度法，因当时的资源限制，其在总刚矩阵生成部分的加速效果并不明显.

本文利用GPU强大的并行计算能力和CPU的高效逻辑处理能力，将有限元法中计算量很大的单刚矩阵计算和总刚矩阵生成及线性方程组的求解交给GPU运算，CPU负责相应数据的前处理和后处理及整个分析过程中的逻辑关系.应用CUDA平台完成基于GPU的有限元分析程序，在本文的计算平台上运行，并用弹性力学平面问题和空间问题的有限元求解测试其运行效果.在总刚矩阵生成部分，为避免计算时线程冲突，采用文献[9]和[11-12]中都提到的对结构单元进行着色的方法.在求解大规模线性方程组时充分利用现有资源，在完全调用现有库函数的情况下实现基于GPU的预条件共轭梯度法.

1 CUDA编程模型

NVIDIA的CUDA并行计算模型使用C语言作为开发语言，同时也支持其他编程语言或应用程序接口.一个完整的CUDA程序包含连续运行在CPU端的主程序及用CUDA编写的运行在并行GPU设备上的内核程序.由图1可知，在CUDA架构下，应用程序由CPU端的串行程序（host端程序）和GPU端的并行程序（device端程序或Kernel程序）组成.运行流程：在CPU端准备数据，然后传到GPU端进行并行计算，最后将计算好的数据再传回CPU端.

CPU端的程序主要负责数据的准备和一些逻辑运算，其中Kernel函数是整个应用程序的关键.Kernel函数以线程网格（Grid）的形式组织，每个线程网格由若干个线程块（Block）组成，而每个线程块又由若干个线程（Thread）组成.在执行时，Kernel函数以线程块为单位，各线程块间并行执行，不同线程块间只能通过全局显存进行数据共享，同一线程块内的线程之间可以通过共享内存通信，即在Kernel函数中存在2个层次的并行：在线程网格中的线程块间并行和在线程块中的线程间并行.

在Kernel函数的程序中，每个线程拥有自己的私有寄存器（register）和局部存储器（local memory）；每个线程块拥有一个共享存储器（shared memory）；线程网格中所有线程都可以访问全局存储器（global memory）.还有2种程序中所有线程都可以访问的只读存储器：常数存储器（constant memory）和纹理存储器（texture memory），它们分别为不同的应用进行优化.其中，全局存储器、常数存储器和纹理存储器中的值在一个内核函数执行完成后被继续保持，可以被同一程序中的其他内核函数调用.

2 有限元总刚矩阵

2.1 单元着色

为能够并行计算单元刚度矩阵且不产生线程冲突，通过单元着色方法生成刚度矩阵.单元着色原则是：对于共享同一个节点的所有单元都着色为不同的颜色，即对于任一自由度中的单元没有2个单元的颜色是相同的.

完成对结构中的所有单元着色后，再进行同一颜色的分组.

2.2 总刚矩阵的压缩存储

有限元法生成的总刚矩阵为大型稀疏矩阵，如果使用与稠密矩阵一样的满阵存储法存储，不仅存储量和计算量大，而且会浪费很多不必要的内存空间，需采用压缩存储的方式存储该稀疏矩阵，即只存储稀疏矩阵中的非零元素.在进行单元刚度计算之前，先计算出单元中相互有贡献的节点，对于相互没有贡献的节点，不给予存储空间.压缩存储法中的行压缩存储（Compressed Sparse Row，CSR）法对矩阵的结构没有要求，而且以CSR格式存储的稀疏矩阵能够更好地满足GPU的并行计算.有限元法生成的总刚矩阵通常是无规则的稀疏矩阵，因此采用CSR格式存储.