（计算机系统结构专业论文）基于gpu平台的klu并行算法的研究：对角线块的lu分解.pdf

i 原创性声明 本人声明：所呈交的学位论文是本人在导师的指导下进行的研究工作及取得的研究 成果。除本文已经注明引用的内容外，论文中不包含其他人已经发表或撰写过的研究成果， 也不包含为获得凼苤直太堂及其他教育机构的学位或证书而使用过的材料。与我一同工作 的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意。 学砬论支作者签名：j i 丕玺箜五狸 日 期：业lz ：纶 在学期间研究成果使用承诺书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，即：内蒙古大学有权 将学位论文的全部内容或部分保留并向国家有关机构、部门送交学位论文的复印件和磁盘， 允许编入有关数据库进行检索，也可以采用影印、缩印或其他复制手段保存、汇编学位论 文。为保护学院和导师的知识产权，作者在学期间取得的研究成果属于内蒙古大学。作者 今后使用涉及在学期间主要研究内容或研究成果，须征得内蒙古大学就读期间导师的同意； 若用于发表论文，版权单位必须署名为内蒙古大学方可投稿或公开发表。 学位论文作者签名： 日期： z 警 名期 签师 教 导 指日 内蒙古大学硕士学位论文 基于g p u 平台的k l u 并行算法的研究：对角线块的l u 分解 摘要 随着电子科学技术的发展，为设计高质量的电路，电路模拟必不可少。在 电路模拟的过程中要涉及到稀疏线性方程组的求解。而随着电路矩阵规模的不 断增大，对电路矩阵的求解已经成为电路模拟过程的一个瓶颈。针对电路模拟 过程中产生的电路矩阵韵特点，通常均采用直接法来进行此类线性方程组的求 解。目前常用的求解器有s p a r s e l 3 、s u p e r l u , k l u ，其中以t i m o t h yd a v i s 教 授开发的k l u 最为有效率。 k l u 主要由预处理，首次l u 分解，再分解和回代求解这几个部分组成。 再分解部分是算法的重要组成部分。在一次电路模拟过程中，就是通过多次调 用该部分来完成对稀疏矩阵的数值l u 分解。因此，本文主要对这部分的算法 进行了并行研究和探索，并提出了基于g p u 平台的可行的并行算法。 k l u 算法在l u 分解过程中，采用的算法是基于高斯消去法的 g i l b e r t p e i e r l s 算法。我们通过研究串行算法和程序提出了两种不同的并行设 想，并在g p u 平台上设计和实现了四种不同的并行算法pl i e n 算法、pu l e n 算法、pn k 算法和ps t r e a m 算法。我们在实验平台i 上对这四种并行算法进 行了性能测试和分析，通过分析我们发现ps t r e a m 算法较前三种并行算法在性 能上有较大的优势，但由于该并行算法受限于实验平台i 中g p u 显存的限制， 导致并行度较低，性能较串行算法有所下降。为提高ps t r e a m 算法的并行度， 我们在g p u 显存容量更大的实验平台i i 上对其进行了性能的测试和分析。通 过分析，我们发现随着并行度的提高，ps t r e a m 并行算法性能也随之得到提升， 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 但还是受限于显存容量的限制，导致性能较串行算法并没得到提升。 由于我们是首次尝试对k l u 算法在g p u 平台上进行并行算法的实现，加 之稀疏矩阵数据的稀疏性、l u 分解数据的前后依赖性、硬件限制以及自身编 程经验的不足，导致并行算法性能较原有串行算法略有下降，但在本文中我们 提出的一些并行设想以及尝试也能为同方向的研究者提供很好的借鉴。 关键词：k l u 算法，数值分解，g p u ，并行算法 i r c u i t s i m u l a t i o ni st h en e c e s s a r yp a r tf o rd e s i g n i n gh i g l aq u a l i t yc i r c u i t i nt h ep r o c e s so f s i m u l a t i n gc i r c u i t s ，t h es o l v i n go ft h es p a r s el i n e a re q u a t i o n si si n v o l v e d w i t ht h e i n c r e a s i n g s c a l eo fc i r c u i tm a t r i x ，t h es o l v i n go fc i r c u i tm a t r i xh a sb e c o m ea 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 a n dp r o g r a m ，p u tf o r w a r dt w od i f f e r e n tp a r a l l e li d e a s ，a n dd e s i g na n dr e a l i z a t i o no f t h ef o u rd i f f e r e n tp a r a l l e l a l g o r i t h mp _ l l e na l g o r i t h m ，p u l e na l g o r i t h m ，p _ n k a l g o r i t h ma n dp s t r e a ma l g o r i t h mo nt h eg p up l a t f o r m w et e s t e da n da n a l y s e dt h e p e r f o r m a n c eo ft h ef o u rp a r a l l e la l g o r i t h m si ne x p e r i m e n t a lp l a t f o r mi t h i sp a p e r f o u n dt h a tp s t r e a ma l g o r i t h mh a sl a r g e ra d v a n t a g ei np e r f o r m a n c et h a np r e v i o u s t h r e ep a r a l l e l a l g o r i t hi np e r f o r m a n c eb ya n a l y z i n g b u tt h ep s t r e a ma l g o r i t h m r e s t r i c t e db yt h eg p um e m o r yl i m i t ，t h ep e r f o m a n c ei sl o w e rt h a ns e r i a la l g o r i t h m f o ri m p r o v i n gt h e p a r a l l e ld e g r e eo ft h ep _ s t r e a m ，t h i sp a p e rt e s t e da n da n a l y s e dt h e p e r f o r m a n c eo ft h ea l g o r i t h mi nt h ee x p e r i m e n t a lp l a t f o r mi i ，w h i c hh a st h el a r g e r g p um e m o r yc a p a c i t y t h r o u g ha n a l y s i s ，t h i sp a p e rf o u n dt h a tp s t r e a ma l g o r i t h m p e r f o r m a n c ea s c e n s i o nw i t ht h ep a r a l l e lr i s e ，b u tt h ep e r f o r m a n c ei s n o tp r o m o t e d s e r i a la l g o r i t h m ，s t i l ll i m i t e db yt h el i m i t a t i o no ft h eg p u m e m o r y s i n c ew ea r ef i r s ta t t e m p tt or e s e a r c ha n dr e a l i z et h ep a r a l l e lk l u a l g o r i t h m o n t h eg p u p l a t f o r m ，a n dt h es p a r s e n e s so fs p a r s em a t r i xd a t a ，d a t ad e p e n d e n c ea tl u f a c t o r i z a t i o n ，h a r d w a r ec o n s t r a i n t sa n dm yp o o rp r o g r a m m i n ge x p e r i e n c e ，r e s u l t i n g i np a r a l l e la l g o r i t h mp e r f o r m a n c el o w e rs l i g h t l yt h a no r i g i n a ls e r i a la l g o r i t h m b u t s o m ep a r a l l e li d e a sa n dt r y i n g ，w h i c hw ep r o p o s e di nt h i sp a p e r , c a na l s op r o v i d e v e r yg o o dr e f e r e n c ef o rt h es a m ed i r e c t i o nr e s e a r c h e r s k e y w o r d s ：k l u a l g o r i t h m ；n u m e r i cf a c t o r i z a t i o n ；g p u ；p a r a l l e la l g o r i t h m i v 内蒙古大学硕士学位论文 目录 第一章绪论1 1 1研究背景l 1 2 国内外研究现状j 2 1 3 研究目标与内容2 1 3 1 研究目标一2 1 3 2 研究内容o 一：。j 3 1 4 论文结构3 第二章g p u 通用计算与c u d a 简介4 2 1g p u 通用计算：一4 2 2c u d a 简介6 2 2 1c u d a 开发6 2 2 2c u d a 编程模型6 2 3 本章小结7 第三章k l u ：8 3 1k l u 简介8 3 2k l u 求解稀疏线性方程组流程：9 3 3k l u 在电路模拟中的应用1 l 3 3 1 电路矩阵的特点1 1 3 3 2 电路模拟中的线性系统一1 1 3 4 本章小结1 2 第四章k l u 之l u 分解1 3 4 1 稠密l u 分解一1 3 4 2 稀疏l u 分解1 4 4 3l e f tl o o k i n g 高斯消去法1 5 4 4g i l b e r t p e i e r l s 算法1 7 4 4 1 符号分析：1 7 4 4 2 数值分解1 9 4 5k l u 再分解机制：。1 9 4 6 本章小结2 0 第五章k l u 分解阶段的并行算法研究2 1 5 1l 江，ur e f a c t o r 算法_ 2 1 v 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 5 1 1 算法描述2 1 5 1 2 算法分析2 3 5 2pl l e n 并行算法2 4 5 2 1 算法描述j 2 4 5 2 1 并行性分析2 5 5 3pu l e n 并行算法2 5 5 3 1pu l e n 算法描述2 5 5 3 2 优化策略2 8 5 4pn k 并行算法_ 2 9 5 4 1p1 1 1 ( 算法描述2 9 5 4 2 并行性分析3 2 5 5ps t r e a m 并行算法3 2 5 5 1ps t r e a m 算法描述：。3 2 5 5 2 并行性分析k 3 4 5 6 本章小结3 6 第六章实验环境和实验方法3 7 6 1 实验环境。3 7 6 2 实验方法3 8 6 2 1 测试矩阵3 8 6 2 2 测试方法：3 9 第七章实验结果和性能分析4 0 7 1k l ur e f a c t o r 串并行算法性能比较4 0 7 2 并行算法性能分析：4 2 第八章全文总结与进一步工作。：4 5 8 1 全文总结。4 5 8 2 迸一步工作。4 6 参考文献4 7 致谢：4 9 v i 图2 1c p u 和g p u 不同的设计理念 4 图2 2g p u 体系结构j 5 图2 3c p u 与g p u 计算性能和存储器带宽5 图2 4c u d a 编程模型，6 图2 5c u d a 设备存储器模型：7 图3 1k l u 算法执行过程9 图3 2k l u 算法流程图。10 图4 1 求解l x = b 中石的非零模式图1 9 图5 1 矩阵块所有列的分解流程3 3 表6 1 平台i 硬件环境3 7 表6 2 平台i 软件环境3 7 表6 3 平台硬件环境3 8 表6 4 平台i i 软件环境3 8 ， 表6 5 实验用例3 9 表7 一实验用例详细信息4 0 表7 2k l ur e f a c t o r 与前三种并行算法在平台i 上的测试结果j 4 0 表7 3pl l e n 与pu l e n 算法不同阈值测试结果4 1 表7 4k l ur e f a c t o r 与ps t r e a m 并行算法的测试结果4 2 表7 5ps t r e a m 算法不同阈值测试结果4 2 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 内蒙古大学硕士学位论文 第一章绪论 1 1研究背景 随着电子科学技术的发展，为设计高质量的电路，电路模拟必不可少。s p i c e 模拟器【l 】 作为最为普遍的电路模拟程序，被广泛应用于v l s i 电路设计与模拟中。用s p i c e 模拟器 进行电路模拟时，需要求解大量结构相同，但数值不同的稀疏非对称线性方程组a x = b ，其 中a 是n x n 阶非奇异方阵，x 和b 是1 1 维列向量。 在许多实际问题中，矩阵a 的元素中含有大量的零元，而且随着矩阵阶数的增加，矩 阵的稀疏程度也越来越大。通常称这种零元占优的矩阵为稀疏矩阵，而相应的线性方程组 为稀疏线性方程组【2 翔。在对大型稀疏线性方程组进行求解时，一般不采用求逆矩阵的方法， 因为这种方法的复杂度较高，不利于硬件实现，因此通常用迭代法或直接法进行稀疏线性 方程组的求解【4 5 】。 迭代法的理论基础相对复杂，并且具有多种不同的具体算法，但其基本形式均为从一 个猜测解出发，把方程变形为x ( k ) = b x ( k - 1 + c 的形式( 上标表示迭代次数) ，通过多次迭代逐 步逼近方程组的精确解的方法【5 】。直接法由高斯消元法发展而来，对于一般矩阵可以使用 l u 分解【6 】，即a = l u ，( 其中l 为下三角矩阵，u 为上三角矩阵) 。分解之后，原来的方程 组可以表述为a x = l u x = b ，如果引入一个新矢量y = u x ，那么y 的数值可以通过求解l y = b 获得确定y 的数值后，x 的值可以用相似的方法求解y = u x 获得。在求解大型稀疏线性方程 组时，由于迭代法具有计算量小，近似解精度高，占用内存单元较少，设计程序简单的特 点而被广泛应用。但是，如果使用直接法时l u 分解过程能够被很多后续计算重复使用，则 后续的三角阵求解可以非常快速实现，此时直接法在性能上具有优势。由于在模拟电路瞬 态仿真中，需要求解大量结构相同，但数值不同的稀疏非对称线性方程组a x - - b 。因此s p i c e 模拟器中的稀疏线性方程组求解器均采用直接法来进行线性方程组的求解，目前常用的求 解器有s p a r s e l 3 、s u p e r l u 、k l u ，其中以t i m o t h yd a v i s 教授开发的k l u 最为有效率 7 , 8 , 9 , 1 0 , 1 1 】。针对电路矩阵的结构特征，k l u 算法对求解此类矩阵方程非常有效。然而，通常 来说，在电路模拟中产生的矩阵规模都比较大，而且相同模式的矩阵要进行上百次的求解。 虽然每次求解用时并不是很长，但整个求解过程在电路模拟过程中占用了相当可观的时间， 已经成为了电路模拟仿真的一个瓶颈，因此探索新的算法或通过并行思想来提高k l u 算法 的性能就显得很有必要。 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 n v i d i a 公司于2 0 0 7 年推出了g p u ( g r a p h i cp r o c e s s i n gu n i t ，图形处理器) 通用计算 产品c u d a ( 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 的推 出，g p u 高性能编程得到了迅速的崛起，现在被广泛应用于各个领域。凭借g p u 强大的计 算能力和存储器带宽，在g p u 上开发的并行程序相对c p u 上的串行程序，大部分都能达 到2 3 个数量级的加速。同时g p u 具有良好的可编程性，较低的成本和较短的开发周期， 已经成为了计算领域发展的新趋势【1 2 1 。因此，在g p u 上探索和研究并行的稀疏线性方程组 k l u 分解算法具有很好的现实意义。 1 2 国内外研究现状 k l u 软件包是t i m o t h yd a v i s 教授近年来新开发的求解稀疏线性方程组的求解器，因其 优良的性能和稳定的求解效果，目前被广泛应用于电路模拟器中。k l u 软件包中求解电路 矩阵的核心算法是k l u 算法，相比其他l u 算法，由于其出现时间并不是太长，所以国内 外对其并行算法的研究相对较少，仅实现了基于f p g a 的并行k l u 算法【l3 1 。目前使用g p u 实现大规模稀疏线性方程组直接法求解方面的工作几乎没有，因此研究基于g p u 平台的面 向大规模电路模拟的稀疏线性方程组求解的k l u 并行设计与实现，从研究角度存在着大量 创新的机会。 1 3 研究目标与内容 1 3 1 研究目标 k l u 是求解大型稀疏线性方程组的一个软件包，有效的应用在电路模拟矩阵的求解中。 然而，随着矩阵规模的不断增大，对电路矩阵的求解已经成为电路模拟全过程的主要瓶颈 之一，因此寻求高效的矩阵求解算法显得尤为重要。目前，解决这一问题主要有两个思路， 一种思路是寻求新的矩阵求解算法，另一种思路则是对已有的串行求解算法进行并行研究。 k l u 是目前最好的求解矩阵软件包，主要包括对矩阵的前期预处理，对角线块l u 分解再 分解，回代求解这些部分，且对角线块l u 分解是k l u 软件包中的主要部分。由于国内外 对k l u 的并行研究并不多，所以在本文中我们准备对其进行并行算法研究的尝试。由于时 间和精力的限制，本论文只对l u 分解部分在g p u 上进行并行的探索和研究。 2 内蒙古大学硕士学位论文 1 3 2 研究内容 全面理解k l u 算法，理清k l u 算法求解a x = b 的流程。 仔细研究k l u 算法中的l u 分解部分。 对k l u 算法中的l u 分解部分进行并行性研究。 阅读k l u 软件包中的源程序，从程序级再次挖掘其可并行的部分。 ， 编写基于g p u 平台的并行l u 分解程序。 测试程序的正确性与效率。 1 4 论文结构 本文按如下方式组织： 第一章介绍本文的研究背景和内容。 第二章介绍g p u ( g r a p h i c sp r o c e s s i n gu n i t ) 以及c u d a ( c o m p u t eu n i f i e dd e v i c e a r c h i t e c t u r e ) 的基本知识。 第三章对k l u 算法进行详细的介绍，明确其求解稀疏线性方程组的流程。 第四章对k l u 算法中的l u 分解阶段的算法进行分析研究。 第五章在g p u 平台上对k l u 中l u 分解中存在的并行性进行探索和研究，并提出了 可行的并行方案。 第六章对本文中测试程序所用的软硬件平台及测试用例进行了介绍。 第七章在g p u 平台上对并行算法进行测试并进行性能分析。 第八章对全文进行总结和讨论。 3 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 第二章g p u 通用计算与c u d a 简介 近年来，随着g p u 硬件技术的发展和高性能计算的兴起，g p u 通用计算技术得到飞速 发展。2 0 0 7 年，n v i d i a 推出了基于g p u 的c u d a 编程技术1 4 , 1 5 1 更加使得g p u 在高性能 计算领域得到重大的应用。c u d a 能够有效的利用g p u 强大的浮点运算能力，将g p u 通 用计算应用于数值计算、天文计算、流体力学计算等领域，并在这些应用中获得了一到两 个数量级的加速，取得了令人瞩目的成纠1 6 1 。本章我们将对n v i d i ag p u 的结构和c u d a 进行介绍。 2 1g p u 通用计算 目前，主流计算机中的处理器主要有c p u 和g p u 。传统上，g p u 只是负责图像渲染 任务，大部分的处理都交给c p u 。 2 l 世纪人类所面临的重大科技问题，如卫星成像数据的处理、全球气候准确预报、基 因工程等，这些问题的数据规模已经达到t b ( t e r ab y t e ，2 4 0 b y t e ) 甚至p b ( p e t ab y t e ， 2 5 0 b y t e ) 量级，没有万亿以上的计算能力是无法解决的。这对计算速度提出了严峻的挑战。 g p u 在处理能力和存储器带宽上相对c p u 有明显优势，在成本和功耗上也不需要付出 太大代价，从而为这些问题提供了新的解决方案。由于图形渲染的高度并行性，使得g p u 可以通过增加并行处理单元和存储器控制单元的方式提高处理能力和存储器带宽。如图2 1 所示，c p u 与g p u 采用了两种完全不同的设计理念，我们从该图可以看出，g p u 将更多 的晶体管用作执行单元，而不像c p u 那样用作复杂的控制逻辑和缓冲，因此，g p u 拥有很 强的并行计算能力。 固囹塑翌霪翌翌翟霎隧圈 圈囵口 】工至互互囡 囹 囫囹习】j 互夏歪互圈囫 豳囵 】豆重正e 匠囫 圈囹 l 互互工互e 匪囹 豳园j 工】互互蓝 日团囫 豳园 了】卫工互至e 园圃 c p ug p u 图2 1c p u 和g p u 不同的设计理念2 1 f i g u r e2 1c p u sa n dg p u sh a v ef u n d a m e n t a l l yd i f f e r e n td e s i g np h i l o s o p h i e s 【1 2 】 随着g p u 设计理念的不断发展，g p u 逐渐被用于越来越多的通用计算中。图2 2 是一 个现代的g p u 体系结构。在这种结构中，基本单元为流处理器组( s t r e a m i n gm u l t i p r o c e s s o r ， 4 内蒙古大学硕士学位论文 s m ) ，每个s m 包含8 个流处理器( s t r e a m i n gp r o c e s s o r ，s p ) ，每个s p 都有一个m a d ( m u l t i p l y - a d d ) 单元和附加的乘法单元，这些s p 共享控制逻辑和指令缓存。同时，我们 从图2 3 中也可以看出g p u 强大的计算能力，如g 8 0 系列的g p u 的浮点计算能力达到 5 0 0 g f l o p s ，g t 2 0 0 的浮点计算能力已经超过1 t f l o p s 。另外，现代的g p u 拥有强大的带宽， 如g 8 0 的存储器带宽可以达到8 6 4 g b s ，g t 2 0 0 的存储器带宽可以达到1 6 0 g b s 。以上的 这些g p u 体系结构的特点都决定了g p u 强大的并行计算能力。 “ 图2 2g p u 体系结构【1 6 】 f i g u r e2 2a r c h i t e c t u r eo fag p u t l 6 】 _ 糍 7 。一_ m 黼渊魄蚺h 懒 z 十撕蝴， _ - i - h 喇执棚 _ t ， 7 潲黼二| 絮 i 基多一 嘲m i 曩妇 错# 掰牺0 _ 例 甜一甜x 鳓 一钠 瓣一揪鼢， 蜮一懒溉， 一j 、夕 一= 雾阳 婶i - 掣叫品咖 硼删硒蹦硼黼j 晰翟酾口 图2 3c p u 与g p u 计算性能和存储器带宽n 们 f i g u r e2 3c o m p u t i n gp e r f o r m a n c ea n dm e m o r yb a n d w i d t hb e t w e e nc p u sa n dg p u s 1 6 1 | 。 5 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 2 2c u d a 简介 2 2 1c u d a 开发 2 0 0 7 年6 月，n v i d i a 推出了c u d a ( c o m p u t eu n i f i e dd e v i c ea r c l l i t e c m r e ，统一计算 设备架构) 1 4 , 1 5 , 1 7 】。c u d a 是一种将g p u 作为数据并行计算设备的软硬件体系。c u d a 不 需要借助于图形学a p i ，并采用了比较容易掌握的类c 语言进行开发。当然，要开发高性 能的g p u 通用计算程序，开发人员仍需掌握并行算法和g p u 架构方面的知识。 2 2 2c u d a 编程模型 c u d a 编程模型将c p u 作为主机，g p u 作为c p u 的协处理器 或称设备) 。在这个模 型中，c p u 与g p u 协同工作。c p u 负责逻辑性强的事务处理和串行计算，g p u 负责高度 并行化的计算任务。c u d a 采用了s i m t ( s i n g l ei n s t r u c t i o n ，m u l t i p l et h r e a d ，单指令多线 程) 【1 6 1 执行模型，s i m t 是对s i m d ( s i n g l ei n s t r u c t i o n ，m u l t i p l ed a t a ，单指令多数据) 【3 】 的一种改进。c u d a 编程模型把在设备上执行的并行计算函数叫做k e r n e l ( 内核函数) 。在 内核中最基本的单位是线程( t h r e a d ) ，多个线程组合成线程块( b l o c k ) ，多个线程块又组合 成网格( g r i d ) ，目前，每个k e r n e l 中只能包含一个鲥d 。在内核中，我们通过t h r e a l d x x ， t h r e a l d x y ，t h r e a l d x z 以及b k l d x x ，b l o c k l d x y 确定线程在总线程组织中的线程号i d 。一个 k e r n e l 函数并不是一个完整的程序，而是整个c u d a 程序中的一个可以被并行执行的步骤。 如图2 4 所示，一个完整的c u d a 程序时由一系列的设备端k e r n e l 和主机端得串行处理步 骤共同组成的。 行代码 井行内攘 l m e l ( ( ) ) ) 厂 i ，：_ j 丽焉孵舅臻嘲 8 l 钟k f o 。0 3i i戤呻“1 0 7i 瞰“2 ， 蔷嚣嚣蒋餐载嚣器磊蓠蔷辅蓉羲毳蔷毳 b i k t o i f b l t “，1 )l 叫o l k l 2 。1 ) 嚣祷嚣篱饕 l 嚣黼蓍蔷韬嚣篱嚣蓊鞯 图2 4c u d a 编程模型【1 刁 f i g u r e2 4c u d ap r o g r a m m i n gm o d e l t l 2 1 6 间。设备有自己的d r a m 空问， c u d a 的多层存储结构如图2 5 读写，从而实现了主机和设备之 间的数据传递。在g p u 存储模型中，主要有以下的几种存储器： 全局存储器，拥有很大的空间，但是访问速度较慢，每次访问有2 0 0 3 0 0 周期的延迟； 本地3 2 位寄存器，g 2 0 0g p u 上每个s m1 6 3 8 4 个寄存器； 共享存储器，每个b l o c k 共享，g 2 0 0g p u 上每个s m 拥有1 6 k b 的空间； 常量存储器，所有缸a d 只能对其进行读操作，拥有自己的c a c h e ，访问速度快。 图2 5c u d a 设备存储器模型【1 6 1 f i g u r e2 5c u d ad e v i c em e m o r ym o d e l 0 6 在c u d a 程序，真正的执行单位是线程束( w a r p ) ，w a r p 的大小由硬件的计算能力决 定。在g t 2 0 0b 8 1 g p u 中，一个w a r p 有3 2 个连续的线程组成。 c u d a 程序的编译器为n v o 圮，主机端代码由n v c c 编译成c 文件形式输出，而设备端代 码由n v c c 编译成p t x ( p a r a l l e lt h r e a de x e c u t i o n ) 代码【1 9 】或者二进制代码。 在c u d a 编程模型中，由于硬件资源的限制，c u d a 程序也存在众多的限制，我们以 t e s l ac 1 0 6 0 2 0 ，2 1 2 2 1 为例，每个b l o e k 最多包含5 1 2 个t h r e a d ，每个s m 最多同时执行1 0 2 4 个 线程，每个s m 最多8 个b l o c k ，每个s m 的最大寄存器数目为1 6 3 8 4 ，每个s m 的最大共 享存储器的大小为1 6 k b ，以上的这些限制因素决定了我们设计的并行算法的资源占有率， 并且我们只有综合考虑了上面的各个因素才能达到最好的性能。 2 3 本章小结 本章主要对g p u 通用计算的发展和c u d a 编程进行了介绍。通过介绍，我们了解到了 g p u 强大的计算能力和开发g p u 并行程序的c u d a 编程模型。 7 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 第三章k l u 3 1k l u 简介 k l u 1 0 , 1 1 1 ( c l a r kk e n tl u ) 软件包，是专门针对大型电路模拟矩阵设计的一种新型的稀 疏矩阵求解算法，它以l e f t 1 0 0 k i n g 高斯消去法为基础。所谓l e f t 1 0 0 k i n g 高斯消去法，即 对矩阵进行l u 分解时，从原矩阵的左起第- y u 开始，依次向右分解各列的一种l u 分解法。 k l u 软件包主要包括b t f , a m d ，c o l a m d ，k l u 等部分，其中b t f 、a m d 与c o l a m d 均用于线性方程组求解的预处理阶段，即对于原始方程组 出= b ( 3 1 ) 经过b t f 变换后系数矩阵变成了上三角块的形式，a m d 与c o l a m d 主要是针对经过 b t f 变换后得到的新系数矩阵对角线上的各个块进行变换，目的是减少l u 分解阶段的非 零元填充；k l u 是软件包的核心部分，是对对角线上的块进行l u 分解以及对x 的回代求 解。 相应的k l u 算法可分为四个阶段：预处理阶段、符号分析阶段、数值分解再分解阶段、 求解阶段。各阶段的详细说明如下： 在预处理阶段，将系数矩阵转化为上三角块的形式，也即对角线块下面的元素都为零 元( k l u 。_analyze) 在符合分析阶段，将每个对角线块进行符号分解，以确定每个对角线块分解以后l 、u 矩阵的非零模式( k l u _ f a e t o r ) 在数值分解再分解阶段，将每个对角块进行数值分解，确定每个对角块分解以后l 、u 矩阵的非零元素值( k l ur e f a c t o r ) 。 在回代求解阶段，首先求解右下角的对角线块，然后逐个回代，求出最终的解 ( k l u s o l v e ) 。 相应地，k _ l u 软件包也分为三个部分：k l u _ _ a n a l y z e 、k l u f a c t o r k l u r e f a c t o r 、k l u s o l v e 。 例如，对于线性方程组a x = b ，其中 8 内蒙古大学硕士学位论文 a = 20 o一1 3 o 203o 03604 00l42 00o01 则该线性方程组的系数矩阵的变化过程如图3 - 1 所示： 23o0 0 30406 0一l一32o o0loo 0420l? 腻訇暑 黼、钞 f ；：旷0 0 - 弛1 - 3 i v1 褂6o r 4 i 【-010一一141巳20 000 i1 一 l：i f | ：垦。了o - 了i 川- 3 雕6o ：4l 【0 10一一14000 0 j 童il ii l 蕊匝弼 数值分解阶段 - - - - - 卜 对每个块进行数 值l u 分解 。”表示非零元，这一阶段只确定非零元的位置 俸i褥-0-ifa2j渊fz叁_4j 团b 2 0 0 俸 一 。 i ： i li o 蜘l a 3 - jl 。3 j 图3 1k l u 算法执行过程 f i g u r e3 1c o u r s eo ft h ek l ua l g o r i t h m 3 2k l u 求解稀疏线性方程组流程 对于稀疏线性方程组 令r 为对角线阵，则a x = b 可写为 令p ，q 分别为初等行、列矩阵，则 彳j = b 胁= r b p r a q a r x = p r b 9 ( 3 - 2 ) ( 3 - 3 ) 基于g p u 平台的k l u 并行算法研究：对角线块的l u 分解 令p ，q 为初等矩阵p ，q 再进行初等变换得到的矩阵，则 p r a q q r z = p r b ( 3 4 ) 矩阵( p 翩q ) 包含对角线块和非对角线块元素两部分，其中对角线上块是需要进行l u 分解的部分，而非对角线块上的元素是不需要进行l u 分解的。 ( p 肌q ) = l u + f ，其中l u 为对角线上所有块分解后的集合，f 表示非对角线块上的 元素集合。方程( 3 - 4 ) 就变为： ( l u + f ) q 7 x = p r b ( 3 5 ) 此时，便可得到原方程组的解为t 石= q ( l u + f ) 。1 ( p r b ) ( 3 6 ) k l u 算法流程可用图3 - 2 表示： a 伫 对角线上各块厂 对角线上的各块 p “p a q q “ a m d 算法 ll 进行按比例缩放 为减少l u 分解中填充进行的前期处理 p p r a q i q l 原方程组a x = b 变为( p 舢q ) q ，x = p r b 厕线上畚块l u 硪 、进对角线上元素为呈 l l z ：q ( l u + f ) - l ( 腑) 图3 2k l u 算法流程图 f i g u r e3 2f l o wc h a r to ft h ek l ua l g o r i t h m 对于一个3 x 3 块形式的方程组，方程( 6 ) 可表示为 互2 厶： o剁驰 1 0 ， 工0 o 厶 。l 为程方的应相 内蒙古大学硕士学位论文 厶。u 。宰五+ e ：木五+ 互，木墨= 骂 厶：宰五+ 最。木置= 岛 厶，宰托= 马 ( 3 - 8 ) ( 3 - 9 ) ( 3 - 1 0 ) 在回代求解阶段，我们先求解方程( 3 - 1 0 ) 得到五，然后将墨代入到方程( 3 - 8 ) 与 一一 ( 3 - 9 ) 中消去五，接着求解方程( 3 - 9 ) ，得到五，然后再代入方程( 3 - 8 ) 消去五