基于GPU的直流电阻法三维有限元数值模拟.pdf

摘要 摘要 直流电阻率三维正演通过数值模拟获得三维地电模型中的电场分布 是三维 反演和资料解释的基础 快速准确的电阻率三维数值模拟一直是直流电法勘探中 的重要问题 有限单元法能够适应于任意复杂的三维电性结构 在电阻率三维数值模拟中 应用广泛 其主要问题是 对于大型细网格 最后形成的大型稀疏线性方程组的 求解计算量巨大 非常费时 严重阻碍了电阻率三维资料的反演解释 利用逐次 对称超松弛 s s o r 预处理共轭梯度法求解三维有限单元法形成的线性系统 计算速度提高很多 但是在处理大型精细网格时 计算效率依然不够 并行计算是解决大量密集计算的有效途径 近年来 基于g p u 的c u d a 并 行计算机体系获得了越来越多的关注 相对于传统的c p u 并行架构 c u d a 架 构在浮点运算能力和内存带宽上都占有优势 开发语言是标准c 语言的扩展 容易接受 同时还提供大量实用高效的库函数协助程序设计 这些特性让c u d a 在面世后很快就得到了广泛应用 但应用于地球电磁三维数值模拟还非常少见 本文针对大型电阻率三维有限元数值模拟问题 提出了基于s s o r 预处理共 轭梯度法解大型线性系统的并行方案 并在c u d a 平台上实现 电阻率三维数 值模拟并行计算结果表明 基于g p u 的c 切 a 并行计算结果正确 算法高效 可获得1 0 倍以上的加速比 关键词 c a 电阻率三维正演有限单元预处理共轭梯度并行计算 摘要 i i a b s t r a c t a b s t r a c t t h f e e d i m e n s i o n a l 3 d f o r w a r dm o d e l i n go fd cr e s i s t i v i t vm e t h o dj su s e dt o o b t a i nt h ed i s t r i b u t i o no ft h ec u r r e n tn e l di n3 dg e o e l e c t r i c a lm o d e l w h i c hi st h e b a s e o f3 一dr e s i s t i v i t yi n v e r s i o na n dd a t a e x p l a n a t i o n ar 印i da 1 1 da c c u r a c v3 d r e s i s t i v i t rm o d e l i n gi sa l w a y sa i li m p o r t a n tp r o b l e mi ne l e c t r i c a le x p l o r a t i o nm e t h o d s i n c et h e3 df i n i t ee l e m e n t f e m e m o dc a nb ea d a p t e dt oa n y c o m p l i c a t e d3 一d g e o e l e c t r i c a ls t m c t u r e i ti sw i d e l y 印p l i e di n3 一dr e s i s t i v i t yn u m e n c a ls i m u l a t i o n t h em 旬o rp r o b l 锄si tc o n f r o n t e dw i t h h o w e v e r a r eh u g e 狮o u n to f c o i n p u t a t i o n a n dt i m ec o n s 眦l i n gt os o l v et h el a 唱e1 i n e a rs y s t e ma r i s e n 行o mt 1 1 e3 一df e m o d e l i n g f o rav e 巧f i n e3 ds n u c t l l r em o d e l w h i c bs e 矗o u s l yi m p e d e st h e3 一di n v e r s i o na n d d a t ai n t e 印r e t a t i o n s y m m e t r i cs u c c e s s i v eo v e r r e l a x a t i o np r e c o n d i t i o n e dc o n i u g a t e 孕a d i e n tm e t 圭l o d s s o r p c g h a sb e e n e m p l o y e d t o s p e e d u p t h e3 df e c o h l p u t a t i o n h o w e v e r t h ee 街c i e n c yi ss t i l ll o ww h i l eh a n d l i n ga l a i 售ef i n eg i r d p a r a l l e l c o m p u t i n g i sa ne f r e c t i v e w a yt o s o l v em i sp r o b l e m r e c e n t l y g p u b a s e dc u d a p a r a l l e la r c h i t e c t u r er e c e i v e dm o r ea 1 1 dm o r ea t t e n t i o n u n l 浓et 1 1 e t r a d i t i o n a lc p up a r a l l e la r c h i t e c t u r e c u d ah a so b v i o u sa d v a l l t a g eo fc o m p u t a t i o n a l h o r s 印o w e ra n dm e m o r yb a n d w i d t ho v e rc p u w i t has l i 曲t l ye x t e n d e dc l a n g u a g e w h l c hl se a s yt om a s t e ra n dn u m e r o u sl i b r a r ys u p p o n n l e s ef e a t u r e sm a k e i tq u i c k l y w i d es p r e a da f t e ri t s a p p e 跏1 c e b u tt h e r ei sl i t t l eg p u b a s e dp a r a l l e lc o m p u t i n g c 锄e do u ti 1 1g e o e l e c 仃o m a g n e t i c3 d m o d e l i n gs of h t h i sp 印e rd i s c u s s e st h ep a r a l l e l i z a t i o no fs s o rm e t h o de v o l v e di nl a r g es c a l e 3 一dr e s i s t i v i t y 姗m e r i c a ls i i n u l a t i o n w h i c hw o r k so u to nac u d a p l a t f 0 肌 t h e r e s u l to f3 dr e s i s t i v i t yp a r a l l e ln u m e r i c a ls i m u l a t j 锄s h o w st h a tg p u b a s e dc u d a p a r a l l e lc o m p u t a t i o ni so fm g ha c c u r a c ya j l de 伍c i e n c y a n dg a i n sa no v e rt e nt i m e s s p e e d u pb e n e f i t k e yw o r d s c o m p u t eu n i f i e dd e v i c ea r c h i t e c t u r e c u d a 3 一df o r w a r dm o d e l i n g f i n i t ee l e m e n t p r e c o n d i t i o n e dc o n j u g a t eg r a d i e n t p a r a l l e lc o m p u t i n g 订i a b s t r a c t v 目录 强录 第一章 第二章 2 1 2 2 第三章 3 1 3 2 3 3 3 4 第四章 4 1 4 2 引言 1 直流电法有限单元三维正演 3 地中稳定电流场的基本性质 3 2 1 1 基本性质 3 2 1 2 边界条件 4 有限单元法求解 4 s s o r 预条件共轭梯度算法 7 共轭梯度法 c g 7 预处理共轭梯度法 p c g 8 稀疏矩阵压缩存储 1 0 s s o r 方法伪代码和流程优化 1 l c u d a 并行算法 15 g p u 通用计算 1 6 c u d a 编程基础 1 8 4 2 1 主机和设备异构执行模式 1 8 4 2 2 线程组织模式 18 4 2 3 c 切 a 存储模型 2 0 4 3c u d a 库函数 2 1 4 4 k e m e l 函数设计 2 3 第五章基于g p u 的电阻率三维有限元并行计算实现 2 5 5 1 正确性验证 2 5 5 2 并行程序性能评测 2 8 第六章结论 3 3 参考文献 3 5 致谢 3 7 在读期间发表的学术论文与取得的其他研究成果 3 9 目录 一 一 v t 第一章引言 第一章引言 直流电法是电法勘探的重要分支之一 把直流电源和大地相连 产生稳态电 流场 由于地下岩石和矿物导电性质存在差异 使得这一电流场呈现出各种不同 的分布 通过研究这一电流场 获取地下岩层的电性结构 从而达到找矿目的和 解决其他地质问题 随着地质勘探目标越来越复杂 要求越来越精细 三维电法勘探及资料解释 也越来越获得关注 然而 众所周知 快速 精确的电阻率三维正演是三维反演 解释的基础 直流电法三维正演的方法一般没有解析法 只能通过物理模拟方法 或数值模拟方法获得三维电流场的分布 物理模拟方法是建立物理模拟系统进行 实验观测 对实验设备和场地有较高要求 同时对某些模型 还面临着施工困难 成本也高 数值模拟方法是通过电子计算机 利用偏微方程数值方法计算出电场 分布的数值解 随着计算机技术的快速发展 数值模拟技术成为研究复杂模型电 磁场分布的主要手段 目前数值模拟中常见的方法有边界单元 积分方程 有限差分和有限单元等 方法 三维边界单元 积分方程法由于其仅需在异常体表面进行剖分 实际是降 维计算 计算简单 速度 在电磁三维数值模拟中应用广泛 但其仅限于简单 有限个异常体模型的数值模拟 当地下结构复杂 异常体很多时 边界单元 积 分方程法就很繁杂 因而在三维反演解释中应用不多 有限差分和有限单元均适 应于复杂的结构模拟 因而有广泛的应用 有限差分的思想是把求解区域剖分为 网格 在网格的节点上用差分方程来近似待求微分方程 从而得到微分方程的近 似解 如果加密网格 差分方程对微分方程的近似更好 近似解的精度也会随之 提高 有限差分法由于其实现简单 快速 在很多领域占据主流地位 但是在处 理复杂几何形状的问题时显得力不从心 求解精度降低甚至失效 有限单元法通 过变分原理或者伽辽金法将偏微分方程转换为代数方程 然后对线性方程组求 解 有限元单法在适应复杂电性结构 特别是地形模拟方面较有限差分法有明显 优势 吴小平 2 0 0 3 因而有限元法在电阻率数值模拟中得到了更广泛的应用 z h o ua n dg r e e l l l l a l g h 2 0 0 1 吴小平 2 0 0 3 王威 2 0 1 0 不过其计算量大 特别在处理复杂地形和大尺度问题时巨大的计算量成为其应用的明显障碍 吴小 平 2 0 0 3 提出了快速的算法 使得计算速度较之前的方法极大提高 1 0 倍 内存的需求也大大减小 然而 由于实际问题的复杂性 计算量仍然十分巨大 并且随着对问题的研究不断深入 还产生了更加要求计算能力的新课题 如对精 细结构的研究等 为了适应研究课题的新要求 迸一步提高数值模拟的速度 将 密集计算并行化是行之有效的方法 也是该领域发展的趋势 1 第一一章引言 c u d a c o m p u t eu n i n e dd e v i c ea r c h i t e c m r e 是n v i d i a 公司 英伟达公司 于2 0 0 6 年1 1 月正式推出的并行软硬件体系 作为一种新型的并行架构 获得了 广泛关注 它区别于以往的并行计算机的显著特点是 不再采用c p u c e n t r a l p r o c e s s i n gu n j t 作为主要计算部件 而是用g p u g r a p h i cp r o c e s s i n gu n i t 取而代 之 相对于更侧重于处理复杂指令流的c p u 专门设计用来进行图形图像渲染 的g p u 在浮点运算能力和内存访问带宽等方面均占有优势 这让以g p u 为运算 部件的c u d a 架构具备了强劲的计算能力 c a 平台采用一种类似于标准c 语言的编程语言 c u d ac 进行开发 同时很好的掩盖了硬件结构对并行算法设 计的限制和要求 具有很强的通用性和友好性 另外 组建c u d a 平台的费用 要少很多 具有明显的经济性 性能出色 易于使用 成本节约的c u d a 架构 很快就得到了广泛应用 包括分子动力学 a n d e l s o ne t a 1 2 0 0 8 计算生物学 l i u e t a 1 2 0 0 7 量子蒙特卡罗 a n d e r s o ne t a 1 2 0 7 流体力学模拟 r a n d v i ka n d p u l l a l l 2 0 0 7 c t 图像重组 e t a 1 2 0 0 8 地震分析 k o m a t i t s c h ae t a 1 2 0 0 9 光线追踪 g a r a l l z h aa 1 dl 0 0 p 2 0 1o 等诸多大规模计算领域 目前 代表世界尖端 水平的超级计算机中也开始大量采用g p u 提高计算能力 如我国的 天河一号 超级计算机系统 等 本文也选取g p u 作为计算工具 探讨对应的并行算法 希望以此提高直流电阻率法三维正演的计算效率 为快速 精细的三维反演奠定 基础 第二章直流电法有限单元三维正演 第二章直流电阻法有限单元三维正演 由前章的讨论可知 采用直流电法进行勘探工作时 在地下建立起了稳定电 流场 分析这一电流场的分布特点和规律 是进行三维正演的理论基础 2 1 地中稳定电流场的基本性质 2 1 1基本性质 地中稳定电流场具有以下三个主要性质 由微分形式表示 1 地中电流密度与电场强度成正比 以j 记电流密度矢量 以e 记电场强度矢量 那么有 j 栅 里 1 1 p 在各向同性介质中 p 是标量 2 地中电流连续 以 记电流源的电流强度 那么包含这一电流源的任意闭合曲面s 的通量为 中j n 嬲 1 2 s 其中n 是面元搬的单位法线矢量 如果闭合曲面s 不包含电流源 那么上 述积分等于0 这便是电荷守恒定律 说明任意的闭合曲面内 不存在正电荷或 负电荷的不断积累 也就是说 在稳定电流场中电流是连续的 其微分形式如下 d i v j o 1 3 3 稳定电流场的势场性 稳定电流场不随时间变化 是一种势场 定义其中任意一点p 的电位u 等 于把单位正电荷从该点移至无穷远处时 电场力所做的功 即有 u ie d l 1 4 o 那么有关系式 e 一g r a d u 1 5 同时 该势场是无旋场 绕任意闭合回路 电流场所作的功恒等于零 即有 爹e d l o 1 6 第二章直流电法有限单元三维正演 其微分形式为 2 1 2 边界条件 有以下三类边界条件 第一类 第二类 第三类 2 2 有限单元法求解 r o t e 0 u o r u 篆p o 2 万尺 k 一土型 o n i i 烈 n 咖 u u 2 j n j 2 n e t t e 2 旦 咝 仍t a l l 鼠 1 7 1 8 1 9 1 1 0 1 1 1 1 1 2 1 1 3 1 1 4 有稳态电流场的三个基本性质可以得到下式 v o v v x y z 一 万 r r 0 1 1 5 边界条件可以改写为 其中b r r 0 r p r r o 由迦辽金方法 此边值问题可以化为等价的积分形式如下 l 万v v 6 v v 万 r b d q 一 r 万v 6 v v n 二三旦掣v d i 1 1 7 事 第二章直流电法有限单元三维正演 分邵枳分得剑 l 万v v 6 v v 矗q 一l v 万v 7 a v d q l 万v 6 v v n d 1 1 1 1 8 将 1 1 8 式代入 1 1 7 得到 l v 西 丁 6 v 谢q 一上西 万 r r 0 d q r 万v 竺掣v d r o 1 1 9 由 1 1 9 式得到泛函极值问题 风 v v r 咿v v p 一上v 瞰r r 0 牌 t 掣v 2 d r 1 2 卯 v o 离散化得到 f d 喜l v v 6 v v d q 一喜l v 2 万 r r o d q 莓l 又有 1 2 1 1 m v n r v v r n 1 2 2 其中 i 言 1 参考 1 碾7 7 1 六f 磊 碾 缶 是结点f 的坐标 对于单元的中心 y c z c 有 把 1 2 2 式代入 1 2 1 式后 有 芒2 言2 一l x t 以 2 7 72 y 一以 d p 2 2 一 z z c c 1 2 3 第二章直 赢电法有限单元三维正演 f v l v n7 v 6 v n 7 v d q 一2 如 r 0 r 掣 n r n r w r 芝v 丁 l v n r r 6 n r d q v 刮r o 莓v r 但鼍型川 i 1 2 4 羔v e k v 一2 如 r o v r k v 艺v r 西一2 r 0 v r 两 p 1 r v 丁k 1 v 一2 v r p v 丁k 2 v v l v 一2 v r p 其中p 0 o 0 o o 0 7 再由卯 v o 得到 k v p 求解此方程 即可求得每个节点的电位 又由 k j 二 v n 丁 1 6 v n r d q 和 1 2 5 1 2 6 砼 k q 警警 呸 等警 吒 警警 慨尝掣 吒 譬譬蝎 警譬 1 2 7 q 2 上 吒2 j 上 a 3 2 i 土 1 z o x0 v鲫鲫0 z0 y 帆警警 吒 警警怕 警警f c j co z c o zo zo z 可知 k 是8 8 的对称矩阵 其下三角矩阵可以表示如下 趔u 琏 趔 1 趔 趔 霹 丁 m c 1 2 8 第三章s s o r 预条件共轭梯度算法 第三章s s o r 预条件共轭梯度算法 通过有限单元方法 最后得到线性系统 么x 6 3 1 并且随着网格格点的增加 彳的维数也在增加 如何正确快速的求解这一大 型线性系统是一个重大课题 关于这一问题现在已经有了很多求解方法 其中共 轭梯度法 c o n j u g a t eg r a d i e n t 由于其高度的计算效率 良好的计算稳定性而得 到了越来越广泛的应用 本文沿用共轭梯度法的基本思路 结合计算三维电导率 的实际情况 引入对称超松弛预条件共轭梯度 s y m m e 仃i cs u c c e s s i v e o v e r r e l a x a t i o np r e c o n d i t i o n e dc o n j u g a t eg r a d i e n t 迭代算法 p 印a 1 r a k a k i s l9 8 6 并 发展分别与之对应的基于g p u 的并行算法 在保证计算精度的同时 大大提高 了计算速度 3 1 共轭梯度法 c g 共轭梯度法是由h e n s t e n s e 和s t i e f e l 1 9 5 2 提出的迭代算法 几十年问不 断发展改善 在包括地球物理等众多研究领域有着广泛的应用 其核心思想是通 过目标函数的梯度 构造共轭方向 并将其作为搜索方向 逐步到达目标函数的 极值点 这里的目标函数为 妒 x 去x 丁血一6 丁x 3 2 其极值点满足 v 烈x 血一6 o 3 3 也就是线性系统 3 1 的解 其具体迭代过程如下 首先初始化 26 一瓴 3 4 p o2 然后进行迭代 第三章s s o r 预条件共轭梯度算法 l o r z 2 u l 2 口 兰 l b 却i 2 z i 口i 3 5 i 12l g 却i 屈 鼎 当系数矩阵爿是对称正定矩阵时 共轭梯度法一定是收敛的 这正是其稳定 性的保证 但是在处理一般的矩阵时 由于爿的条件数很大而导致收敛变慢 严 重影响计算速度 为了解决这一问题 一种方法是把么进行完全c h o l e s k y 分解 即 么 口 3 6 其中工是下三角矩阵 这 做法带来新的问题 那就是这样的分解计算同样 耗时 同时 即使在么是稀疏矩阵的情况下 分解得到的三也不再是稀疏矩阵 非零元素的增加将会带来额外的计算操作和内存开销 3 2 预处理共轭梯度法 p c g 为了克服共轭梯度法对一般矩阵收敛缓慢的缺点 通常的做法是对系数矩阵 进行预处理 降价其条件数 然后再对等价问题使用共轭梯度法求解 这种方法 称为预处理共轭梯度法 p r e c o n d i t i o n e dc o n j u g a t eg r a d i e n t 预处理方法灵活多 样 对不同的问题有不同的思路 本文采用对称超松弛预条件共轭梯度 即s s o r 方法 其核心思想是选取合适的预优矩阵m c c r 使得c 1 么 c 丁 1 的条件数较 小 m 1 么近似于单位矩阵 把待求的线性方程组改写成m 血 m 1 6 再运用 c g 方法求解 s s o r 方法有多种预优矩阵的选l 汉办法 如s a a d 2 0 0 0 的分解如 下 m d 加 d 一1 d 劬f c c l 3 7 其中d 是对角阵 由么的对角线组成 e 和f 为三角阵 分别对应彳的严 格下三角部分和严格上三角部分 缈是松弛因亍 又如k l a u s 1 9 9 5 采取的分解 方式 m e j 纪 f 3 8 其中 是单位矩阵 e 和f 满足下列关系 第三章s s o r 预条件共轭梯度算法 彳 e f e f 3 9 确定预条件矩阵m c c 丁 一1 后 代入c g 方法的迭代流程 得到s s o r 的 计算方法 如下 首先初始化 嘉 b 呐 风 c c l 1 然后进行迭代 声r 待0 1 2 倔 丛竺 盛2 1 p 钯 篡i 二麓 b l 12 一 鸽 层 锗 a 1 c c r 1 层只 与c g 方法比较 s s o r 方法的迭代过程中出现了 c c r 一一项 而在大型 矩阵运算中 求逆的运算是困难和极为费时的 所幸这里需要计算的仅仅是 c c r 1 与一个向量的乘积 于是可以通过求解两个线性方程组来避免矩阵求逆 的运算 计算过程如下 令 g f a 伊 巧 3 1 2 两边同时乘以c 伊 得到 c 吼 i 3 1 3 再令 e 吼 3 1 4 代入 3 13 式得到 颐 l 3 1 5 首先求解 3 1 5 式 得到 然后再求解 3 1 4 式 得到g i 即为所需盼向量 虽然求解线性方程组的计算量也较大 但是与求逆运算比起来还是小量 同时这 里的两个线性方程组的系数矩阵都是三角阵 在一定程度上也减少了运算量 简 化了计算过程 第三章s s o r 预条件共轭梯度算法 可以看出 s s o r 方法的单次迭代时间是长于c g 方法的 其优势在于减少 了迭代所需的步数 所以能否取得较好的速度提升 关键在于么 1 的近似 c c 7 一 的好坏 3 3 稀疏矩阵压缩存储 在本文的三维电导率计算实例中 得到的系数矩阵彳均为大型对称稀疏矩 阵 显然完全存储么是不可取的 不仅极大的浪费内存资源 造成内存溢出的风 险 而且增加大量的对零元素的操作 这些操作对计算结果毫无影响 属于无用 的运算 严重影响计算的效率 因此通常都会对大型稀疏矩阵进行压缩存储 只 记录非零元素的值和位置信息 压缩存储的方式多种多样 各有优劣 本文结合 g p u 并行计算的实际情况 采用行压缩格式 c o m p r e s s e ds p a r s er o wf o 衄a t c s r 这一格式具有压缩存储的一般性质 同时也是n v m 认公司提供的 c u s p a r s e 程序库的数据接口类型之一 可以为程序实现带来很大的便利 节 省开发时间 提高开发效率 下面对其进行说巨目 一个咒 2 的稀疏矩阵彳的c s r 格式由以下的参数构成 注意 按c 语言的 惯例 下列数组和矩阵的下标都是从0 开始计数的 1 整型变量1 1 i l z 等于彳的非零元素个数 2 实型数组v a l 维数是如z 的数组 按行优先顺序存储么的非零元素的值 实际程序中 彳的对角线上的零元素也存储在其中 换而言之 对角线上的所有 元素都是作为非零元素来看待的 3 整型数组r o wp 仃 维数是甩 1 的数组 其中最后一个元素的值等于姗z 前面的咒个元素 存取的是对应行的首个非零元素在v a l 数组中的偏移量 也就 是下标 例如r o wp 仃 2 的值是第3 行 矩阵的行号从0 开始计数 的首个非零 元素在v a l 数组中的位置 4 整型数组c o l i n d 维数是1 1 1 1 z 的数组 每个元素是v a l 数组中对应元素 的列标 下面用一个4 4 的矩阵彳为例子进行说明 么 5 17 7 0 oo o 7 7o o3 60 o 0 03 6 2 72 3 0 oo 02 3一o 9 那么n n z 的值为1 0 对角线中的零元素也计算在内 三个数组的值如表3 1 第三章s s o r 预条件共轭梯度算法 表3 1按c s r 格式存储矩阵月 v a l5 17 77 7o o3 63 6 2 72 32 3 0 9 r o w p t r o2581o c o li n do1o1212323 在程序设计中 对么的对角线的操作很多 快速的对么的对角元素寻址是非 常重要的 对于一般情况 定位按照c s r 格式存储的矩阵的对角元素是比较麻 烦的 每定位一个对角元素 需要在c o l j n d 数组的一部分中进行 次搜索操作 同时考虑到么的对称性 本文在实际的程序设计中 仅存储了彳的下三角部分 包 括对角线 以上面的矩阵为例 实际保存的矩阵e 如下 e 5 10 00 00 0 7 7 0 0 0 o0 0 o o3 6 2 7 o o o 0o 02 3 0 9 此时n n z 的值为7 三个数组的值如表3 2 表3 2 按c s r 格式存储矩阵e v a l5 17 70 o3 6 2 72 3 0 9 r o w j t r 01357 c o li n d 001l22 3 对于这样的三角阵 对角元素的定位变得简单 有如下关系式 研司 珏删一p 纱p 1 一1 3 1 6 由此看出 选取合适的数据存储方式 可以为程序设计带来方便 在后文中 还将看到 这一存储方式还可以简化部分迭代流程 3 4 s s o r 方法伪代码和流程优化 结合c s r 压缩存储方式 现在可以得到s s o r 预处理共轭梯度法的伪代码 这里的预优矩阵按以前面所述的s a a d y 的方法为例 首先要决定数据存储结构 在s s o r 的迭代流程 3 1 1 中 既有对称矩阵彳出 现 也有下三角矩阵c 出现 按前面的讨论 对于彳只需要存储其下三角部分e 即可 有 彳 e d e 3 1 7 同时注意到c 的定义 3 7 式中 c 的形式稍显复杂 不妨令 u d e 3 18 用己 来代替c 参与计算流程 现在有了两个下三角矩阵 和u 出于节省 1 1 第三罩s s o r 预条件共轭梯度算法 内存开销的考虑 同时结合库函数调用的数据接口的实际情况 采取只存储u 的 方案 那么现在需要用u 来表示e 这很容易得到 由 3 1 8 式 有 e 旦一旦 f 3 1 9 1 缈 把 3 19 式代入 3 17 式得到 么 旦 1 一三 d 斗堡 3 2 0 倒甜国 另外在计算m 1 一项的时候 增加了一步对角矩阵和向量的乘法 按下面 的步骤进行 仍然令 两边同乘矩阵m 得 令 代入 3 2 1 式 得 g m r 岣 凹一1 u r g r d 1 矿g 3 2 1 3 2 2 3 2 3 3 2 4 求解这一线性方程组 得到 然后在 3 2 3 式两边同乘对角矩阵d 得 砚 扩g 3 2 5 不妨令f d f 求解 3 2 5 便得到了g 通过上面的讨论 确定了只存储下三角矩阵u 的存储方案 并得到了其它矩 阵与u 的关系 那么可以得到s s o r 方法的伪代码 如图3 1 现在对s s o r 方法的流程进行必要的优化 庄意在循环体中对向量p 的更新 p f l g 州 a 3 2 6 那么后面计算却的一步可以写成 却 彳g f l 4 3 2 7 把爿和u 的关系式 3 2 0 代入 3 2 7 得 却 旦 1 一三 观 竺 却i 3 2 8 在 3 2 8 式中 后两项在前面已经算出来了 在计算吼 m 巧 的过程中求 11 第三章s s 0 r 预条件共轭梯度算法 得了u 丁阢 z 而爿p 在上一次循环中已经算出 采用 3 2 8 式代替原来对却川 的计算 减少了一次矩阵和向量的乘法运算 而代价是增加了一次向量与向量的 加法 从而减少了运算的总量 而且省掉的u g 这一运算 是比较费时的一步 因为只存储了u 所以在矩阵乘以向量的运算之外 还有额外的对转置矩阵的元 素定位的操作 再加上u 是以c s r 格式压缩存储的 其元素与转置矩阵的元素 之间对应关系也变得复杂 也增加了程序编写的难度 3 2 8 式这一替代 既减 少了计算量 又节省了开发时间 为整个工作带来了很大的方便 第三章s s o r 预条件共毫巨梯度算法 x x 跳i n l i z a t i o n r 6 一 旦 1 一三 陕 叫叫 p r 三2 r j d p 矾l r 2 d s o v eu tq t 2 ig e tq m r 2g 却 旦p 1 一三 印 竺p 国 缈 谢6 幻 厂 g r d o t q p r e r d o t q r n o t a 口 o p 却 x x o p r 一口却 胛 厶 厂 w 五f 彪 里 占 p j d 舰地 r 乞 d s o l v eu tq t 2 嚏e tq m r 崩6 幻 g 9 r d o t q r d o t q p r e r d o t q p r e r d o t q p q pp 却 旦p 1 一三 印 堡p 国 国 j r d o t q 口 o p 和 x x a p r r 一口却 7 6 一彳x o 图3 1s s o r 方法伪代码 1 4 e臼 c口cx 旷i 第四章c u d a 并行算法 第四章c u d a 并行算法 n v i d i a 公司开发的c u d a 架构 虽然在并行计算领域的出现得较晚 但 是凭借强劲的性能和丰富的技术支持获得了越来越多的关注 并且在很多领域有 着广泛的应用 如磁流体力学 w o n ge ta 1 2 0 1 1 材料破坏研究 融h aa 1 1 ds m i d 2 0 1 1 化学反应扩散数值模拟 m o l n 缸e ta 1 2 0 1 1 计算生物学 l i ue ta 1 2 0 l o 空气污染传播模式研究 m o l n 缸e ta 1 2 0 1 0 量子计算 g u t i 6 r r e ze ta 1 2 0 1 0 云计算 z h a n g 2 0 1 2 等 在地球物理领域 也涌现出很多采用g p u 并 行计算而获得性能提高的例子 如地震波传播数值模拟 k o m a t i t s c he ta 1 2 0 0 9 地震成像 c h e ne ta l 2 0 1 2 叠前时间偏移 s h ie ta l 2 0 1 1 地震三维有限差 分模拟研究 z h o ue ta 1 2 0 1 2 等 基于g p u 的并行计算开始变得普及 首先源自于g p u 并行系统强大的计算 性能 在计算机图像处理领域 像高分辨率三维实时图像演算这种强烈的市场需 求 推动g p u 的发展成为了载有多核 能够进行多线程操作的并行系统 从图 4 1 和图4 2 n v 认 2 0 1 1 中可以看出 g p u 相对于c p u 在浮点运算能力 和内存带宽方面都占有很大优势 翻鞫婚潸翳蛹l 翻麓糖海 图4 1g p u 和c p u 每秒浮点运算操作数量对比 第四章c u d a 并行算法 图4 2g p u 和c p u 内存带宽对比 另一方面 在为通用计算目的而设计的c l d a 架构上面可以使用标准c 语 言进行开发 而标准c 语言广为人知 于是在熟悉c 语言的各行各业人员中推 广c u d a 并不困难 在此基础上 c u d a 还通过第三方的支持 实现了与f o n r a n c j a v a o p e n c l 和d i r e c t c o m p u t e 等其它高级语言和开发工具的对接 c u d a 的可推广人群更为广泛 更进一步 为了满足来自不同领域的客户需求 n v d n 和第三方公司还推出了很多的库函数和工具程序包 例如c u f f t c u b l a s n p p c u d p p p h v s xp h v s i c s 多样的语言支持和丰富的库函数 让c u d a 成 为一个方便易用的平台 并且还有越来越多的第三方公司加入进来 本文利用g p u 在并行计算方面的技术优势 提出直流电法三维电阻率的 c u d a 并行方案 配合使用内核函数和库函数完成并行程序 在保证计算精度的 同时取得了较好的加速效果 4 1g p u 通用计算 g p u 的设计初衷并不是用于通用计算 而是用来辅助c p u 完成复杂的图形 图像渲染等计算机图形学领域的应用 g p u 从 叻能单一的计算机组成部件 发 展到现在的通用计算平台 与c p u 的发展和并行计算的兴起密切相关 在过去很长的时期内 对计算机性能的提高的主要手段一直是提升处理器的 运算速度 随着c p u 的主频从m h z 量级飞跃到g h z 量级 其性能也有飞跃式 第四章c u d a 并行算法 的进步 然而这一提速手段受到物理法则的限制 当半导体芯片制造工艺已经逼 近芯片尺寸的物理极限时 再也不能通过在单位面积的电路板上增加晶体管的手 段来提高性能 芯片制造厂商不得不寻求另外的方式来增强芯片的运算能力 既 然单个芯片的处理能力已经趋近于极限 那么只能从增加芯片数量着手 这同样 也是一个有效的途径 在2 0 0 5 年出现了双核的c p u 这表明 对计算机性能的 提升 开始走上了以数量取胜的道路 这一发展趋势越来越明显 在双核c p u 之后 涌现出了三核 四核乃至八核的c p u 迅速的占据了主流地位 到现在 单核的c p u 已经不多见了 很显然 多核c p u 的优势就在于可以把任务分配到 各个核心上并行处理 采用并行化设计的软件也大量涌现出来 并行计算开始变 得普及 随着c p u 的发展 很多应用成为可能 这些应用又反过来对计算机性能提 出更高的要求 也正是这样 才推动了整个计算机行业的不断前进 在这些新的 需求中 对于图形图像处理的需求非常突出 其本质根源为视觉信号是我们人类 获取信息的最有效的手段 这种强烈的需求推动g p u 经历更为迅猛的发展 快 速 高质量的图形图像处理意味着巨大的计算量 所以g p u 更注重提高计算性 能 从图4 1 可以看出 g p u 的计算能力 一直领先于同时代的c p u 和c p u 的发展趋势一样 g p u 在提升运算速度的同时 也开始增加处理核心数量来增 强性能 同样也支持并行计算 g p u 的计算性能十分优秀 于是变有人想利用这一强大的计算资源 让g p u 进行非图形图像渲染运算 来解决其它类型的问题 即让g p u 用于通用计算 这一想法很有创造力 但是实现起来相当困难 g p u 并不是为通用计算设计的 要利用其中的芯片进行运算 必须先把待计算的问题转化为图像处理的问题 然 后通过d i r e c 或o p e n g l 等图像处理的应用程序接口 在g p u 上处理转换得 到的图像处理问题 最后还要将结果转换到原问题所需的形式 这对非图形处理 背景的研究人员来讲 不仅要解决自己本行的问题 还要成为一个图形处理的专 家 无疑是过于沉重的负担 因此 在c u d a 问世之前 一直只有少数人在利 用g p u 进行通用计算 这一状况在c u d a 出现之后有了彻底的改变 这一崭新的g p u 架构内嵌了 指令通道 可以让算术逻辑运算单元 撕t h m e t i c1 0 9 i cu 1 1 i t a l u 执行非图形图像 的指令 同时开放了g p u 自有的内存 也就是通常所说的显存 的存取权限 这为g p u 通用计算带来了极大的方便 再也不需要借助额外的应用程序接口 用户能直接使用g p u 来解决问题 同时 如前面所提到 n v i d i a 采用了标准c 语言的扩展来作为c u d a 平台的开发语言 在让许多熟悉c 语言的用户 只需 要进行不多的额外学习 便能使用g p u 强大的并行计算能力 而且 在n v i d i a 第四章c u d a 并行算法 的市场营销策略推动下 c u d a 架构的g p u 容易得到 即使是用于传统目的的 显卡 也支持c u d a 基于g p u 的并行计算变得随处可见 4 2c u d a 编程基础 尽管是一种全新的架构 c u d a 却很好的沿用了为大众熟知的开发方式 对 于c 语言的用户来说 可以用过去所熟悉的方式来工作 从而把精力集中在解 决问题上 而不是花在适应新环境方面 4 2 1主机和设备异构执行模式 在c u d a 程序设计中 用主机 h o s t 来表示c p u 用设备 d e v i c e 来表示 g p u 所谓异构执行 是指一个程序 部分在主机 部分在设备 这两种不同的 架构上执行的模式 借用标准c 语言中的函数口u n c t i o n 概念 把在设备上运行 的函数称为内核函数 k e m e l 这样可以把一个c 切 a 程序看成许多主机上的函 数和设备上的内核函数的顺序组合 其执行方式与普通的c 程序一致 在这种执行方式中 有一点是c u d a 的特点 内核函数在设备上是并行执 行的 g p u 强大的计算性能 就来自于其并行计算的能力 那么一个c u d a 程 序的流程大致是这样 c p u 处理串行的部分 到了需要并行处理的部分 先把 数据传送到g p u 的显存里 然后g p u 开始工作 g p u 在运行好并行的部分后 把结果传送回主机端的内存 程序的执行再回到c p u 见图4 3 n v i d n 2 0 1 1 图中左侧即为程序执行流程 右侧部分则刻画了 并行和串行执行的不同情况 每 个箭头就表示一个线程 异构执行这一方式还有一个很好的运用 当g p u 开始 工作以后 通常都是大量的密集计算任务 所以要花上一定时间 c p u 可以利 用程序流程还没有返回的这段空闲去处理额外的任务 这样便实现了主机和设备 的任务并行 4 2 2 线程组织模式 线程 t h r e a d 是并行程序设计中的常用概念 在并行算法设计阶段 可以 把一个线程看成是一个虚拟的处理器 在这一阶段 最主要考虑的就是如何把任 务分配给各个线程 并合理安排线程结构 以达到较好的效率 在c u d a 架构中 线程有着两种层次的组 织结构 首先是由线程构成线程 块 b l o c k 然后由线程块构成线程网 g r i d b l o c k 和g r i d 的维度都是三维 分别对应不同维度的问题 它们的具体大小是由用户定义 不同的组织结构对并 第四章c u d a 并行算法 行程序的效率有所影响 这也是算法设计阶段需要考虑的重要问题 描述一个线 程块要用到几个c u d a 内建的参数 首先是b l o c k d i m 描述线程块的尺寸 由 于线程块是三维的 所以有三个分量b l o c k d i m x b l o c k d i m y b l o c k d i m z 其 次是b l o c k i d x 描述线程块在网格里的编号 由于网格也是三维的 所以也有 b l o c k i d x x b l o c k i d x v b l o c k i d x z 三分量 一个二维网格见图4 4 n v i d n 2 0 1 1 对线程来讲 不存在维度的描述 但是也同样有三个分量的编号 即t 1 1 r e a d i d x x t h r e a d i d x y 缸e a d i d x z o 对网格 则和线程块一样 有描述尺寸的鲥 1 d i m x 酣d d i m y g r i d d i m z 和描述编号的舒d i d x x 酊d i d y 酣d i d z o 网格的尺寸 是有限制的 随硬件设备变化 图4 3 异构执行示意图 图4 4 线程结构示意图 1 9 第四章c u d a 并行算法 这种多维的组织网络对设计算法很有帮助 例如 想要实行矩阵加法的并行 算法 由于矩阵每个元素在参与加法运算时都是独立的 那么一种很自然的思路 是用一个独立的线程 也就是一个独立的虚拟处理器来完成一个元素的加法运 算 然后同样自然的采用二维的线程块和网格来组织这些线程 如果只能是一维 的线程块 对线程的定位等操作都变得麻烦 在图4 3 中可以看到 k e m e l 函数在g p u 上的执行是多个线程同时进行的 更细节的情况是 c u d a 架构中是按线程块为单位来执行线程的 这里只停留在 逻辑执行的阶段 映射到硬件 则是以更小的单位来执行 所以程序设计时要 考虑如何合理分配线程块的大小等因素 而不用考虑具体执行的硬件有多少处理 核心等细节问题 这一方式一方面降低了程序设计的难度 另一