（计算机应用技术专业论文）并行代数多重网格算法的优化及应用.pdf

燃剃 u n i v e r s i t yo f sc i e n c ea n d t e c h n o l o g yo fc h i n a ad i s s e r t a t i o nf o rm a s t e r sd e g r e e o p t i mi z a t iona nd a pp lic a t iono f a l g ebr aicmul t igr ida l g or i t hmf or p a r a ll e i a u t h o r sn a m e ：x i a o f e n gz h u s p e c i a l i t y ：c o m p u t e ra p p l i c a t i o nt e c h n o l o g y s u p e r v i s o r ： a s s o c i a t ep r o f l a n g f a n gd o n g f i n i s h e dt i m e ：j u n e10 m ，2 010 中国科学技术大学学位论文原创性声明 本人声明所呈交的学位论文，是本人在导师指导下进行研究工作所取得的成 果。除已特别加以标注和致谢的地方外，论文中不包含任何他人已经发表或撰写 过的研究成果。与我一同工作的同志对本研究所做的贡献均已在论文中作了明确 的说明。 作者签名： 擗 签字日期：丝l ! ：至 中国科学技术大学学位论文授权使用声明 作为申请学位的条件之一，学位论文著作权拥有者授权中国科学技术大学拥 有学位论文的部分使用权，即：学校有权按有关规定向国家有关部门或机构送交 论文的复印件和电子版，允许论文被查阅和借阅，可以将学位论文编入有关数据 库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。本人 提交的电子文档的内容和纸质论文的内容相一致。 保密的学位论文在解密后也遵守此规定。 签字日期：训。q 弓 导师签名：盔蜴导师签名：! 龟兰左 - 答字bg q ：礁b 皇：兰： 摘要 摘要 近年来，随着一些实际应用领域中大规模稀疏矩阵求解问题的推动，代数多 重网格a m g ( a l g e b r a i cm u l t i g r i da l g o r i t h m ) 算法及其并行化的研究成为了数值 计算领域的热点。 本文在原始a m g 算法和m p r s 算法( m u l t i p l ep h a s eg r i d c o a r s e n i n g a l g o r i t h m ) 的基础上，设计了一种优化的动态阈值r s 算法d v r s ( d y n a m i cv a l u e g r i d c o a r s e n i n ga l g o r i t h m ) 。在v i s u a ls t u d i o2 0 0 8 环境下进行了实验，结果表明， 本算法适用于更广泛的领域，与原有的并行a m g 算法相比，在并行效率上有着 不错的提升，主要体现在迭代次数和迭代时间上。本文主要工作如下： 1 总结了目前常用的代数多重网格算法及其并行化方法，介绍了这些算法 的原理，并分析了这些算法各自的优缺点。 2 在已有的并行a m g 算法( r s 3 算法和m p r s 算法) 的基础上，引入了 动态阈值的概念，它每一轮都会按照阈值变更系数的大小发生变化，用于对单个 网格点进行检测。由此提出了优化的d v r s 算法来处理并行a m g 问题。 3 对于不同类型的稀疏矩阵，以7 * 7 规模的特殊矩阵模拟，并用不同的初 始阈值和阈值变更系数加上不同的函数计算机制进行并行化实现，研究和验证对 于各种不同类型的特殊矩阵采用哪个范围的初始阈值和阈值变更系数更能发挥 d v r s 算法的优势。 4 研究了基于代数多重网格算法的并行化实现，在对二维浅水波的简化 n a v i e r - s t o k e s 方程离散和分析数据之后，用优化的d v r s 算法和传统算法进行实 验并分析结果。实验表明，优化的d v r s 算法在网格点规模比较大的时候，能 够发挥比传统k s 3 算法更好的迭代速度和加速比。 关键词：代数多重网格算法并行计算d v r s 算法r s 3 算法二维水波动画 动态阈值阈值变更系数 一 a b s t r a c t a b s t r a c t r e c e n ty e a r s ，t h en e e dt os o l v et h es o l u t i o no fl a r g e - s c a l es p a r s em a t r i xi ns o m e p r a c t i c a la p p l i c a t i o n sh a ss p a r k e dg r e a ti n t e r e s ti nt h er e s e a r c ho fa l g e b r a i cm u l t i g r i d a l g o r i t h m ( a m g ) f o rp a r a l l e l w ep r o p o s e do n ek i n do fo p t i m i z e dd y n a m i ct h r e s h o l d v a l u er sa l g o r i t h m ( d v r s ) b a s e do nt h ee x i s t i n gp a r a l l e la m ga l g o r i t h m ss u c h a so r i g i n a la m g a l g o r i t h ma n dm u l t i p l ep h a s eg r i d c o a r s e n i n ga l g o r i t h m u n d e rt h ee n v i r o n m e n to f v i s u a ls t u d i o2 0 0 8 t h ee x p e r i m e n tr e s u l t si n d i c a t et h a tt h en e wa l g o r i t h mi ss u i t a b l e f o rm o r ew i d e s p r e a dd o m a i n sa n di sa b l et oi m p r o v ee x t e n d i b i l i t yo ft h ea m g p a r a l l e lc o m p u t i n g t h em a i nr e s e a r c hc o n t e n ta n dt h ea c h i e v e m e n t sa r ea sf o l l o w 1 w r es u m m a r i z et h ep r e s e n tc o m m o n l yu s e da l g e b r a i cm u l t i g r i da l g o r i t h ma n d t h ep a r a l l e lc o m p u t i n ga l g o r i t h ma n da n a l y z et h ep r i n c i p l e so ft h e s ea l g o r i t h m sa r e i n t r o d u c e da n dr e s p e c t i v eg o o da n db a dp o i n t so ft h e s ea l g o r i t h m s 2 b a s e do nt h ee x i s t i n gp a r a l l e la m ga l g o r i t h m s ( r s 3a l g o r i t h ma n dm p r s a l g o r i t h m ) ，w ei n t r o d u c et h ec o n c e p to fd y n a m i ct h r e s h o l d v a l u e i tc h a n g e si t s s i z e d e f e rt ot h et h r e s h o l dv a l u ec o e f f i c i e n ta n di su s e dt oc a r d e do nt h ee x a m i n a t i o nt o e v e r ys i n g l eg r i dp o i n t t h u sw ei n t r o d u c et h eo p t i m i z e dd v r sa l g o r i t h mt o d e a l w i t ht h ep a r a l l e la m gi s s u e s 3 r e g a r d i n gt h ed i f f e r e n tt y p e so fl a r g e - s c a l es p a r s em a t r i x e s ，w eu s et h e 7 幸7 s c a l e s s p e c i a lm a t r i xt o s i m u l a t ea n du s ed i f f e r e n ti n i t i a lt h r e s h o l dv a l u ea n d t h r e s h o l dv a l u ec o e f f i c i e n tt or e a l i z ep a r a l l e lc o m p u t i n g w es t u d ya n dv a l i d a t ew h i c h s c a l eo fi n i t i a lt h r e s h o l dv a l u ea n dt h r e s h o l dv a l u ec o e f f i c i e n ti ss u i tt od i f f e r e n tt y p e o fs p e c i a lm a t r i xf o rd v r sa l g o r i t h m 4 w es t u d yt h ep a r a l l e lc o m p u t i n gb a s e do na l g e b r a i cm u l t i g r i da l g o r i t h ma n d u s et h es i m p l i f i e dn a v i e r - s t o k e se q u a t i o no ft w o - d i m e n s i o n a ls h a l l o w w a t e rw a v et o e x p l a i na n dv a l i d a t ew h a tw ea r es t u d y i n g t h ee x p e r i m e n tr e s u l ts h o w st h a tt h e o p t i m i z e dd v r sa l g o r i t h mc a n d ob e t t e rt h a nt h et r a d i t i o n a lr s 3a l g o r i t h ma t i t e r a t i o ns p e e da n da c c e l e r a t i o nr a t i ow h e nt h eg r i ds c a l ei sl a r g e k e yw o r d s ：a l g e b r a i cm u l t i g r i da l g o r i t h m ，p a r a l l e lc o m p u t i n g ，d v r s a l g o r i t h m ， r s 3a l g o r i t h m ，t w o - d i m e n s i o n a lw a v ea n i m a t i o n ， d y n a m i ct h r e s h o l dv a l u e ，t h r e s h o l dv a l u ec o e f f i c i e n t i i i 一 2 5 本章小结2 0 第三章动态阈值和d v r s 算法2 1 3 1 动态阂值思想的引入2 1 3 2 动态阈值算法的引入2 2 3 3 动态阈值算法( d v r s ) 2 5 3 4 本章小结2 8 第四章动态阈值算法中的关键变量和函数2 9 4 1 函数f ( s c o p e ，d y n a m i c ) 的计算选择机制2 9 4 2 大规模稀疏矩阵的特性差异。3 2 4 3 代表性矩阵中关键系数和函数机制的选择3 4 4 4 关键变量和函数计算机制选择分析4 0 4 5 本章小结。4 0 第五章基于d v r s 算法的水波模拟4 1 5 1 二维浅水波方程及其离散4 l v 目录 5 1 1 水波类和a m g 类的定义及作用4 1 5 2 二维浅水波方程的求解及水波模拟4 3 5 2 1 二维浅水波方程的求解步骤。4 4 5 2 2d v r s 算法求解方程的步骤4 4 5 2 3 水波模拟的图形绘制4 6 5 3 基于三种a m g 算法的二维浅水波方程求解4 7 5 4 各类算法实验结果分析5 0 5 5 本章小结5 2 第六章结束语5 3 6 1 总结5 3 6 2 工作展望5 3 参考文献5 5 致谢5 9 在读期间发表的学术论文与取得的研究成果6 1 v i 绪论 第一章绪论 流体力学与人类日常生活和生产事业密切相关。流体力学的基本方程 n a v i e r s t o k e s 方程可以精确地描述生产生活中流体现象的物理过程。但是流体 力学方程的高效求解一直是影响其应用的关键，以前多用直接法求解，规模稍大 的问题需要耗费几个小时甚至几天的时间，不能做到实时应用。现在的实时应用 多采用迭代法，可以在合理的时间内求得问题的近似解。但是一般迭代法收敛速 度受网格尺度影响较大，随着网格剖分的细化收敛速度减慢，因而难以高效求解 大型流体力学问题。 在计算机图形学中，要逼真的模拟水波运动状态是一个非常有挑战性的课 题，水波作为我们生活的这个世界里最普通的一种自然现象自然是随处可见，所 以水波模拟的课题在计算机应用领域，比如游戏，广告，影视等各个领域中有着 非常广泛的应用，它的研究既有理论意义，也有实际意义。 1 1 二维浅水波方程的离散 关于水波的模拟，目前大致可分为以下两类： 1 基于波的分析，通过直接构造参数曲面来表示水面 1 2 。但是，这种 让波的振幅等随着时间变化而产生的动画，效果并不是很好，因为让他们随着时 间而变化是很难设计的。因而此类方法只适合于模拟单一，静态的水体运动。 2 基于物理模型 3 】 4 】用流体力学方程来描述水流的运动。可以用二维的 n a v i e r - s t o k e s 方程建立模型，用计算流体力学中的数值分析工具来求解方程，最 后由数值来得到水流的形态，这样的水波是由方程的初始条件和边界条件自动产 生，并容易收到控制。 在计算流体动力学领域中，n a v i e r - s t o k e s 方程属于最常用物理模型。但是要 直接求解n a v i e r s t o k e s 方程非常的困难，而且效率很低，很难达到实时应用的 要求，所以可以将之简化n - 维空间，得到相对更容易求解的二维浅水波方程： 要+ 1 , l 罢+ v 罢+ g 銎= o + + 1 ，+ 2 = u a t a xa vu 瓠 跏a v挑a h 、 瓦+ “瓦+ v 瓦+ g o y = o ( 1 1 ) 钟舐却。 l l l ， i o h + 鱼u ( h 一6 ) 】+ 昙 v ( h 一6 ) 】：o o l o x o y 式中，玑1 ，分别为x 、y 方向的速度，h 是水位高度，b 是河底高度，h - b 是 绪论 水深。方程的简化带来了一定程度的细节损失，但是因为降低了一维会减少很多 的计算量，比原方程更具有实时使用价值，并且通过求解方程，可以直接得到水 面数据的高度场，避免了求解三维方程带来的流体运动边界跟踪和几何表示等问 题。 基于物理的流体动画主要分为两种： 1 欧拉法：从研究流体所占据的空间中各个固定点处的物理量入手，分析 流体内部空间中每个固定点上流体速度、压强、密度等参数随时间的变化，欧拉 法是基于网格的方法，保持了网格的规则性，只能使用小步长时间离散； 2 拉格朗日法：从分析流体各个微团的运动着手，即研究流体中某一指定 微团的速度、压强、密度等描述流体运动的参数随时间的变化，拉格朗日法是基 于粒子的方法，对稳定性限制较少，可以使用大步长时间离散。 采用一种介于欧拉法和拉格朗日法之间的算法，即半拉格朗日法【5 】。它结 合了欧拉法的规则性和拉格朗日法的稳定性，在网格划分很细的同时，也能采用 大步长的时间离散，因此能加速流体的造型速度。 先考虑一维浅水波方程，再将其扩展到二维情况，最后可得到二维情况下的 离散形式 h , 9 + a t 2 9 c 訾警十譬$ 础硼攀+ 学， m 2 ， = 瓦+ r ( 砭鱼生乏云妻r 二丝+ 砭：! 訾 。 一硝，( 氇二氇+ 竭 1 2 求解思路 本文对流体力学方程推导出的二维浅水波方程的求解进行了研究，针对此类 方程的特点一大规模性、不规则性和稀疏性，采用代数多重网格法作为求解所导 出的用于水波模拟的二维浅水波方程的迭代方法，它能递归地将细网格的问题放 到粗网格上求解，特点是收敛速度与网格尺度无关。 大规模性是指问题的规模一般都很大；稀疏性是指矩阵中非零元素的个数偏 少，零元素的个数很多；不规则性是指矩阵中非零元素的分布没有规律。 目前已有的求解大规模不规则稀疏矩阵方程的迭代方法，多用如j a c o b i 和 g a u s s - s e i d e l 等一般的迭代方法，计算量大计算周期长，他们在问题规模增大 时，求解的迭代速度和迭代次数明显增大，不符合实时应用的要求，所以选用求 2 绪论 解效率更高的代数多重网格算法来求解此类方程。 多重网格方法m g 【6 】 7 】是在2 0 世纪6 0 年代，由苏联的计算数学家 f e d o r a n k o 提出了m g 算法最初的思想，当时没有引起其他人的注意，7 0 年代多 重网格方法开始获得迅速的发展，近年来，它已经成为求解大规模科学工程计算 问题最有效的方法之一。在m g 算法思想基础上，优于传统的j a c o b i 迭代和 g a u s s 。s e i d e l 迭代法的一种新的大规模稀疏矩阵求解方法代数多重网格a m g ( a l g e b r a i cm u l t i g r i da l g o r i t h m ) 算法就在上世纪八十年代初，由a b r a n d t 和 s f m e c o m i c k 等人提了出来 8 】 9 】。 代数多重网格算法的大致思路粗略描述如下，具体算法参见2 1 小节： 1 构造一系列的嵌套网格层，并求解各层的虚拟网格q ”，插值算子职， 限制算子r “，网格算子么”和光滑算子g ”。在这里假设矩阵方程是n * n 规模 的矩阵，乘以n 宰1 规模的未知量向量，等于n * i 规模的结果向量。则未知量向 量中的n 个未知量就是n 个网格点，这n 个网格点组成了原始的网格；n * n 规 模的矩阵就是原始的网格算子： 1 1 粗化过程：在原始的网格中，求解各个网格点的影响值，找出具有最大 影响值的网格点设为粗网格点，根据原则把另一部分和粗网格点相关的网格点设 置为细网格点，然后继续寻找尚未划分的网格点中具有最大影响值的点设置为粗 网格点，如此重复直到所有网格点都被划分为粗网格点或者细网格点，则所有粗 网格点组成了粗网格； 1 2 生成插值算子：由求解的粗网格，按照固定计算方法求此层的插值算子； 1 3 求解限制算子，粗网格算子：把插值算子转置即为限制算子，然后由同 一层的限制算子木网格算子宰插值算子，得到这一层的粗网格算子，这个粗网格算 子即为下一层的网格算子，如此重复构造一系列嵌套的网格和每层对应的系数， 最粗网格层的规模一般设置为2 或者3 ，这是为了方便直接求解； 2 由构造的嵌套网格层和各层系数进行迭代求解； 2 1 给原始的未知量向量一个估计解，这个估计解与精确解之间有一个误差， 这就是原始误差： 2 2 从最原始的网格层开始计算，如果还未到能直接求解的最粗网格层，则 进行一次光滑之后，进入下一层网格层，这里光滑成为前光滑，是为了尽量减少 网格粗化带来的原始网格特性的损失； 2 3 在最粗网格层直接求解得到一个亏量，亏量是指新误差与原始误差之差， 由这个亏量加上原始误差得到一个新的误差； 2 4 进行一次后光滑，这里后光滑的作用与前光滑相同，然后把这个新的 误差插值到上一层，如此重复直到插值回到原始网格层，得到原始网格层的 一个新的误差，由新的误差得到一个新的估计解，它会比刚开始估计的解更接近 3 为了提高a m g 算法的并行效率和性能；在国外，研究a m g 并行算法的文献稍 微多一些，但和国内的一样，其研究的目的都是为了提高a m g 算法实现并行的 效果。 综合国内国外的已有的关于a m g 算法的文献 1 0 1 1 1 1 可以看出：首先，a m g 算法的效率提高近年来得到了人们的重视，尤其是在具有大规模稀疏特性的矩阵 方程求解方面有着重要的应用；其次，已有a m g 方法的主要思路多数是，先对 基础a m g 算法实现并行化，然后对算法中各个用到并行方法的阶段，如粗网格 的选取、插值算子的构造等阶段，进行算法改进以达到提高并行效率的目的。 并行a m g 算法的研究，目前的研究成果有： 1 c l e a r y ，a j ，f a l g o u t ，r d ，h e n s o n ，v e ，j o n e s ，j e 等人于19 9 8 年， 基于一个并行极大独力集，提出了一种完全并行的粗化算法【l o 】即c l j p 算法， 极大地促进了a m g 的并行计算研究： 2 徐小文，莫则尧，在2 0 0 5 年，提出了m p r s 算法 1 1 】，算法主要是各处 理机按r s 算法对所含网格点进行粗化，通过设定同步条件，然后在拟边界上进 行通信，更新自己边界上的粗点，再继续重复粗化的一种并行计算机制； 3 h e n s o n ，v e ，y a n g ，u m ，等人于2 0 0 2 年，提出了f a l g o u t 等各不相同， 性能效率也不同的算法【1 2 】。 由于a m g 算法的内在串行顺序决定了其并行化的程度并不高，目前只能在 粗化环节、限制算子的计算，对较多的矩阵矩阵、矩阵向量运算实现并行化。 目前已有的研究成果中，代数多重网格法的并行算法的优化，大多是从以下两个 方向考虑的：其一是通过对粗化环节算法的改变来达到改善并行效率的目的，其 优点是原始粗化算法已经得到实现，并行化的实现相对简单，但缺点是可改进的 细节相对较少；其二是通过修改插值算子的构造方法，来达到改善迭代过程收敛 4 绪论 性的目的，但目前插值算子的构造方法相对成熟，要实现更简单有效的改动比较 不易。今后的工作将主要在这两个方面进行研究和探索，实现并行a m g 算法效 率的提高。关于a m g 算法的已有研究成果，本文将在第二章中给出更多更详细 的说明和介绍。 r s 3 算法除了在对较多的矩阵矩阵、矩阵向量运算实现并行化之外，在网 格粗化阶段，先把原始网格中的网格点分配给各个处理器，可以用片分割也可以 用交替分割，然后各个处理器分别在范围内找出一个具有最大影响值的网格点设 置为粗网格点，也即是说r s 3 算法相当于把原始网格分成几个部分，在各个处 理器上分别串行实现a m g 算法，这样容易产生一个问题，各处理器找出的粗网 格点，它在整个网格中未必具有很大的影响值，如果把它作为粗网格点保留到粗 网格中，会导致原始网格的特性无法很好的继承到粗网格中，原始矩阵的特性损 失较大。 同样的，c l j p 算法是把各处理器中具有最大影响值的点放到一个并行极大 独立集中，到最后在这个独立集中的网格点都会被设置为粗网格点，它和r s 3 算法存在同样的问题，没有对挑选出的粗网格点在整个网格中的影响值大小进行 比较考虑。 m p r s 算法在这方面加入了一个固定的阈值，各处理器中最大影响值的网格 点只有在大于这个闽值的时候才会被选为粗网格点，否则就将同步等待。这样的 改进相对之前的两个算法已经有了改进，但是固定的阈值在有些时候并不好，首 先不同的问题肯定要用不同的固定阂值大小，这需要大量的实验和经验判断，如 果阈值大小不合理，很可能导致阈值太大使得只有一个处理器工作，其他处理器 始终在同步等待，也可能阈值太小，使得每次判断都通过，各处理器挑选的网格 点仍然都被选为粗网格点，变相的简化成了r s 3 算法。 必须注意的是有时候粗网格往往不能把原始网格中的网格特性很好的保留 下来，只因为一些被选的粗网格点，在整个网格点范围内并不具有很大的影响值， 本文优化之后的并行粗化算法将很好的处理这个问题。 1 4 内容概述及本文贡献 参看了已有关于并行a m g 算法的文献后，使我对a m g 算法本身和并行实 现a m g 算法的理解更为深入。本文实现了代数多重网格算法及其并行化，特别 是其中粗网格选取这个阶段的并行化，在已有的一些并行粗化算法的基础上，加 入了动态阈值的思想，优化了粗网格选取这个阶段的并行化实现过程，在一定程 度上提高了算法的整体求解效率： 本文主要完成了以下工作： 1 介绍说明了代数多重网格算法( a m g ) 的思想和详细算法实现流程，重 5 3 介绍了当今流体力学方程中广泛应用到的n a v i e r - s t o k e 方程，介绍了n s 方程的离散形式( 二维浅水波方程) ，分析水波类数据结构和水波模拟时图形绘 制方面的工作。 4 在已有的并行a m g 算法( r s 3 算法和m p r s 算法) 的基础上，引入了 动态阈值的概念，它每一轮都会按照阈值变更系数的大小发生变化，用于对单个 网格点进行检测。由此提出了优化的动态阈值r s 算法d v r s ( d y n a m i cv a l u e g r i d c o a r s e n i n ga l g o r i t h m ) 来处理并行a m g 问题。 5 用d v r s 算法对二维浅水波方程迭代求解，实现二维水波动画的模拟并 给出实验结果和效果截图。 d v r s 算法最主要的思想是为各处理器中具有最大影响值的网格点设置一 个动态的阂值，这个阈值会随着迭代次数的变化而变化，如果各处理器选出的粗 网格点太少，动态阈值会变小，使得下一次闽值判断时，更多的网格点能被选为 粗网格点。实验表明，优化的d v r s 算法在网格点规模比较大的时候，能够发 挥比传统r s 3 算法更好的迭代速度和加速比。 对于不同类型的稀疏矩阵，以固定规模的特殊矩阵模拟，并用不同的初始阂 值和阈值变更系数进行并行化实现，再配合不同的计算选择机制，使得动态阈值 可以变化更多，再考虑对于各种不同类型的特殊矩阵采用哪个范围的初始阈值和 阂值变更系数更能发挥d v r s 算法的优势。 本文主要贡献如下： 1 加入了动态阂值的思想，对已有的并行代数多重网格算法做了优化并详 细介绍了优化后的动态阈值r s 算法； 2 对动态阈值r s 算法中的初始阈值系数和阈值变更系数这两个关键变量 以及计算选择机制分析并实验，选择不同的关键变量和计算选择机制可以使动态 阈值r s 算法表现出不同的特性，对可能出现的一些组合以及它们在算法中不同 的作用做了说明； 3 用优化的动态阈值r s 算法实现水波模拟，详细说明了实验系统的中计 算、绘图等设置情况，对实验结果做了分析。 1 5 论文组织 6 绪论 全文共分六章，章节安排如下： 第一章，对二维浅水波方程进行离散和数据分析，介绍了a m g 算法的历史， 研究现状和已有论文对本文的启发提示作用。 第二章，介绍了a m g 算法的实现原理，主要是介绍粗网格的选取和插值算 子的构造，并对内部数据进行分析，介绍了已有的几个主要的a m g 并行算法， 详细阐述了这几种算法( 包括内在串行r s 3 算法、完全并行c l j p 算法以及写下 本文灵感由来的m p r s 算法) 的原理及实现，分析了各个算法的优缺点。 第三章，介绍了动态阂值思想的由来，该思想在现实生活中适用性的分析以 及基于动态阈值思想的d v r s 算法，详细阐述了d v r s 算法的原理及实现。 第四章，总结了以上两章中所述算法的实现结果，由实验结果分析和验证各 个算法的优劣和不同的适用范围，再配以二维浅水波方程的实现结果，效果截图 来绘制曲线图并对其进行分析。 第五章，介绍了水波模拟的图形绘制工作，分析并用实验实现d v r s 算法 中两个最关键变量( 初始阂值和阈值变更系数) 在具有不同特殊性的矩阵( 如上 三角矩阵，对称矩阵等) 中不同的最佳取值范围，给出了程序实现的流程图以及 实验结果。 第六章，对本文所做工作进行了总结，并指出了今后的研究方向。 7 绪论 8 并行a i d g 算法 第二章并行a m g 算法 现今大规模稀疏矩阵都具有一些特性：边界越来越复杂，要求高精度。每个 网格块周围临近的网格数目不同，在此基础上离散得到的稀疏矩阵具有不规则的 特性。代数多重网格法是一种求解大型不规则稀疏矩阵方程的迭代方法。它递归 地将细网格的问题放到粗网格上求解，特点是收敛速度与网格尺度无关。 多重网格方法m o ( m u l t i g f i da l g o r i t h m ) l 约最大特点是，当网格剖分更精细时， 收敛速度不会变慢。代数多重网格法a m g ( a l g e b r a i cm u l t i g r i da l g o r i t h m ) 是 利用几何多重网格方法的思想建立起来的一种自动求解矩阵方程的迭代方法，与 几何多重网格法依赖求解问题的几何结构不同，它只需要代数方程组系数矩阵的 信息，因而具有很好的通用性，可以求解许多不同类型的应用问题。 代数多重网格a m g 算法是一种求解大规模不规则稀疏矩阵方程的迭代方 法，是数值计算领域中一种重要的加速迭代收敛技术。a m g 算法由建立和求解 两个阶段构成，可以递归地将细网格的问题放到粗网格上求解，除了建立阶段的 网格粗化部分以外，其余的大多数计算都属于矩阵与矩阵或矩阵与向量之间的计 算，可以方便地实现并行，难点在于网格粗化过程的并行化实现。目前相对成熟 和应用较多的并行a m g 算法主要有以下几种：r s 3 算法、c l j p 算法、m p r s 算法等。 2 1 代数多重网格( a m g ) 算法 在研究a m g 并行化的问题之前，先来看一下串行的a m g 算法的详细原理 及实现流程。 代数多重网格法( a m g ) 1 3 i 是利用几何多重网格方法的思想建立起来的一 种自动求解矩阵方程的迭代方法，它不利用所求解问题的几何和物理性质，仅仅 利用代数方程中的系数矩阵构造多重网格算法所需的五个分量，光滑误差的过程 基本是固定的，而且能自动给出一套虚构的粗细网格，然后按照通常的多重网格 循环过程进行求解。由于a m g 采用一种纯代数的层次方法直接处理矩阵方程， 因而对许多不同类型的应用问题都能表现出很好的算法稳健性，特别是能有效求 解二维或三维大规模( 稀疏) 无结构方程组。 一个代数多重网格算法主要分为两个过程： 1 预备过程：这个过程从系数矩阵构造出多重网格的五个分量：虚拟网格 q ”，插值算子，二，限制算子露”，网格算子彳”和光滑算子g 卅，这是a m g 算法中前期构造的主要部分： 9 并行a m g 算法 标准的多重网格循环迭代求解过程。 面讨论具体的代数多重网格方法，重点考虑如何定义和构造a m g 的五个 先考虑线性代数方程组 a x = b ，_ - - - b ，( i = i ，行) ( 2 1 ) 霉l 其中a = ( ) ，x = ( 五，) 7 ，b = ( b l ，也) 7 。构造一串越来越小的代数 方程组 j 4 ”x ”= b 册，瞄彳= ( f - 1 ，) ( 2 2 ) j = i 其中a ”= ( ”) 舢，彳”= ( “l ”，u n ”) 7 ，b ”= ( 6 l ”，也”) 7 ，m = l ，m ， 玎= 门川a 1 = a ，u 1 = u ，f 1 = f 。这一串方程组在形式上具有几何多 重网格方法中粗网格的作用，特别是一个网格点为与已知矩阵有关的有向图结点 的虚构网格q ”，它能够简单地看作是未知数x ? ( 1 ，) 的集合，我们可以形 象的把未知数当成是一个几何图形中的结点，他们之间的关系则为结点之间的有 向边，这样的图就是一个网格，所谓的粗网格就是在这个图中挑选一些有向边连 接较多的结点保留下来，连接他们之间的有向边需要重新来画，这样就挑选出了 粗网格。此外，还必须定义插值算子1 l 。：g ( q ”1 ) 寸g ( q ”) ，限制算子 名“：g ( f 2 ) 一g ( q ”1 ) 和松弛的方法g ” g ( f 2 ”) 专g ( f 2 ”) ，其中g ( f 2 ”) 是q ”上 网格函数的线性空间。只要这些量给出，就能够利用通常的多重网格法的循环过 程来求解代数方程组( 2 1 ) 。一般地，采用g a l e r k i n 原则来求限制算子和粗网格算 子a 肿1 ：g ( f 2 ”1 ) 专g ( f 2 册+ 1 ) 为 名= ( 职，) 7 ，a 卅+ 1 = 名“a ”i l 。 ( 2 3 ) 这样能保证粗网格算子保持细网格算子的许多有用的性质，如对称正定等， 且粗网格方程精确解的校正是插值算子值域里好的近似，即粗网格校正算子满足 变分原理，而松弛过程在代数多重网格法中是固定的，通常取g a u s s - s e i d e l 型的 松弛法，这样能够看出，只需对粗网格的选取和插值算子的构造给出定义，就可 以形成一个a m g 循环过程。 2 1 1 粗网格的选取 首先构造一串网格q 1 ，q m ，其中最细网格q 1 必须选得足够细能够保 证所需的精度，而最粗网格q m 的选取必须使得其上问题精确求解的花费，相比 最细网格上的一个松弛步的花费可以忽略不计。在几何多重网格方法中，一般地 q ”取为q 辨+ 1 的一致细化( 把多重网格方法在几何空间画出，q “1 就是在q ”中 1 0 并行a m g 算法 去掉一些结点之后的新的图) ，而在代数多重网格方法中，将网格q ”分为两部分 c 册和f ”( c ”即为粗网格，f ”则为被去掉的那些结点) ，更粗的网格q ”1 = c ”， f 肿- - q 册一c 脚，虚拟网格q 肘的划分只需利用代数方程组的信息。 定义q ”中个点f 是强连接到点j ，如果 ( 一硝) o o i 警a x ( 一) ，0 o o 1 ，f j ( 2 4 ) 这个定义实际上是针对m 阵的，即矩阵对称正定，且非对角元素非正，通 常取o o = 0 2 5 。 令掣是i 点的所有强连接点之集，e = c ”i 卵，划分q ”为集合c 和，” 的原则是： ( c 9 1 ) 对每个点f f ”，任何点f 应该或者是在掣中，或者至少强连接 到一个口中的点； ( c r 2 ) c ”应该是具有这种性质的所有点的最大子集，即c ”中没有两个c 点 是相互强连接的。 实际上，对许多方程组来说，同时严格满足这两个原则是不可能的，所以要 求( c r l ) 严格地成立，( c p 2 ) 作为一种构造粗网格的指导，利用这两个准则就可以 构造出一种选择c ”和f ”的具体算法，即可以构造粗网格了。 现在定义所有强连接到点f 的点之集为$ = j ：f s 7 ) ，且对一个集合尸， 假定l 尸| 定义为集p 中元素个数，j r u g e 和k s t u e b e n 给出了如下过程，说明影响 集、依赖集的定义和求法。首先，c 点的选择可如下进行： 1 设c “= ，f ”= ，u ”= q ”，且对所有f ，五= $ ) ： 2 选耿当前网格) 中有着最大丑值的f ，设c ”= c ”y i ) ，u = u - i ) ； 3 对所有j s t , iu 执行第4 步和第5 步； 4 设f ”= f ”y j ) ，且u = u - i ) ； 5 对所有f s j i 【，令丑= 丑+ 1 ； 6 对所有j s 7iu ，令兄，= 名，- 1 ： 7 若u = 矽，则停止；否则，回到第2 步。 上文试图通过在网格上均匀分布c 点以满足准n ( c r 2 ) ，然后结合插值权重 的计算，将生成的测试f 点用准贝l j ( c r l ) 检验，以确保( c r l ) 成立，当然如必需时 可增加新的c 点，应该注意到如果上述1 7 步被有效的执行，仅需o ( n ) 次运算 即可。 2 1 2 插值算子的构造 当粗网格给定后，就可以定义插值算子，设误差p 册= “? 一u ”，甜? 为离散方 1 1 并行a m g 算法 程的精确解，1 , 1 ”为离散方程的估计解，每个c ”中的变量也同时在粗网格中，这 样这些点的值直接有q 胛+ 1 中相应的量以权重1 给出，所以插值公式的关键实际 上是如何用q 卅+ 1 = c ”中的变量来计算f 聊中的变量，即每个点i f ”的值是由一 个较小的插值变量之集e ”互c 肼以e ”= 职。e 肿1 的形式给出，即 纠1 ) j _ 埘。善 ( 2 5 ) 很明显，为了计算效率的缘故，只肿应该取为点i 附近的c 点的集合，它必 须是c 朋的一个合理小的子集，当然任何这样的插值应为满秩的。 取只”= 口，得到一个针对m 阵的插值公式如下 = e c ，硝( o ) ，v i e f ” ( 2 6 ) 这里 哪) - _ 忐b + 叫醒赡。印鲻) ( 2 7 ) 其中w = j q ”，0 ) 表示点i 的所有邻点，秽= 坍一口， d ：= d ? is ? ，d ? = d ? 一d ： 2 1 3 光滑误差 考虑下面的p o i s s o n 方程，不考虑边界条件 一一2 厂( 工，y ) f 2 8 ) 其阻尼j a c o b i ( 雅克比) 迭代公式为 z + 1 ( x ，y ) 2 百1 【办2 厂( x ，j ，) + 甜( x 一矗，y ) f 2 9 ) + “七( x + ，y ) + “七( 工，y 一| i z ) + “七( x ，y + 矗) 】 “”= ( 1 一e a ) u + c oz + 1 若阻尼系数取；，得到 “+ 1 ( x ，y ) 2 了1 矗2 厂( z ，y ) + “+ 材( x 一 ，y ) ( 2 i o ) + “( x + ，y ) + “( x ，y - h ) + u ( 工，y + 向) 】 误差，= u 一, t ，贝l j v k + l 2 孝【，( x ，y ) + v ( x 一厅，y ) ( 2 1 1 ) + v ( x + ，y ) + v 。( x ，j ，一j z ) + 1 ，( x ，y + j i z ) 】 可以看到，第胁1 次迭代后每点的误差为第k 次迭代后该点及其周围4 点的 1 2 并行a k l g 算法 误差的平均。( 误差平均之后插值的效果就比较好) 许多经典的迭代方法( j a c o b i ，g a u s s s e i d e l 等等) ，如果被适当地应用到离散椭圆 问题，都有一个好的光滑误差的效果。 综上所述，可以让我们对代数多重网格( a m g ) 算法已经有了一个大概的印 象，下面详细看一下a m ( 3 的具体算法。 考虑线性代数方程组： 彳u = f ( 2 1 2 ) 其中，彳= ( 吩) 。( f = l ，以) ，u = ( ，“。y ，